MODELING OF CHARGE INJECTION AND TRANSPORT IN
ORGANIC SEMICONDUCTORS WITH APPLICATIONS TO

CONDUCTING ATOMIC FORCE MICROSCOPY

By

Kanokkorn Pimcharoen

A DISSERTATION

Submitted to

Michigan State University

in partial fulﬁllment of the requirements

for the degree of

Physics — Doctor of Philosophy

2018

ABSTRACT

MODELING OF CHARGE INJECTION AND TRANSPORT IN ORGANIC

SEMICONDUCTORS WITH APPLICATIONS TO CONDUCTING ATOMIC FORCE

MICROSCOPY

By

Kanokkorn Pimcharoen

Charge injection and transport in organic semiconductors are key factors controlling the

device performance, and have been intensively investigated by conductive atomic force mi-

croscope (c-AFM) experiments in the space-charge-limited current (SCLC) regime. The

simpliﬁed SCLC theory, despite being widely used to describe the unipolar SCLC, has lim-

itations in explaining the current-voltage responses of c-AFM measurements due to two

major reasons. First, the conventional planar model does not include the eﬀect of current

spreading commonly found beneath the conducting tip. Secondly, the theory only considers

drift transport, and assumes that charge diﬀusion can be neglected, causing discrepancies

in its predictions of transport behaviors that will be discussed thoroughly here. The focus

of this thesis is on developing numerical models for hole-only devices with the full descrip-

tion of drift and diﬀusion transport mechanisms, which is called the drift-diﬀusion (DD-)

SCLC model. The applications of the models in the analysis of c-AFM experimental data

are presented.

We generalize the theory which takes both drift and diﬀusion currents into account, lead-

ing to more realistic DD-SCLC models for several applications. We then develop numerical

approaches that eﬃciently simulate the hole-only SCLCs for one-, two-, and three- dimen-

sional systems. In the case of fully 3-D calculations, the DD-SCLC model is able to treat

inhomogeneous systems including spatially varying trap distributions, nanoscale morpholo-

gies, and the tip-plane (c-AFM) geometry. In the theoretical studies, the device simulations

elucidate a number of crucial factors that aﬀect the charge transport in the SCLC regime, in-

cluding charge diﬀusion, traps, as well as, nanoscale morphology. We introduce the method-

ology of characterizing the current-voltage responses from c-AFM measurements, which has

been used in elucidating the experiments on semiconductor poly(3-hexylthiophene) (P3HT)

thin ﬁlms that develop ﬁbrous morphologies after thermal annealing.

Vita brevis,
ars longa,

occasio praeceps,

experimentum periculosum,

iudicium diﬃcile.

— Hippocrates; c.460-c.370 BC

iv

ACKNOWLEDGMENTS

My PhD journey was nothing but one eventful ride that I was so fortunate to experience.

For that I would like to express my gratitude to all the people who taught me how science

worked in both simple and elegant ways.

First and foremost, I would like to express my deepest gratitude to my advisor, Professor

Phillip M. Duxbury, for guiding me through my graduate program with his insight, support

and great patience. Your broad spectrum from an awesome physicist, an excellent sportsman

as well as an article (a/an/the) modulator, will lesson me for life. Also, I would like to extend

my gratitude to the rest of my committee members, Professor Andrew J. Christlieb, Professor

Richard R. Lunt, Professor Vladimir G. Zelevinsky and Professor Pengpeng Zhang. It was

valuable to learn from you all.

It will always be my honor to have been a part of Professor Duxbury’s group. I would

like to express my gratitude to the former and the current members. To Jiwu, Dan, Luke,

Saurabh, Jenni, Corey, Non, Bin, Xukun, Tim and Brandon, thank you for all productive

discussions and amazing friendships. I cannot say enough thanks to Dan for simulating plenty

of model morphologies I can play with, and to Tim for proofreading my ﬁnal draft. The

thanks are also due to Jibing and Sean from Professor Zhang’s group for their contributions

in the experiments of semiconducting poly(thiophene).

Last but not least, I could not have ended my PhD years without my friends and family.

Thanks all my friends for being there for me through my ups and downs. My dissertation

would not have been possible without the very ﬁrst proofread from Jiwu. I do not think I

would make it out in one piece. It has been the best time of my life because you were there

to keep me sane, to remind me to live to the fullest, to toughen up, to have faith, to laugh,

v

to cry, and much more. I saved the largest thanks for last. To our big family in China and

Thailand, your invincible loves and ultimate supports meant the world to me.

I acknowledge the ﬁnancial support during my ﬁrst ﬁve years at MSU from the Thai

government scholarship.

vi

TABLE OF CONTENTS

LIST OF TABLES . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

LIST OF FIGURES . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

ix

xi

KEY TO SYMBOLS . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xx

Chapter 1

Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
1.1 Motivation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
1.2 Overview . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

Chapter 2 Theory of Charge Carrier Injection - Background . . . . . . . .
2.1 Charge Injection in Organic Semiconductors . . . . . . . . . . . . . . . . . .
2.2 Fundamental Equations to Describe SCLC in Hole Injection System . . . . .
2.3 The Simpliﬁed Theory of SCLC: Planar (Plane-Parallel) Geometry . . . . .
2.3.1 Trap-Free SCLC . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
2.3.2
. . . . . . .
2.4 The Simpliﬁed Model of SCLC: Tip-Plane (c-AFM) Geometry . . . . . . . .

SCLC with The Exponential Distribution of Trap States

Scaling Factors

Chapter 3 The Drift-Diﬀusion SCLC Transport Models - New Method .
3.1 Modeling the Hole Injection System . . . . . . . . . . . . . . . . . . . . . . .
3.1.1 Governing Equations . . . . . . . . . . . . . . . . . . . . . . . . . . .
3.1.2 Electrode Boundary Conditions . . . . . . . . . . . . . . . . . . . . .
3.1.3
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
3.2 Numerical Method . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
3.2.1 Numerical Derivatives
. . . . . . . . . . . . . . . . . . . . . . . . . .
3.2.2 Numerical Integrations . . . . . . . . . . . . . . . . . . . . . . . . . .
3.2.3
Scharfetter-Gummel Discretization . . . . . . . . . . . . . . . . . . .
3.2.4 Numerical Calculation Scheme . . . . . . . . . . . . . . . . . . . . . .
1-D DD-SCLC Model: Homogeneous Planar Geometry . . . . . . . . . . . .
3.3.1 Discretization . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
3.3.2 Potential Solver . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
3.3.3 Hole Density Solver . . . . . . . . . . . . . . . . . . . . . . . . . . . .
3.3.4 Veriﬁcation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
2-D DD-SCLC Model: Semi-Homogeneous Tip-Plane Geometry . . . . . . .
3.4.1 Discretization . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
3.4.2 Potential Solver . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
3.4.3 Hole Density Solver . . . . . . . . . . . . . . . . . . . . . . . . . . . .
3.4.4 Veriﬁcation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
3.5 The Fully 3-D DD-SCLC Model in a Cartesian Coordinate System . . . . . .
3.5.1 Discretization . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

3.3

3.4

1
1
4

6
6
8
10
12
14
18

22
22
23
26
27
28
29
32
34
37
40
41
42
43
44
50
52
55
59
63
72
75

vii

3.5.2 Potential Solver . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
3.5.3 Hole Density Solver . . . . . . . . . . . . . . . . . . . . . . . . . . . .
3.5.4 Veriﬁcation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
3.6 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

77
80
84
89

Chapter 4 The Drift-Diﬀusion SCLC Transport Models - New Results .
4.1 Diﬀusion Eﬀects . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
4.1.1 Planar Geometry . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
4.1.2 Tip-Plane Geometry . . . . . . . . . . . . . . . . . . . . . . . . . . .

91
91
96
99
4.2 Trap Eﬀects . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 103
4.2.1 Planar Geometry . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 104
4.2.2 Tip-Plane Geometry . . . . . . . . . . . . . . . . . . . . . . . . . . . 110
4.3 Morphological Eﬀects: Semi-Crystaline P3HT Thin Films . . . . . . . . . . . 116
4.3.1 An Isolated Fiber: ‘Face-On’ Fiber and ‘Edge-On’ Fiber
. . . . . . . 119
4.3.2 An ‘Edge-On’ Fiber Network . . . . . . . . . . . . . . . . . . . . . . 132
4.4 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 136

Chapter 5 The Drift-Diﬀusion SCLC Transport Models - Application to

an Experimental System . . . . . . . . . . . . . . . . . . . . . . . . 140
5.1 Experiments and Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 141
5.2 Theoretical Simulations and Analyses . . . . . . . . . . . . . . . . . . . . . . 144
5.2.1 Trap and Mobility Anisotropic Eﬀects
. . . . . . . . . . . . . . . . . 147
5.2.2 Traps and Morphological Eﬀects . . . . . . . . . . . . . . . . . . . . . 156
5.3 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 160

Chapter 6 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 162

APPENDIX . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 168

BIBLIOGRAPHY . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 175

viii

LIST OF TABLES

Table 2.1: Unipolar space-charge-limited current (SCLC) for hole injection of a planar
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

device.

21

Table 4.1: Summary of parameters used in the device simulations . . . . . . . . . . .

92

Table 5.1: Summary of parameters used in the device simulations to analyze c-AFM

data . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 147

Table 5.2: The eﬀective mobilities in the vertical direction found from ﬁtting a repre-
sentative examples of c-AFM current-voltage responses measured from an
annealed P3HT sample for a range of anisotropic ratios.

. . . . . . . . . . 155

Table 5.3: An analysis of point-to-point variation in transport properties across the
ﬁlm is presented for c-AFM data taken at 48 locations on the surface of
an annealed sample. The average, standard deviation, minimum and max-
imum of the eﬀective mobility in the vertical direction found from ﬁtting
the experimental data for a range of anisotropic ratio.

. . . . . . . . . . . 155

Table 5.4: The eﬀective mobilities in the vertical direction found from ﬁtting a rep-
resentative example of c-AFM current-voltage responses measured from a
non-annealed P3HT sample . . . . . . . . . . . . . . . . . . . . . . . . . . 155

Table A.1: Summary of parameters m and kJ estimated by ﬁtting theoretical current-
voltage characteristics of planar devices in the SCLC regime (see section
4.2.1 for details of device simulations) to the semi-empirical expression (4.3)
for a range of trap parameters ft and rt. Note that R-square of regression
is high (R2 > 0.9999).

. . . . . . . . . . . . . . . . . . . . . . . . . . . . 170

Table A.2: Summary of parameters m and kJ estimated by ﬁtting theoretical current-
voltage characteristics in the tip-plane geometry in the SCLC regime (see
section 4.2.2 for details of device simulations) to the semi-empirical expres-
sion (4.4) for a range of trap parameters ft and rt and a range of mobility
anisotropy µρ/µz. Note that R2 > 0.9999 for all regressions. Table contin-
ued on the next page.

. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 171

Table A.2: (cont’d) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 172

ix

Table A.3: Summary of parameters m and kJ estimated by ﬁtting theoretical current-
voltage characteristics in the tip-plane geometry in the SCLC regime (see
section 5.2.1 for details of device simulations) to the semi-empirical expres-
sion (4.4) for a range of trap parameters ft and rt and a range of mobility
anisotropy µρ/µz. Note that R2 > 0.9999 for all regressions. Table contin-
ued on the next page.

. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 173

Table A.3: (cont’d) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 174

x

LIST OF FIGURES

Figure 3.1: The schematic representation of hole-only devices The path of hole injec-
tion and hole transport are indicated . . . . . . . . . . . . . . . . . . . .

23

Figure 3.2: Numerical method:

at the mid point xi, (c) numerical integration (cid:82) xi+1

(a) illustration of spatial discretization in a one-
dimensional Cartesian coordinate system,(b) numerical derivative df /dx
f (x)dx around the
mid point xi+1/2 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

xi

Figure 3.3: A ﬂow diagram used to calculate the current-voltage characteristic of hole-
only devices. The algorithm is based on the self-consistent scheme. We
ﬁrst establish the device setup with its morphology and parameter values,
and then generate initial trial values for electric potential ψo and free
hole carrier po
f . Subsequently, Poisson’s equation is solved numerically to
ﬁnd the estimated values of potential ψn. Utilizing the new trial values
of potential ψ, the estimated values of free hole carrier density pn
f are
evaluated from the continuity equation and the drift-diﬀusion equation,
leading to the new trial values pf . This entire process is repeated until
convergence is reached, when the relative correction of free hole carrier
density is less than the value of the tolerance δ.
. . . . . . . . . . . . . .

Figure 3.4: Schematic illustration of a planar hole-only device made of a thin-ﬁlm
of organic semiconductor. The arrows illustrate the hole current path
responding to hole injection at the top electrode. When the material
is homogeneous, the vertical conductance is uniform. A 1-D DD-SCLC
model is often suﬃcient to describe the hole carrier transport in such a
device. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

Figure 3.5:

Illustration of the spatial discretization in 1-dimensional Cartesian coor-
dinate system. The device thickness L is discretized into the grid of N + 2
mesh points. Each mesh point is equally separated by ∆x = L/(N + 1).

Figure 3.6: Mesh size Veriﬁcation – Convergence of the simulated hole current den-
sity with respect to the number of mesh points, when the applied voltage
Va = 0.1, 1, 10V . Hole current density is calculated from the 1-D DD-
SCLC model of the trap-free (ft = 0) hole-only devices of 80nm thickness
(Nv = 1.25 × 10−21cm−3, µp = 10−4cm2V −1s−1). Discretizing the sys-
tem into 41 mesh points with the constant 2nm mesh interval provides a
well converged solution, and results in a maximum 7% error compared to
the result from 801 mesh points (0.1nm mesh interval). . . . . . . . . . .

xi

29

39

40

41

47

Figure 3.7: Nv Veriﬁcation – Convergence of the simulated hole current density with
respect to the eﬀective density of states in HOMO (Nv), when the applied
voltage Va = 0.1, 1, 10V . Hole current density is calculated from the 1-
D DD-SCLC model of the trap-free (ft = 0) hole-only devices of 80nm
thickness (µp = 10−4cm2V −1s−1). The Nv value of 1.25× 1021cm−3 will
. . . . . . . . . . . . . . . . . . . .
be used in all following calculations.

Figure 3.8: Mott-Gurney convergence – The theoretical current-voltage characteristic
calculated from the 1-D DD-SCLC model of the trap-free (ft = 0) hole-
only device obeys the power law J ∝ V m, demonstrating the transition
from the diﬀusion aﬀected regime of the ﬁtting exponent m = 1.52 to
the Mott-Gurney regime of the ﬁtting exponent m = 2 with increasing
applied voltage.
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

Figure 3.9: Schematic illustration of a three-dimensional hole-only device of tip-plane
geometry. The numerical calculation is performed using a cylindrical co-
ordinate system with the cylindrical symmetry, reducing the calculation
domain to the two-dimensional system as shown. Hole transport is then
. . . . . . . . . . . . . . . . . . .
described by a 2-D DD-SCLC model.

Figure 3.10: Schematic illustration of the spatial discretization in a cylindrical coor-
dinate system. With cylindrical symmetry, the two-dimensional device is
divided into the sub-area of size ∆ρ × ∆z with mesh points located at
each corner. Each sub-area carries the single values of dielectric constant,
hole mobility, and trap parameters. While it is homogeneous within sub-
area, each sub-area can have diﬀerent parameter values, producing the
inhomogeneity in the system.
. . . . . . . . . . . . . . . . . . . . . . . .

48

49

51

53

Figure 3.11: Schematic illustration of Gaussian Surface of (a) axis cell for the centered

mesh at z-axis (i = 0), and (b) regular cell for the centered mesh elsewhere 53

Figure 3.12: BiCGSTAB VS Inversion – The theoretical current-voltage characteris-
tic of the trap-free (ft = 0) hole-only device solved iteratively by using
BiCGSTAB module is compared with that solved directly by inversion of
a tridiagonal matrix. In both calculations, we use a 1-D DD-SCLC model
with L = 80nm, Nv = 1.25 × 1021cm−3,r = 3, µp = 10−4cm2V −1s−1.
.

Figure 3.13: The comparison of the theoretical current-voltage characteristics of the
trap-free (ft = 0) hole-only device calculated from 1-D and 2-D DD-SCLC
models.
In both models, the transport layer (orange) is homogeneous
with the constant parameters: r = 3, µp = 10−4cm2V −1s−1, and Nv =
1.25 × 1021cm−3.
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

64

65

xii

Figure 3.14: Finite-size eﬀect–Convergence of the simulated hole current with respect
to device dimensional ratio of radius to thickness, when the applied volt-
age Va = 0.1, 1, 10V . By using a 2-D DD-SCLC model, the trap-free
(ft = 0) hole-only device is modeled in tip-plane geometry (Rtip ∼ 1nm),
assuming homogeneous transport layer (r = 3, Nv = 1.25 × 1021cm−3)
with isotropic hole mobility (µp,ρ = µp,z = 10−4cm2V −1s−1). Using the
dimensional ratio R/L = 3.75 provides a well-converged solutions while
the ﬁnite-size eﬀect is negligible. All the following calculations of tip-plane
devices with an isotropic mobility will be performed with this dimensional
ratio. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

Figure 3.15: Finite-size eﬀect–Convergence of the simulated hole current with respect
to device dimensional ratio of radius to thickness, when the applied voltage
Va = 0.1, 1, 10V . By using a 2-D DD-SCLC model, the trap-free (ft = 0)
holse-only device is modeled in tip-plane geometry (Rtip ∼ 1nm), assum-
ing homogeneous transport layer (r = 3, Nv = 1.25 × 1021cm−3) with
anisotropic hole mobility (µp,ρ/µp,z = 25, µp,z = 10−4cm2V −1s−1). Us-
ing the dimensional ratio R/L = 6.25 provides a well-converged solutions
while the ﬁnite-size eﬀect is negligible. All the following calculations of
tip-plane devices with anisotropic mobility will be performed with this
dimensional ratio. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

Figure 3.16: Current as a function of vacuum layer thickness simulated using the given
device structure, compared to the current calculated from a 2-D DD-
SCLC model with insulating boundary conditions described by eq.(3.59),
when Va = 0.5, 1, 10V . At high vacuum layer thickness, the consistency
of the numerical current between the two device geometries is observed
with approximately a maximum of 1% diﬀerence. For the transport layer
(orange), we use L = 80nm, r = 3, µp,ρ = µp,z = 10−4cm2V −1s−1, an5d
Nv = 1.25 × 1021cm−3. . . . . . . . . . . . . . . . . . . . . . . . . . . . .

Figure 3.17: Illustration of a three-dimensional hole-only device of tip-plane geometry.
The numerical calculation is performed using a three-dimensional coordi-
nate system, allowing full customization of trap distribution and nanoscale
morphology in the transport layer.
. . . . . . . . . . . . . . . . . . . . .

Figure 3.18: Illustration of the spatial discretization in a three-dimensional Cartesian
coordinate system; (a) the device is divided into cube cells (h×h×h) with
mesh points located at each corner; (b) each mesh point is surrounded by
six nearest neighbors, and is adjacent to eight cells; (c) each cell containss
constant values of dielectric constant, hole mobility, and trap parameters.
The inhomogeneity is generated by varying the values of these parameters
from one cell to another. . . . . . . . . . . . . . . . . . . . . . . . . . . .

68

69

71

73

76

xiii

Figure 3.19: Gaussian surface SG as a cube of size h × h × h surrounding mesh point
(i, j, k) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

76

Figure 3.20: The comparison of the theoretical current-voltage characteristics of trap-
free (ft = 0) hole-only devices calculated from 1-D and 3-D DD-SCLC
models.
In both models, the transport layer (orange) is homogeneous
with the constant parameters: r = 3, µp = 10−4cm2V −1s−1, and Nv =
1.25 × 1021cm−3. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

Figure 3.21: The schematic representation of tip contacts in a 3-D DD-SCLC model,
constructed on the 2nm×2nm square grid as (a) an approximately circular
geometry of 6nm radius and (b) a square geometry of 10nm side-length. In
a 2-D DD-SCLC model, the tip contact is assigned as the circular area of
6nm radius (red dashed circle). Due to the diﬀerence in the discretization,
the tip models in the 3-D model are slightly smaller that in the 2-D model,
such that A3−D(a)/A2−D = 0.85 and A3−D(b)/A2−D = 0.88, respectively.

Figure 3.22: The comparison of the theoretical current-voltage characteristics of the
trap-free (ft = 0) hole-only devices calculated from 2-D and 3-D DD-
SCLC models, in which the tip electrode is approximated by a circular
contact of 6nm radius, as illustrated in ﬁgure 3.21(a). In both models,
the transport layer (orange) is homogeneous with the constant parameters:
r = 3, µp = 10−4cm2V −1s−1, and Nv = 1.25 × 1021cm−3. . . . . . . . .

Figure 3.23: The comparison of the theoretical current-voltage characteristics of the
trap-free (ft = 0) hole-only devices calculated from 2-D and 3-D DD-
SCLC models, in which the tip electrode is approximated by a circular con-
tact of 6nm radius and a square contact of 10nm side length, respectively.
(See ﬁgure 3.21(b)) In both models, the transport layer (orange) is homo-
geneous with the constant parameters: r = 3, µp = 10−4cm2V −1s−1,
and Nv = 1.25 × 1021cm−3.
. . . . . . . . . . . . . . . . . . . . . . . . .

Figure 4.1: Based on the simpliﬁed SCLC model of a trap-free hole-only device in the
planar geometry, the analytical proﬁles of electric potential (solid black
line) and hole carrier density (dashed red line) are demonstrated using
. . . . . . . . . . . . . . . . . . . . . . . . . . .
non-dimensional scales.

85

86

87

88

93

xiv

Figure 4.2: Plot of the analytical ratio of diﬀusion current density to drift current
density of a planar ﬂow, trap-free, hole-only device, Jp,di/Jp,dr, versus
the dimensionless coordinate, x/L. While the drift current density is the
immediate analytic solution of the simpliﬁed SCLC model as known as
Mott-Gurney equation (2.13), the diﬀusion current density is obtained
from the further derivation of the analytic solutions of the model. The
extremely large value of diﬀusion current near the injecting electrode (x =
0) suggests the diﬀusion assists hole injection at the injecting contact.
.

Figure 4.3: The theoretical proﬁle of electric potential and hole carrier density as a
function of the distance from the injecting electrode for planar geometry,
generated using the 1-D DD-SCLC model described in section 3.3, using
parameters given in table 4.1. The proﬁles exhibit three distinct transport
regimes, plotted for three diﬀerent applied voltages Va = 0.1, 1, 10V , from
top to bottom.
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

95

98

Figure 4.4: The theoretical proﬁle of electric potential and hole carrier density at
the center of a circular tip (ρ = 0) as a function of the distance from
the injecting electrode for tip-plane geometry, generated by the 2-D DD-
SCLC model described in section 3.4, using parameters given in table 4.1.
The proﬁles exhibits three distinct transport regimes, plotted for three
diﬀerent applied voltages Va = 0.1, 1, 10V .

. . . . . . . . . . . . . . . . . 100

Figure 4.5: Position of virtual anode (dashed line) and virtual cathode (solid line)
as a function of applied voltage for planar (black) and tip-plane (red)
geometries, as seen in ﬁgure 4.3 and ﬁgure 4.4 respectively. The inset
shows the convergence of the virtual tip position.

. . . . . . . . . . . . . 102

Figure 4.6: The trap inﬂuenced current-voltage characteristic of a hole-only device in
planar geometry obeys the power law J ∝ V m, exhibiting three regimes:
Ohm’s law m = 1 (red), trap-limited SCLC m > 2 (green), and trap-ﬁlled
SCLC m = 2 (blue).

. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 106

Figure 4.7: The theoretical current-voltage characteristics of planar hole-only devices
in the SCLC regime obey the power law (I ∝ V m) when (a) ft = 0.005
with varying rt yielding the exponent in the range m = 2.23 − 4.71,
and (b) rt = 0.30 with varying ft yielding the exponent in the range
m = 1.89 − 3.72

. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 108

Figure 4.8: Dependence of the exponent m on the two trap parameters; trap den-
sity characterized by parameter ft and trap temperature characterized
by parameter rt for hole-only devices in planar geometry based on (a)
the simpliﬁed SCLC model as described by eq.(2.27), (b) 1-D DD-SCLC
model as seen in ﬁgure 4.7 and described by eq.(4.3).

. . . . . . . . . . 109

xv

Figure 4.9: The trap inﬂuenced current-voltage characteristic of a hole-only device in
the tip-plane geometry obeys the power law J ∝ V m, exhibiting three
regimes: Ohm’s law m = 1 (red), trap limited SCLC m > 2 (green), and
trap ﬁlled SCLC m = 2 (blue). Hole mobility is assumed to be isotropic.

111

Figure 4.10: The theoretical current-voltage characteristics of tip-plane devices in the
SCLC regime obey the power law (I ∝ V m) when (a) ft = 0.005 with
varying rt yielding the exponent in the range m = 2.89−5.52, and (b) rt =
0.30 with varying ft yielding the exponent in the range m = 2.00 − 4.19.
Hole mobility is assumed to be isotropic.

. . . . . . . . . . . . . . . . . 113

Figure 4.11: The theoretical current-voltage characteristics of tip-plane devices in the
SCLC regime for a range of mobility anisotropic ratios, µρ/µz, with ﬁxed
trap parameters ft = 0.005, rt = 0.30.

. . . . . . . . . . . . . . . . . . . 114

Figure 4.12: Dependence of the exponent m on the exponentially distributed trap DOSs
characterized by parameters ft and rt for hole-only devices in tip-plane
geometry in the SCLC regime, resulting from extensive calculations for
theoretical current-voltage characteristics as demonstrated in ﬁgure 4.10
and described by the semi-empirical expression (4.4).

. . . . . . . . . . 115

Figure 4.13: Schematic illustration of crystallized P3HT nanoﬁbers in three possible
orientation schemes, including (a) ‘face-on’, (b) ‘edge-on’ and (c) ‘end-on’
orientations, corresponding to the parallel and perpendicular alignments
of self-organized P3HT lamellae to the substrate, respectively.

. . . . . . 117

Figure 4.14: Schematic illustration of device structures with the tip-plane electrode
conﬁguration to study the morphological eﬀect of an isolated ﬁber with
(a) ‘face-on’ and (b) ‘edge-on’ orientations. . . . . . . . . . . . . . . . . . 120

Figure 4.15: Based on the hole-only devices with an isolated P3HT ﬁber shown in ﬁg-
ure 4.14, the theoretical current-voltage characteristics when the tip is in
contact with the ﬁber are plotted with varying the ﬁber mobility along the
polythiophene backbone direction (µB). The theoretical current-voltage
characteristic of an amorphous device is plotted as the reference to demon-
strate the current enhancement in the presence of a single nanoﬁber. The
inset shows the same data plotted in the double logarithmic scale.

. . . 122

xvi

Figure 4.16: The vertical (a-e) and lateral (f-j) current of the semiconducting P3HT
devices when the biasing voltage Va = 1V is applied to the conduct-
ing tip (red dot) and the bottom electrode is grounded. Three diﬀerent
types of morphologies have been used in the device simulation: (a)&(f)
an amorphous (no-ﬁber) domain; (b-c)&(g-h) an isolated ‘face-on’ ﬁber
embedded in an amorphous domain; and (d-e)&(i-j) an isolated ‘edge-on’
ﬁber embedded in an amorphous domain. In the last two morphologies,
the comparison of current ﬂows when tip is and is not in contact with the
ﬁber is demonstrated.

. . . . . . . . . . . . . . . . . . . . . . . . . . . . 127

Figure 4.17: Simulating the c-AFM measurements, the dependence of the vertical cur-
rent as a function of tip position on the proﬁle trace crossing the ‘face-on’
ﬁber along; (a) the polythiophene backbone direction and (b) the alkyl
side chain direction, at the applied voltages Va = 0.1, 1, 10V . The proﬁle
trace is indicated by the dashed line shown in the inset. . . . . . . . . . . 129

Figure 4.18: Simulating c-AFM measurements, (a) the dependence of the vertical cur-
rent as a function of tip position on the proﬁle trace crossing the ‘edge-on’
ﬁber along the polythiophene backbone direction and (b) the dependence
of the vertical current as a function of tip-electrode separation, at the
applied voltage Va = 0.1, 1, 10V . . . . . . . . . . . . . . . . . . . . . . . . 131

Figure 4.19: The 600nm×600nm cross-sections of three model morphologies consisting
of ‘edge-on’ ﬁbers aligned in-plane with (a) random spacing (ﬁber fraction
= 0.40), (b) random spacing and random orientation (ﬁber fraction =
0.40) and (c) random spacing and random length (ﬁber fraction = 0.32).

133

Figure 4.20: Theoretical current-voltage characteristics of hole-only devices in the tip-
plane geometry with three model morphologies (see ﬁgure 4.19) simulated
when the conducting tip is in contact with the ﬁber and when it is not,
as labeled by ‘ON’ and ‘OFF’, respectively. The inset shows the same
data plotted in the double logarithmic scale. Despite the diﬀerences in
model morphologies, the total currents simulated from the devices with
an ‘edge-on’ ﬁber network increase by approximately 50% compared to
those simulated from the devices with a single With the presence of an
‘edge-on’ ﬁber network, the total currents, despite the diﬀerences in model
morphologies and tip locations, increases by approximately 50% compared
to those simulated from the devices with a single ‘edge-on’ ﬁber.

. . . . 134

xvii

Figure 5.1: Reported in ref. [1], (a) an image of the topography of the annealed sample
for ﬁxed-spot current-voltage (I-V ) measurement; (b) I-V spectra with
the tip located on/oﬀ the ﬁbers at locations labeled in (a). The inset of
(b) shows log-log plots of the I-V spectra (Log|I| vs. Log|V |) from the
negative bias regime. (c) Typical I-V curves in the negative sample bias
for annealed samples (blue) and non-annealed samples (green), and the
ﬁtting to power laws IαV m (red). For annealed samples, the I-V spectra
have three distinct segments: Ohms law ( m = 1), space-charge limited
transport with shallow traps (m = 2 − 2.25), and a higher slope regime
(m > 5) with increase of the applied bias. In non-annealed samples, trap
eﬀects are pronounced in the SCLC regime as the exponent is signiﬁcantly
higher (m = 4 − 5).

. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 142

Figure 5.2: The schematic representation of the device structure (a) and the corre-
sponding energy alignment (b) of the hole-only device. This device be-
haves similarly to a semiconductor diode. When the device is positive bi-
ased, no charge carrier transport presents. Conversely, under the negative
bias, holes are injected from Pt probe into the system, traveling through
the layer of P3HT and PEDOT:PSS, and collected at the ITO contact;
while electrons injected from ITO cannot transport into the system due to
the layer of PEDOT:PSS. Therefore, this is the hole only injection system. 145

Figure 5.3: Based on the setup of c-AFM measurements [1], the theoretical current-
voltage characteristics in the SCLC regime are well-described by the semi-
empirical expression (4.4) when (a) trap parameters is ﬁxed (ft = 0.001
rt = 0.30) with varying mobility anisotropic ratios µρ/µz; (b) mobility
is isotropic and ft = 0.005 with the varying rt yielding the exponent in
range m = 2.31 − 5.19; and (c) mobility is isotropic and rt = 0.35 with
varying ft yielding the exponent in range m = 1.96 − 3.61.

. . . . . . . 149

Figure 5.4: Dependence of the exponent m on the parameters ft and rt for hole-
only devices in the SCLC regime, resulting from extensive calculations of
theoretical current-voltage characteristics demonstrated in ﬁgure 5.3 [1].

150

Figure 5.5: The representative current-voltage (I-V ) response from the annealed sam-
ple (solid circle) is ﬁtted to the theoretical I-V relation (solid red line) in
the low trap regime with two pairs of trap parameters that result in the
same slope: (a) ft = 0.0005 and rt = 0.35, (b) ft = 0.001 and rt = 0.40.
Shown in the inset are the same plot in linear scale. (c) A histogram of
the normalized vertical hole mobility extracted from 48 I-V datasets on
the annealed P3HT ﬁlm, with the comparison of an isotropic model and
an anisotropic model with µr/µz = 25. The mobility is normalized to the
maximum value of each case.

. . . . . . . . . . . . . . . . . . . . . . . . 152

xviii

Figure 5.6: The representative current-voltage response from the non-annealed sample
(solid circle) is ﬁtted to the theoretical current-voltage relation (solid red
line) in the low trap regime with two pairs of trap parameters that result
in the same slope: (a) ft = 0.0075 and rt = 0.25, (b) ft = 0.01 and
rt = 0.25. Shown in the inset are the same plot in linear scale. . . . . . . 154

Figure 5.7: The 600nm × 600nm cross-sections of the representative morphology at
the vertical position z in the device thickness direction demonstrating the
crystallized P3HT ﬁbers of 20nm width (green wires) aligned with random
spacing in the amorphous region (brown area). The ﬁber fraction of cross-
section surface area is 0.33. The bottom electrode is placed at z = 0nm,
and the locations of tip contact are marked on the z = 80nm cross-section
for the on-ﬁber and oﬀ-ﬁber c-AFM simulations.

. . . . . . . . . . . . . 158

Figure 5.8: Adapted from [1], the on-ﬁber and oﬀ-ﬁber current-voltage responses of
the annealed P3HT sample are well described by the theoretical current-
voltage characteristics simulated from the representative 3-D morphology
of an ‘edge-on’ ﬁber network shown in the inset on the bottom-right corner.
(see ﬁgure 5.7 for cross-section morphology) To described the anisotropic
transport in ﬁbers (green), the in-plane mobility is 100 times larger than
the vertical mobility and the background mobility of amorphous region
(brown). The background mobility is ﬁtted to be 7.78× 10−5cm2V −1s−1
with trap parameters ft = 1.00 × 10−3 and rt = 0.30. In the inset on the
top-left corner, the vertical current ﬂow when the c-AFM tip is located
on the ﬁber is visualized.

. . . . . . . . . . . . . . . . . . . . . . . . . . 159

xix

KEY TO SYMBOLS

Boltzmann’s constant

The fundamental charge

Temperature
Thermal potential, deﬁned as Vt = kbT /q ≈ 0.026V at room temperature.
Externally applied voltage

Electric potential

Electric ﬁeld

Dielectric constant

The vacuum permittivity

Total intrinsic hole carrier density, consisting of free and trap hole densities, i.e.,
po = po,f + po,t.

Total injected hole carrier density, consisting of free and trap hole densities, i.e.,
p = pf + pt

Total holes

Hole mobility

The eﬀective hole mobility in the vertical direction

Hole diﬀusion coeﬃcient

Total hole current, consisting of drift and diﬀusion currents, i.e., Ip = Ip,dr +Ip,di

Total hole current density, consisting of drift and diﬀusion current densities, i.e.,
Jp = Jp,dr + Jp,di
The Bernoulli function, deﬁned as B(x) =

x

ex − 1

.

The eﬀective density of states in the HOMO

Total trap density

kb

q

T

Vt

V, Va

ψ
E

r

o

po

p

Qp

µp

µo

Dp

I, Ip

J, Jp

B

Nv

Nt

Tt, Et

Trap temperature and trap energy (Et = kbTt)

ft, rt
S

Normalized trap parameters, deﬁned as ft = Nt/Nv and rt = T /Tt, respectively.

The spatial distribution of traps

xx

SG

L, Lz

Lx, Ly

Ltip

R

Gaussian surface

Device thickness
Device side lengths in the x− and y− directions in Guassian coordinate systems

Square tip side length

Device radius in the cylindrical coordinate systems

Circular tip radius

Generalized tip domain

Rtip
Dtip
xc
kJ , kI , m The ﬁtting parameters of the power relation J = µokJ V m and I = µokI V m.

Characteristic coordinate

xxi

Chapter 1

Introduction

1.1 Motivation

With the drastic growth of global energy consumption as well as worldwide environmental

concern, renewable energy has drawn an extensive interest in the past decades. Solar energy

is widely recognized as one of the most essential energy sources for the future sustainable

economy. Photovoltaic energy conversion has been developing to harvest solar energy and to

eﬃciently generate electricity. As well as well-known technologies, such as the conventional

inorganic and the emerging perovskite photovoltaic devices, organic photovoltaic (OPV)

devices are promising, owing to several attractive features [2, 3, 4, 5, 6, 7, 8]. The fabrication

processes of OPV devices involve low-temperature and low-cost deposition techniques that

easily translate to large-scale production. Since their manufacturing speed is very high and

their thermal budget is very low, the OPV devices possess advantages over their inorganic

counterparts in the aspect of ecological footprint and energy payback time, meaning that

organic materials are more environment-friendly and they need shorter operation time to

recover the energy spent in their fabrications, respectively. Furthermore, OPV devices are

light weight, ﬂexible and semi-transparent, oﬀering versatile applications. Over the course

of time, the power conversion eﬃciency of OPVs has vastly improved.

Since the basic structure of the OPVs was established as early as the mid-1950s consisted

1

of a single layer of organic semiconductor with electrodes attached to its top and bottom

surfaces [9], we observed the improvement of the device performance with respect to the

large variation of the organic materials and the electrodes [10, 11, 12, 13, 14, 15, 16]. After

that, double layered OPV devices were introduced in 1986 [17], and the bulk heterojunction

(BHJ) design came about in 1995 [18]. It was found that the new structures can boost the

performance of OPV devices greatly. Many other works have been done in this direction

to achieve ever better optimization in OPV structures and morphologies, yielding over 10%

power conversion eﬃciency recently [19, 20]. However, the underlying mechanisms regarding

how the device operates is not yet well understood due to the large variation in the device

materials and complex interactions among them. A clear understanding of this matter

will not only push forward the study of OPV materials, but will also shed light on the

investigation of other related devices, such as organic light emitting diodes (OLEDs) [21, 22,

23] and organic ﬁeld eﬀect transistors (OFETs) [24, 25].

Charge carrier transport is crucial in organic semiconductors, and has been investigated

broadly via various methods such as time-of-ﬂight (TOF) methods, ﬁeld-eﬀect transistor

(FET) conﬁgurations, pulse-radiolysis time-resolved microwave conductivity (PR-TRMC),

as well as the measurement of current-voltage characteristics in the space-charge limited

current (SCLC) regime [26]. Among them, the method involving the SCLC is particularly

interesting because of its simplicity and ﬂexibility. The experiment is based on the measure-

ment of SCLC in a diode conﬁguration, the so called planar (plane-parallel) geometry, in

which an organic layer is fabricated on the substrate in contact with top and bottom elec-

trodes [26, 27, 28]. With the appropriate analytical and/or numerical model, these current-

voltage data in the SCLC regime can be analyzed to determine the eﬀective carrier mobility,

as well as some other physical parameters, such as energy barriers, trap density, and trap

2

energy states.

Theoretically, the drift-diﬀusion (DD) model describing bulk transport of charge carriers

in the SCLC regime is governed by three basic equations, including Poisson’s equation, the

continuity equation at steady state, and the drift-diﬀusion equations of the charge carriers.

The early theoretical work is vastly simpliﬁed by ignoring the contribution of charge diﬀusion

and only considering the drift driven transport, such that the theoretical model can be solved

analytically. The analytical expressions describing the dependences of trap-free SCLC and

trap-limited SCLC on applied voltage in a planar geometry are rather simple, and have been

widely used to determine charge carrier mobilities in experiments [26, 27, 28]. The simpliﬁed

SCLC theory for the planar-ﬂow current will be reviewed in section 2.3.

Conductive atomic force microscope (c-AFM) has been used extensively to study organic

semiconductor devices in the past decade [29, 30, 31, 32, 33, 34]. One of the eﬀorts is to use

the local current-voltage data measured by c-AFM to determine the local charge carrier

mobility at speciﬁc locations of the ﬁlm. Many studies indicated the importance of c-AFM

tip geometry, and the existing analytical solutions for the conventional planar geometry, that

is the Mott-Gurney equation, can no longer be applied [29, 35]. Reid et al. [32] introduced

a semi-empirical formula that made a correction due to the tip geometry, which will be

reviewed in section 2.4.

However, the topic of SCLC in an organic semiconductor is more complicated than that

described by the semi-empirical formula of [32] and its predecessor, the Mott-Gurney equa-

tion. This is due to the fact that in the results of the simpliﬁed SCLC model, the eﬀect

of charge diﬀusion is excluded from the theory. In many situations ignoring the diﬀusion

current can lead to the misinterpretation of the experiment results. With this consideration

in mind, the full description of the DD-SCLC model is not only crucial to characterize SCLC

3

in experiments, but also to understand the eﬀect of charge diﬀusion on the simpliﬁed SCLC

model.

1.2 Overview

The central focus of this dissertation is to develop the computational methods of DD-SCLC

model to describe charge carrier transport in hole injection systems. We will extend the

simpliﬁed model to consider the mechanism of charge diﬀusion and the inhomogeneity of

devices including the c-AFM tip geometry, nanoscale morphologies, and spatially-varying

trap distributions. This can be seen as a building block to gain more insight into the

SCLC transport and to better characterize the experiments. To do so, this dissertation is

constructed as follows.

Chapter 2 explains the theory of charge carrier injection in an organic semiconductor,

emphasizing the SCLC model for hole-only devices. To gain insight into the current stage of

the SCLC model, this chapter will give an overview of the theoretical works on the simpliﬁed

SCLC model that have been done in the past, including the simpliﬁed model of trap-free

SCLC and trap-limited SCLC in a planar geometry, as well as the semi-empirical model of

trap-free SCLC in a tip-plane geometry.

In chapter 3, we present our development in the DD-SCLC model for hole-only devices.

The model will be described in detail, along with the numerical methods used to achieve

numerical solutions of the models in one-, two- and three- dimensional systems. All one-, two-

and three- dimensional models are veriﬁed to be consistent with each other, and validated

against the simpliﬁed SCLC theory in the drift-dominated SCLC regime.

In chapter 4, we proceed to the simulation of hole-only devices in the SCLC regime

4

using the numerical tools we successfully developed. Theoretical study on various topics of

SCLC will be presented, including the diﬀusion eﬀect in planar and tip-plane geometries, the

trap eﬀect in planar and tip-plane geometries, as well as morphological eﬀects, particularly in

the tip-plane geometry. The ﬁbrous morphology of semi-crystaline poly(3-hexylthiophene)

(P3HT) thin ﬁlms will be used as a model system to demonstrate the mechanism enabling

c-AFM to visualize the surface morphology in experiments [1, 31, 36, 37, 38]. The model

morphologies used in the studies are generated by Dr. Daniel Olds.

In chapter 5, two- and three- dimensional DD-SCLC models will be used to elucidate the

c-AFM experiments carried out by Professor Zhang’s group at Michigan State University [1].

The experiment aims to investigate the eﬀect of thermal annealing on the morphology of

P3HT thin ﬁlms. The statistical variation in local electrical conductivity of P3HT ﬁlms was

characterized by measuring the c-AFM current-voltage responses across the ﬁlm surface.

These current-voltage data are analyzed using the two- and three- dimensional DD-SCLC

models, in which the eﬀect of current spreading beneath the c-AFM tip, the eﬀect of traps,

the mobility anisotropy and the nanoscale morphology are incorporated.

Chapter 6 presents the conclusions of this dissertation and an outlook on future research

directions.

5

Chapter 2

Theory of Charge Carrier Injection -

Background

2.1 Charge Injection in Organic Semiconductors

Organic semiconductors are typically wide band gap semiconductors, where the charge carrier

concentration in thermal equilibrium at room temperature is very low, and so is its electrical

conductivity. The shortage of charge carriers can be overcome by incorporating an extrinsic

source of free charges, including doping, dissociation of photogenerated electron-hole pairs

(excitons), and charge injection at the electrodes. The topic of charge injection has drawn

a great deal of attention, due to its ability to probe the intrinsic properties of the material,

such as charge carrier mobility, defects, and traps.

Charge injection experiments are often done in the dark environment to eliminate the

exciton eﬀect. Dark charge injection is widely used to investigate electronic properties of

organic semiconductors [27, 29, 31, 35, 39, 40, 41, 42]. The experiments involve the measure-

ment of current-voltage response in a diode device conﬁguration (planar or plane-parallel

geometry), in which an organic layer is fabricated on the substrate in contact with top and

bottom electrodes. When voltage is applied, the amount of electric current ﬂowing through

the organic semiconductor device primarily depends on two basic processes: (i) the charge

6

injection at the electrode/organic-semiconductor interfaces, and (ii) the charge transport in

the bulk organic semiconductor. Two mechanisms thus control the electric current in the two

distinct regimes: injection-limited current (ILC), and space-charge-limited current (SCLC),

respectively.

In the ILC regime, the dark injection current is controlled by the eﬃciency of charge

injection at the electrodes. An injection barrier often occurs when a metal electrode comes

into contact with an organic semiconductor, preventing charges from being injected into the

organic material, hence suppressing the current. In order to resolve this limitation, one can

select proper electrodes so that the injection barrier at the electrode/organic-semiconductor

contact is negligible. The ideal injecting contact is essentially an unlimited charge reservoir.

This contact, known as ohmic contact, can eliminate the eﬀect of ILC.

On the other hand, in the SCLC regime, the dark injection current is governed by charge

carrier transport.

In general, organic semiconductors are highly disordered, especially in

the amorphous state, and hence the charge carrier mobility is much smaller than that of

conventional inorganic semiconductors. When charges are injected into the material, they

do not spread out eﬃciently. Instead, they form a cloud of space charge around the contact.

The presence of space charge distorts the local electric ﬁeld strength, and particularly reduces

the electric ﬁeld at the injecting electrode. When the electric ﬁeld strength generated by

the induced space charges is in equilibrium with the electric ﬁeld strength from the applied

voltage, the device reaches steady state with constant electric current, i.e., the so-called

SCLC.

In addition to space charge, the current in the SCLC regime can also be aﬀected by other

bulk properties of the organic semiconductor, including spatially-varying charge mobility,

nanoscale morphologies, as well as trap states. A complete model including all the relevant

7

bulk properties can provide more insights into the underlying charge transport mechanism in

organic semiconductors. In the rest of this chapter, we will review the simpliﬁed theoretical

model considering only the space charge eﬀect. An extended theoretical model will be the

focus of the later chapter.

2.2 Fundamental Equations to Describe SCLC in Hole

Injection System

To eliminate the interaction of electrons and holes, charge injection experiments are mostly

unipolar, where one carrier type is dominant. Whether it is the injection of electrons or holes

depends on the energy level alignment of the devices, the energy level mismatch between the

device material and the metal electrode in particular. The ground-state electronic structure

of organic semiconductors is described by the sequence of molecular orbitals ﬁlled by elec-

trons up to the highest occupied molecular orbital (HOMO), and from the lowest unoccupied

molecular orbital (LUMO) onward are empty. In analogy to the traditional band structure

of inorganic semiconductors, the HOMO and LUMO energies take the role of valance and

conduction band edges, respectively. When an electron is injected to the system, it will

occupy that empty orbital at the lowest energy level, i.e., the so called LUMO. Similarly,

the injection of a hole is equivalent to an electron being removed from the HOMO and ex-

tracted at the electrode. Since the topic of hole injection system is the central focus of this

dissertation, the following derivations and notations will be emphasized on hole transport,

and note that the reverse will be equally true for electron injection devices. A thorough

discussion of an electron injection device has been provided by Lampert and Mark [43].

The behavior of hole SCLC at steady state is governed by Poisson’s equation, the conti-

8

nuity equation in steady state, and the drift-diﬀusion equation:

(/q)((cid:126)∇ · (cid:126)E) = p − po,

(cid:126)∇ · (cid:126)Jp = 0,

(cid:126)Jp = (cid:126)Jp,dr + (cid:126)Jp,di = qµppf

(cid:126)E − qDp (cid:126)∇pf ,

(2.1)

(2.2)

(2.3)

where  is the dielectric constant, q is the fundamental charge, µp is the eﬀective hole mobility,

and Dp is the hole diﬀusion coeﬃcient.

Poisson’s equation (2.1) describes the instantaneous electric ﬁeld produced by space

charges in the system as shown on the right hand side of the equation. po = pf o + pto

is the total hole density in thermal equilibrium, consisting of the free hole density pf o and

trapped hole density pto. After voltage is applied, a number of holes are injected into the

system, thus increasing the total hole density p. Similarly, p((cid:126)r) = pf ((cid:126)r) + pt((cid:126)r), where pf

and pt are the actual, spatially varying free and trapped hole density in the steady-state

injection. As discussed in section 2.1, the pristine (undoped) organic semiconductor has a

relatively low free carrier concentration at thermal equilibrium, and most of the charges that

carry the current must be injected. Therefore, we can safely omit the thermal equilibrium
hole density p0 from Poisson’s equation, such that p − po (cid:39) p.

In addition, the transport of charges can be described by the hole current density Jp.

There are two types of hole current, drift current Jp,dr and diﬀusion current Jp,di, driven by
electric ﬁeld (cid:126)E and the gradient of free hole density, respectively. Both of them are considered

in the drift-diﬀusion equation (2.3), while the total current density is conserved according

to the continuity equation (2.2).

9

By simultaneously solving the set of three governing equations (2.1)-(2.3) together with

appropriate boundary conditions, the theoretical current-voltage characteristics of hole-only

devices can be determined, as well as the theoretical proﬁles of the electric ﬁeld and hole

carrier density which can not be measured directly from experiments. We will ﬁrst intro-

duce the closed form solution to the simpliﬁed equations, and dedicate the next chapter to

extensive computational studies of the problem in its complete unsimpliﬁed form. As we

shall see, the comparison between the two methods demonstrates the importance of the full

computation.

2.3 The Simpliﬁed Theory of SCLC: Planar (Plane-

Parallel) Geometry

The simpliﬁed SCLC theory of unipolar injection is based on two fundamental assumptions

[43]: (i) the injecting contact is ohmic, and (ii) the diﬀusion current is negligible.

In the ﬁrst assumption, the ohmic contact is deﬁned as a low-resistance contact, behaving

as an inﬁnite source of charges available for injection. As we discussed in section 2.1, the

electric current of the dark charge injection is controlled by either the injection barrier at the

electrode, or space charges in bulk organic semiconductors. By applying the approximation

(i), the simpliﬁed theory is independent of the contact properties at the electrode/organic-

semiconductor interface, and only the bulk transport properties are considered.

Next, the assumption (ii) is directly related to the charge transport. When the electric

ﬁeld is suﬃciently strong, the drift current dominates, and we can safely omit the contribu-

tion of charge diﬀusion. It is also important to note that the two assumptions do not hold

near contacts where the diﬀusion current is signiﬁcant. Therefore, while the theory provides

10

a good approximation to the charge transport in bulk, namely the current-voltage charac-

teristics of organic semiconductor devices, one shall bear in mind that the result contains

approximations that are invalid near the electrodes.

In this section, the current-voltage characteristics of the hole-only devices in planar ge-

ometries are investigated, where the fully three-dimensional partial diﬀerential equations

(PDEs) can be reduced to one-dimensional ordinary diﬀerential equations (ODEs). We will

consider an organic semiconductor thin ﬁlm of thickness L, in which holes are injected into

the device at the anode x = 0, and extracted from the device at the cathode x = L. In order

to do so, the device is grounded at the anode ψ(0) = 0, and negative bias is applied at the
cathode ψ(L) = −V . Based on the assumptions above, the set of governing equations can

be written as

(/q)(dE/dx) = pf + pt,

dJp/dx = 0 or Jp = constant,

Jp = qµppfE,

(2.4)

(2.5)

(2.6)

where (cid:126)E = E ˆx and (cid:126)Jp = Jp ˆx. Note an additional boundary condition arises from assumption

(i). In order to have ﬁnite current yet inﬁnite charge density at the injecting electrode, the

electric ﬁeld must vanish there, so that

E = 0 at the injecting electrode.

(2.7)

This is the mathematical representation of an ohmic contact for the simpliﬁed SCLC theory.

Next, the problem of SCLC in the planar device geometry will be solved using the simplest

11

model where the device is assumed to be free of traps. More realistic solution will be given

later which takes into account device traps that assume exponentially distributed density of

trap states.

2.3.1 Trap-Free SCLC

The simplest model of unipolar SCLC does not consider traps, such that pt = 0. Substituting

the free hole density pf from Poisson’s equation (2.4) into equation (2.6), we can derive the

drift current as:

Jp = µpE(dE/dx) = (µp/2)(dE 2/dx) ⇒ dE 2/dx = 2Jp/µp.

(2.8)

The general solution has the form

E 2(x) = ax + b.

(2.9)

Applying the ohmic boundary condition at the injecting electrode (in eq.(2.7)), we get

√

E(x) =

ax1/2.

(2.10)

To calculate the electric potential across the device, this equation may be integrated as

√
E(x)dx = −(2/3)

ax3/2 + c

(2.11)

follows.

ψ(x) = −

(cid:90)

It is important to note that this electric potential is based on the ﬂat band condition, where

there are no energy level shifts near the interface at which metal and organic semiconductor

12

come into contact. However, the validity of this ﬂat band picture needs to be thoroughly

examined, especially when the injecting electrode is ohmic.

In analogy to an inorganic

semiconductor picture, when forming an ohmic contact at short-circuit condition (V = 0),

charge transfer from the ohmic electrode into the organic semiconductor occurs, resulting in

the formation of a space-charge layer near the contact at equilibrium. The contact potential

then inﬂuences the bending of the valance band edge, i.e., the so called HOMO band edge in

the case of organic semiconductors. The simpliﬁed theory clearly neglects the band bending

condition. Nevertheless it remains consistent with the initial assumption to exclude contact
behavior. Applying the boundary conditions ψ(x = 0) = 0 and ψ(x = L) = −V , the integral

constants can be determined as

a = (9/8)(V 2/L3) and c = 0.

(2.12)

Substituting the constant a into the solution (2.10) and the combined diﬀerential equation

(2.8), we have the hole current

Jp =

9
8

µp

V 2
L3 ,

(2.13)

which is known as Mott-Gurney equation, or the Child Law for unipolar SCLC [43]. The trap-
free SCLC can be observed as a slope 2 regime in a LogJp − LogV plot. This characteristic

is clearly observed in various unipolar injection experiments, and can be used to evaluate

the eﬀective charge carrier mobility of organic semiconductors [27, 39, 40, 41, 42].

Often the anode and cathode are made of diﬀerent materials with diﬀerent work functions,

leading to a built-in ﬁeld across the device. Consider a hole-only device with asymmetric

contacts, such that the anode is ohmic while the cathode is non-ohmic with barrier potential

equivalent to built-in potential Vbi. In the high voltage regime V > Vbi, as suggested by

13

Lampert and Mark [43], equation (2.13) shall be modiﬁed to handle this case by replacing
V with V − Vbi.

The equations for the electric ﬁeld, the electric potential, as well as the hole density can

be derived as follows,

E(x)
EΩ

=

3
2

(cid:16) x

(cid:17)1/2

L

,

ψ(x)
−V

=

(cid:16) x

(cid:17)3/2

L

(cid:16) x

(cid:17)−1/2

L

,

pf (x)

¯pf

=

1
2

,

(2.14)

where EΩ is the spatial average of the electric ﬁeld which is in the form of the ohmic ﬁeld,

EΩ = (1/L)

E(x)dx = V /L,

(2.15)

(cid:90) L

0

(cid:90) L

0

and ¯pf is the spatial average of the hole carrier density,

¯pf = (1/L)

pf (x)dx = (3/2)(V /qL2),

(2.16)

2.3.2 SCLC with The Exponential Distribution of Trap States

Trap states are commonly observed in pristine organic semiconductors with poor crystallinity,

and charge transport is often hindered by the presence of such traps. Because of the struc-

tural disorder in organic semiconductors, traps do not have a single well-deﬁned energy state.

Rather a Gaussian distribution is employed to describe a single set of traps, whereas an ex-

ponential distribution is often used for multiple sets of traps [43]. Recently, Campbell et al.

proposed an alternative trap model for electrons and holes as the tail states of the approx-

imately Gaussian distributed LUMO and HOMO density of states respectively [44]. These

tail states are pseudo-exponential: the further the tail is away from the center of the distri-

14

bution toward the energy gap, the deeper the trap it represents. An exponential distribution

is then a reasonable approximation for the distribution of traps in organic semiconductors.

The ﬁeld-independent exponentially distributed traps are characterized by the total trap

density Nt and the trap energy Et, which represents their energetic depth for thermally

activated detrapping, Et = kbTt where kb is Boltzmann’s constant. The distribution of the

trap states, i.e., the density of hole trapping states per unit energy, is deﬁned as

Gt(E) = (Nt/kbTt) e(−E/kbTt).

(2.17)

The density of trapped holes is

(cid:90) ∞

Ef (x)

pt(x) =

Gt(E)dE = Ntexp(−Ef (x)/kbTt),

(2.18)

where Ef (x) is the quasi Fermi level of holes at distance x from the anode. As long as

the temperature is suﬃciently low (T << Ef (x)/kb), all the energy states below Ef (x)

are occupied, while those above Ef (x) are empty. Assuming Ef (x) >> kbT , the free hole

density is approximately

pf (x) = Nvexp(−Ef (x)/kbT ),

and the trapped hole density can be written in terms of the free hole density as

−1/l
pt(x) = NtN
v

pf (x)1/l,

(2.19)

(2.20)

where Nv is the eﬀective DOSs in the HOMO, and l = Tt/T >> 1.

In the particular case when trapped holes are dominant (pt >> pf ), the local electric

15

ﬁeld in bulk organic semiconductors are controlled by the density and the energy distribution

of the trapping states. Then Poisson’s equation (2.4) can be rewritten as

(/q)(dE/dx) = pt(x).

(2.21)

Combining the modiﬁed Poisson’s equation (2.21) with the drift current model in equation

(2.6), we have

(/q)(dE/dx) = NtN

−1/l
v

(Jp/µpqE)1/l ⇒ dE/dx = K[E(x)]−1/l

⇒ dE (l+1)/l/dx =

l + 1

l

K,

(2.22)

where the constant K = (Nt/)q(l−1)/l(Jp/Nvµp)1/l. The general solution is in the form

E (l+1)/l(x) = ax + b.

Applying the ohmic boundary condition (in eq. (2.7)), we ﬁnd the electric ﬁeld

E(x) = al/(l+1)xl/(l+1),

(2.23)

(2.24)

and the electric potential

ψ(x) = −

(cid:90)

E(x)dx = −al/(l+1)

(cid:18) l + 1

(cid:19)

2l + 1

x(2l+1)/(l+1) + c.

(2.25)

Again, this electric potential is valid when the band bending eﬀect near contacts is negligible.
Applying the boundary condition from the device setup: ψ(x = 0) = 0 and ψ(x = L) = −V .

16

The integral constants can be determined as

(cid:18)2l + 1

(cid:19)(l+1)/l V (l+1)/l

l + 1

L(2l+1)/l

a =

and c = 0.

(2.26)

Substituting it into the solution in equation (2.24) and the combined diﬀerential equation

(2.22), we have the hole current in the presence of a ﬁeld-independent exponentially dis-

tributed traps,

Jp,T L = q1−lµpNv

(cid:18) 

(cid:19)l(cid:18)2l + 1

(cid:19)l+1 V l+1

l

Nt

l + 1

l + 1

L2l+1

.

(2.27)

Unlike the Mott-Gurney equation, the trap-limited SCLC gives a steeper slope of (l + 1)
in the LogJp − LogV plot, and hence the characteristic trap temperature/energy can be

determined by slope ﬁtting [28,45]. Often this characteristic slope can be observed in organic

semiconductors with low crystallinity, where traps dominate. The increasing of the applied

voltage eventually allows most trap sites to be ﬁlled. The transition from trap-limited SCLC

to trap-free SCLC occurs at voltage

(cid:34)

(cid:18)

VT F L =

N−1

v

9
8

Nt

l + 1

l

(cid:19)(l+1)(cid:35)1/(l−1)

(cid:19)l(cid:18) l + 1

2l + 1

L2.

q


(2.28)

The corresponding equations for the electric ﬁeld, the electric potential, and the hole

density can be derived as follow,

(cid:16) x

(cid:17)l/(l+1)

E(x)
EΩ

=

2l + 1
l + 1

L

(cid:16) x

(cid:17)(2l+1)/(l+1)

L

,

pf (x)

¯pf

=

(cid:16) x

(cid:17)−l/(l+1)

1

l + 1

L

, (2.29)

,

ψ(x)
−V

=

where EΩ is an ohmic ﬁeld deﬁned in (2.15), and ¯pf is the spatial average of the hole carrier

17

density,

(cid:90) L

0

¯pf = (1/L)

(cid:18) 

l

qNt

l + 1

(cid:19)l V l

L2l

,

2l + 1
l + 1

(2.30)

pf (x)dx = Nv(l + 1)

The simpliﬁed theories of both trap-free SCLC and trap-limited SCLC in the presence

of ﬁeld-independent exponentially distributed traps are summarized in table 2.1.

2.4 The Simpliﬁed Model of SCLC: Tip-Plane (c-AFM)

Geometry

The conductive atomic force microscope (c-AFM) is a type of atomic force microscope

equipped with a conductive cantilever tip. In addition to the ability to acquire topologi-

cal properties, c-AFM is capable of obtaining the local electrical transport properties of the

sample. During its operation, the conductive tip scans over the surface of the sample to

construct an image of its morphology. In the mean time, a voltage is applied between its

tip and the sample, and the induced current is measured to produce a high-resolution cur-

rent map. Due to its unique ability of simultaneous investigation of nanoscale morphology

and electronic properties, c-AFM has been used extensively to study organic semiconduc-

tors in the past decade [29, 30, 31, 32, 33, 34, 40, 41, 42, 46]. One of these eﬀorts is to probe

the local charge transport by measuring the local current-voltage characteristic of organic

semiconductors in the SCLC regime, and quantitatively determining the local charge mobil-

ity [29,31,40,41,42,46]. For studies of the local current-voltage response, the tip is positioned

at a speciﬁc location on the surface. The working range of the applied voltage is typically
between ±10V , and the current can be of the order of picoamperes when being used to study

charge transport in organic semiconductors.

18

When getting beyond the conventional planar device geometry, the theory of SCLC

meets a new challenge.
In the early c-AFM work, the eﬀective hole mobility in poly(3-
hexylthiophene) (P3HT) was found to be 1.4 × 10−2cm2V −1s−1 [31], which is much higher

than what was previously observed in a planar diode conﬁguration [27, 47, 48]. This discrep-

ancy brought attention to the importance of the tip geometry, especially current spread-

ing [29, 35], to which the existing analytical solutions for the conventional device geometry

can no longer be applied. Reid et al. [32] introduced a semi-empirical formula that includes a

tip geometry correction on the traditional Mott-Gurney law equation (2.13), which is written

as

(cid:19)1.6±0.1

(cid:18) L

d

V 2
L3 δ

Jp,c−AF M = αµp

(2.31)

where α and δ are two empirical dimensionless parameters, d is the diameter of the tip/sub-

strate contact area, and L is the ﬁlm thickness.

The authors [32] observed that the electric current density measured by c-AFM was

much higher than that in the conventional diode conﬁguration. The scaling exponent was
found to be 1.6 ± 0.1, which reﬂects the current enhancement due to the altered tip/elec-
trode geometry. The dimensionless prefactor α assumes a value of ∼ 8.2, as determined

from the numerical calculation of the standard drift-diﬀusion model. Although the model

includes charge diﬀusion, its eﬀect is insigniﬁcant since the prefactor is achieved at quite

high applied voltage (Va = 10V ). Another scaling parameter introduced in equation (2.31)

is δ. This parameter scales the numerical solution to ﬁt the experimental results of c-AFM

measurements, and its empirical value depends on the ratio of the tip diameter and the
sample thickness. For 0.01 ≤ d/L ≤ 2, this scaling parameter δ was found to have a value

19

of 7.8± 1. Later, the eﬀective mobilities yielded by using the semi-empirical equation (2.31)

on the c-AFM current-voltage data measured on the ﬁlms of P3HT, poly[2-methoxy-5-(3’,7’-

dimethyloctyloxy)-1,4-phenylene vinylene] (MDMO-PPV), and poly-(9,9-dioctylﬂuorene-co-

bis-N,N-(4-butylphenyl)-bis-N,N-phenyl-1,4-phenylenediamine (PFB) were shown to be con-

sistent with the one reported from analyses of diode device measurements. With growing

popularity of c-AFM technique in the research community, the formula (2.31), is utilized

in many recent studies, evaluating values of charge carrier mobility in organic semiconduc-

tors [33, 34, 49].

Though the semi-empirical equation (2.31) successfully corrects the Mott-Gurney law by
taking into account the tip geometry, the formula still follows the square law (Jp ∝ V 2)

which can only be observed in low trap devices. In the next chapter, we will extend the

work by Reid et al. to cover various bulk properties, including trap states, spatially varying

mobility, and morphological eﬀects, under relatively low applied voltages such that diﬀusion

cannot be neglected.

20

Table 2.1: Unipolar space-charge-limited current (SCLC) for hole injection of a planar device.

No traps

V 2
L3

L

(cid:16) x
(cid:17)1/2
(cid:17)3/2
(cid:16) x
(cid:17)1/2
(cid:16)L

L

1
2

x

µp

Jp =

9
8
E(x) =

3
V
2
L
ψ(x) = −V

pf (x) = ¯pf

¯pf =

3
2

V
qL2

Hole current density

Electric ﬁeld

Electric potential

Hole carrier density

Spatial avg hole density

Field-independent exponential distributed traps
Gt(E) = (Nt/kbTt)exp(−E/kbTt), l = Tt/T

(cid:19)l(cid:18)2l + 1

(cid:19)l+1 V l+1

l + 1

L2l+1

l

Nt

l + 1

(cid:18) 
(cid:16) x
(cid:17)l/(l+1)
(cid:17)(2l+1)/(l+1)
(cid:17)−l/(l+1)
(cid:16) x
(cid:18) 

L

L

l

qNt

l + 1

(cid:19)l V l

L2l

2l + 1
l + 1

Jp = q1−lµpNv
(cid:16) x

2l + 1
l + 1

E(x) =

ψ(x) = −V

V
L

L
1

l + 1

pf (x) = ¯pf

¯pf = Nv(l + 1)

21

Chapter 3

The Drift-Diﬀusion SCLC Transport

Models - New Method

In this chapter, we begin by describing the drift-diﬀusion (DD) SCLC model for hole-only

devices. Our theoretical model improves the simpliﬁed SCLC model, as previously described

in chapter 2, in two aspects: (i) the model will incorporate both drift and diﬀusion transport

mechanisms, and (ii) the model will enable the treatment of inhomogeneous systems including

c-AFM tip geometry, nanoscale morphology, as well as spatially-varying trap distributions.

As a result, this mathematical problem becomes very challenging, and impossible to solve

analytically. We will explain the numerical techniques used to achieve numerical solutions

of the DD-SCLC model. The numerical tools are written in Fortran 90 for one- two- and

three- dimensional models. All models will be veriﬁed to be consistent with each other, and

validated by reproducing the result of the simpliﬁed SCLC model in the drift-dominated

regime.

3.1 Modeling the Hole Injection System

In ﬁgure 3.1, the typical operation of hole-only devices is schematically illustrated in a

simple metal-insulator-metal (MIM) structure with ohmic contact at the electrodes. When

22

Figure 3.1: The schematic representation of hole-only devices The path of hole injection and
hole transport are indicated

a positive bias Va is applied, holes are injected from the anode, travel through the device,

and are extracted at the cathode. This is the case of unipolar SCLC transport, where one

carrier type is dominant. In this particular case, the charge carrier recombination can be

ignored, and the hole transport is solely determined by the bulk properties of the organic

semiconductors.

3.1.1 Governing Equations

As discussed previously in chapter 2, the bulk transport of holes in the SCLC regime can

be described by three governing equations, including Poisson’s equation (2.1), the conti-

nuity equation (2.2), and the drift-diﬀusion equation (2.3). The early theoretical work is

vastly simpliﬁed by ignoring the contribution of charge diﬀusion and only considering the

drift driven transport, such that the current-voltage characteristic of the hole-only device

can be derived analytically. The theoretical model presented in this chapter, however, will

include both types of charge transport, drift and diﬀusion. Furthermore, the material in-

23

homogeneity will also be studied by including the spatially-varying dielectric constant and

the spatially-varying hole mobility. The numerical calculations will be developed for one-

two- and three- dimensional models. The one-dimensional (1-D) model will use a Cartesian

coordinate system, the two-dimensional (2-D) model will use a cylindrical coordinate sys-

tem with cylindrical symmetry, and the three-dimensional (3-D) model will use a Cartesian

coordinate system. Hence, all three governing equations are rewritten in more general form

as follows

Poisson’s equation :

Diﬀerential form,

Integral form,

(cid:73)

S=∂V

(cid:126)∇ · ((cid:126)E) = q(pf + pt)
(cid:90)

(cid:18)(cid:90)

V

(cid:126)E · d (cid:126)A = q

pf dV +

V

(3.1a)

(cid:19)

ptdV

= Qp

(3.1b)

Continuity equation at steady state :

Diﬀerential form,

Integral form,

(cid:126)∇ · (cid:126)Jp = 0

(cid:73)

S=∂V

(cid:126)Jp · d (cid:126)A = 0

Drift-diﬀusion equation for holes :

(cid:126)Jp = qµppf

(cid:126)E − qDp (cid:126)∇pf

(3.2a)

(3.2b)

(3.3)

Poisson’s equation (3.1) describes the relation between local electric ﬁeld and local space

24

charge. Here, we present both diﬀerential and integral forms of the equation which will later

be used in various discretization schemes. The equation is initially written in terms of the

electric ﬁeld. Later the numerical calculations will be performed to determine the electric
potential ψ((cid:126)r), relating to the electric ﬁeld through the relation (cid:126)E = − (cid:126)∇ψ.

Due to the energetic disorder commonly found in organic semiconductors, it is important

to include the trap states into the model. The total holes Qp((cid:126)r) shall combine both trapped

holes and free holes. Although these trapped holes do not contribute to the charge transport

of the device, they do inﬂuence the local electric ﬁeld; together with free holes. The trap

DOSs is assumed to be exponentially distributed,

pt((cid:126)r) = Nt(pf /Nv)T /TtS((cid:126)r)

(3.4)

where S((cid:126)r) is the spatial distribution of traps. If trap density in the material is uniform,
S((cid:126)r) becomes unity. Nv is the eﬀective DOSs in the HOMO, Nt is the total trap density

and and Tt is trap temperature.

While the trapped holes are localized, the free holes travel through the bulk material

as is described by the well-known drift-diﬀusion model, as expressed in eq.(3.3). Unlike

the simpliﬁed theory of SCLC, we conjecture that both drift and diﬀusion mechanisms are

locally signiﬁcant in the low voltage regime. While drift current is driven by the gradient

of the electric potential (local electric ﬁeld), diﬀusion current is driven by the gradient of

the free hole density. The hole diﬀusion coeﬃcient Dp reﬂects the local capability of charge

diﬀusion, which obeys the Einstein relation

Dp((cid:126)r) = µp((cid:126)r)kbT /q

(3.5)

25

where kb is Boltzmann’s constant and q is the fundamental charge. Note that, per eq.(3.5),

the expression of the diﬀusion coeﬃcient has been conﬁrmed to be valid in the regime of low

electric ﬁeld [50,51]. It is also important to note that, due to the fact that the generation and

recombination of charges do not occur away from the electrodes, the total hole current reaches

the steady state and becomes conserved, as described by the Continuity equation(3.2).

3.1.2 Electrode Boundary Conditions

In order for the set of governing equations (3.1)-(3.3) to have a unique solution, boundary

conditions are required. By increasing the number of degrees of freedom ranging from the

one-, two and three- dimensional models, we can simulate more complex geometry and the

suitable boundary conditions inevitably vary from model to model. Therefore, the complete

set of boundary conditions for each model will later be explained in detail. However, the

boundary condition at each electrode is similar to that of the simpliﬁed SCLC model deﬁned

in chapter 2.

Electrode boundary conditions :

ψ = Va and pf = Nv

at the hole injecting electrode (anode)

(3.6a)

ψ = 0 and pf = Nv

at the hole extracting electrode (cathode)

(3.6b)

To inject holes into an organic semiconductor device, the positive bias Va is applied at the

injecting electrode. The injected holes then transport through the bulk material, and leave

the device at the extracting electrode which is grounded. To be consistent with the derivation

of the Mott-Gurney equation, both electrodes are assumed to be an ohmic contact, lining

26

up with the HOMO level of the organic semiconductor where the hole density is enormous
(Nv ∼ 1021cm−3). These boundary conditions at the electrode are similar to those used by

Koster et al. [52].

3.1.3 Scaling Factors

In the previous section, we have introduced the governing equations and the boundary con-

ditions at the electrodes. By performing numerical calculations (to be explained in section

3.2), we will obtain the numerical solutions for electric potential and free hole carrier den-

sity, which are diﬀerent by many orders of magnitude. Therefore, these variables need to

be scaled in order to increase the stability and improve the convergence of the calculations.

The electric potential and the electric ﬁeld are scaled by the thermal potential Vt = kbT /q,

while the hole carrier density is scaled by the eﬀective DOSs in the HOMO (Nv). With these

scaling factors, the three governing equations and the electrode boundary conditions can be

rewritten as

Poisson’s equation :

Diﬀerential form,

Integral form,

(cid:73)

f + ftS · p(cid:48)

f

rt)

(cid:126)∇ · ((cid:126)E(cid:48)) = (qNv/Vt)(p(cid:48)
(cid:90)

(cid:126)E(cid:48) · d (cid:126)A =

S=∂V

qNv
Vt

V

f + ftS · p(cid:48)
(p(cid:48)

f

rt)dV

where trap parameters ft = Nt/Nv and rt = T /Tt < 1.

Continuity equation at steady state :

Diﬀerential form,

(cid:126)∇ · (cid:126)Jp = 0

27

(3.7a)

(3.7b)

(3.8a)

Integral form,

(cid:73)

S=∂V

(cid:126)Jp · d (cid:126)A = 0

Drift-diﬀusion equation for holes :

(cid:126)Jp = qµpNvVt(p(cid:48)

f

(cid:126)E(cid:48) − (cid:126)∇p(cid:48)
f )

Electrode boundary conditions :

(3.8b)

(3.9)

ψ(cid:48) = Va/Vt and p(cid:48)
ψ(cid:48) = 0 and p(cid:48)

f = 1

f = 1

at the hole injecting electrode (anode)

(3.10a)

at the hole extracting electrode (cathode)

(3.10b)

In the following sections, we will develop and apply numerical methods for one-, two- and
three- dimensional DD-SCLC models, and prime((cid:48)) will be dropped to simplify the notation.

3.2 Numerical Method

First, consider the one-dimensional case, the x coordinate is divided into N mesh points,
separated by mesh size ∆xi = xi+1−xi as shown in Figure 3.2(a). Then PDEs are discretized

over the mesh points by using numerical derivative for their diﬀerential form, or numerical

integration for their integral form. These discrete equations are presented as a set of linear

equations with solutions as discrete functions. Unlike the analytical solutions, the values

of the discrete functions are only deﬁned at mesh points. Interpolation techniques are then

applied to approximate the values of these functions away from mesh points. Here, ﬁnite

diﬀerence methods (FDM) are used to approximate derivatives and integrations.

28

integration(cid:82) xi+1

xi

Figure 3.2: Numerical method: (a) illustration of spatial discretization in a one-dimensional
Cartesian coordinate system,(b) numerical derivative df /dx at the mid point xi, (c) numerical

f (x)dx around the mid point xi+1/2

3.2.1 Numerical Derivatives

The numerical approximations to derivatives are based on Taylor’s theorem. For the one-

dimensional case, the truncated Taylor’s series expansions of the function f (xi±1) around

x = xi are

f (xi+1) = f (xi) + ∆xi

+ O(∆x3
i )

(3.11a)

f (xi−1) = f (xi) − ∆xi−1

+ O(∆x3

i−1)

Subtracting these two equations gives

(3.11b)

(3.12)

(cid:12)(cid:12)(cid:12)(cid:12)x=xi
(cid:12)(cid:12)(cid:12)(cid:12)x=xi

+

∆x2
i

2

d2f
dx2

∆x2
i−1
2

d2f
dx2

+

df
dx

df
dx

(cid:12)(cid:12)(cid:12)(cid:12)x=xi
(cid:12)(cid:12)(cid:12)(cid:12)x=xi
(cid:12)(cid:12)(cid:12)(cid:12)x=xi

29

df
dx

≈ f (xi+1) − f (xi−1)
∆xi + ∆xi−1

This expression is the ﬁnite diﬀerence approximation of the ﬁrst-order derivative, as illus-

trated in ﬁgure 3.2(b). It is important to note that the numerical approximation converges
to the analytical derivative as (∆xi + ∆xi−1) → 0. The equation (3.12) is rewritten in terms

of the mid-point values as

(cid:12)(cid:12)(cid:12)(cid:12)x=xi

df
dx

≈ f (xi+1/2) − f (xi−1/2)
(∆xi + ∆xi−1)/2

(3.13)

In a similar way, we can approximate the higher-order derivatives, for example, the second

order derivative is given below

(cid:12)(cid:12)(cid:12)(cid:12)x=xi

d2f
dx2

≈

≈

(cid:12)(cid:12)(cid:12)(cid:12)x=xi+1/2

df
dx

(cid:12)(cid:12)(cid:12)(cid:12)x=xi−1/2

− df
dx

(∆xi + ∆xi−1)/2

f (xi+1) − f (xi)

− f (xi) − f (xi−1)

∆xi−1

∆xi

(∆xi + ∆xi−1)/2

(3.14)

When the domain is uniformly divided, such that ∆xi = ∆xi−1 = ∆x, the numerical approx-

imation of ﬁrst-order and second-order derivatives, as expressed in eq.(3.12) and eq.(3.14)

respectively, are simpliﬁed to

df
dx

(cid:12)(cid:12)(cid:12)(cid:12)x=xi
(cid:12)(cid:12)(cid:12)(cid:12)x=xi

≈ f (xi+1) − f (xi−1)

2∆x

d2f
dx2

≈ f (xi+1) + f (xi−1) − 2f (xi)

∆x2

(3.15)

(3.16)

To treat the inhomogeneity of organic materials, multi-dimensional calculations are re-

30

quired.

In three-dimensional Cartesian coordinate systems, the computational domain is

discretized using constant values of ∆x, ∆y and ∆z in the x, y and z directions respectively.

Each mesh point in the domain is indexed by (xi, yj, zk). In what follows, and unless oth-

erwise stated, the mesh points are assumed to be equally separated on each dimension with
constant distance ∆x = xi+1− xi, ∆y = yj+1− yj, and ∆z = zk+1− zk. Consider the scalar

function F (x, y, z), which is discretized as F (xi, yj, zk).

(cid:126)∇F (x, y, z)

= ˆx

∂F
∂x

+ ˆy

∂F
∂y

(cid:12)(cid:12)(cid:12)(cid:12)(xi,yj ,zk)

(cid:12)(cid:12)(cid:12)(cid:12)(xi,yj ,zk)

(cid:12)(cid:12)(cid:12)(cid:12)(xi,yj ,zk)

+ ˆz

∂F
∂z

(cid:12)(cid:12)(cid:12)(xi,yj ,zk)

≈ F (xi+1, yj, zk) − F (xi−1, yj, zk)
ˆx
F (xi, yj+1, zk) − F (xi, yj−1, zk)

2∆x

2∆y

F (xi, yj, zk+1) − F (xi, yj, zk−1)

+

+

ˆy

ˆz

(3.17)

(cid:12)(cid:12)(cid:12)(xi,yj ,zk)

∇2F (x, y, z)

=

∂2F
∂x2

2∆z

+

∂2F
∂y2

(cid:12)(cid:12)(cid:12)(cid:12)(xi,yj ,zk)

(cid:12)(cid:12)(cid:12)(cid:12)(xi,yj ,zk)

(cid:12)(cid:12)(cid:12)(cid:12)(xi,yj ,zk)

+

∂2F
∂z2

≈ F (xi+1, yj, zk) + F (xi−1, yj, zk) − 2F (xi, yj, zk)
F (xi, yj+1, zk) + F (xi, yj−1, zk) − 2F (xi, yj, zk)

∆x2

+

+

∆y2

F (xi, yj, zk+1) + F (xi, yj, zk−1) − 2F (xi, yj, zk)

∆z2

(3.18)

For simplicity, we will write the discrete function with the index notation, such that F (i, j, k) =

F (xi, yj, zk) and (i, j, k) = (xi, yj, zk). Equation (3.17) and equation (3.18) are rewritten as

31

(cid:12)(cid:12)(cid:12)(i,j,k)

(cid:126)∇F

= ˆx

∂F
∂x

(cid:12)(cid:12)(cid:12)(cid:12)(i,j,k)

+ ˆy

∂F
∂y

(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(i,j,k)

(cid:12)(cid:12)(cid:12)(cid:12)(i,j,k)

+ ˆz

∂F
∂z

≈ F (i + 1, j, k) − F (i − 1, j, k)
F (i, j, k + 1) − F (i, j, k − 1)

2∆x

ˆx +

F (i, j + 1, k) − F (i, j − 1, k)

2∆y

ˆy

+

2∆z

ˆz

(3.19)

(cid:12)(cid:12)(cid:12)(i,j,k)

∇2F (x, y, z)

(cid:12)(cid:12)(cid:12)(cid:12)(i,j,k)

(cid:12)(cid:12)(cid:12)(cid:12)(i,j,k)

+

∂2F
∂y2

(cid:12)(cid:12)(cid:12)(cid:12)(i,j,k)

+

∂2F
∂z2

=

∂2F
∂x2

≈ F (i + 1, j, k) + F (i − 1, j, k) − 2F (i, j, k)
F (i, j + 1, k) + F (i, j − 1, k) − 2F (i, j, k)

∆x2

+

+

∆y2

F (i, j, k + 1) + F (i, j, k − 1) − 2F (i, j, k)

∆z2

(3.20)

3.2.2 Numerical Integrations

We can also apply Taylor’s theorem to obtain the numerical approximation to integration.

Consider a one-dimensional integral of function f (x) from xi to to xi+1. To evaluate the

numerical integration, the function f (x) is substituted by its truncated Taylor expansions

32

around the mid-point x = xi+1/2.

(cid:90) xi+1

xi

f (x(cid:48))dx(cid:48) =

(cid:104)
f (xi+1/2) + (x(cid:48) − xi+1/2)

(cid:90) xi+1
xi
(x(cid:48) − xi+1/2)2

+

2

(cid:12)(cid:12)(cid:12)(cid:12)x(cid:48)=xi+1/2

d2f
dx2

df
dx

(cid:12)(cid:12)(cid:12)(cid:12)x(cid:48)=xi+1/2
(cid:21)

dx(cid:48)

+ O(∆x3
i )

≈ f (xi+1/2)∆xi

(3.21)

Equation (3.21) expresses the numerical integration, as illustrated in ﬁgure 3.2(c). Similarly,

the numerical integral of the function f (x) between the midpoint xi−1/2 and xi+1/2 is given

(cid:90) xi+1/2

xi−1/2

f (x(cid:48))dx(cid:48) =

(cid:90) xi+1/2
(cid:104)
f (xi) + (x(cid:48) − xi)
xi−1/2
(x(cid:48) − xi)2

(cid:12)(cid:12)(cid:12)(cid:12)x(cid:48)=xi

+

d2f
dx2

2
∆xi−1 + ∆xi

≈ f (xi)

2

(cid:12)(cid:12)(cid:12)(cid:12)x(cid:48)=xi
(cid:21)

df
dx

+ O(∆x3
i )

dx(cid:48)

(3.22)

We prefer this form of integral since the approximation only involves the value of the function

at the mesh point f (xi). Choosing the uniform mesh separation ∆xi−1 = ∆xi = ∆x, we

have

(cid:90) xi+1/2

xi−1/2

f (x(cid:48))dx(cid:48) ≈ f (xi)∆x

(3.23)

Notice that it is written in a form analogous to equation (3.21).

Often we need to perform two- and three- dimensional integrals. Consider the uniform

discretization in a three-dimensional domain, i.e., suppose that mesh points are located at

(xi, yj, zk), the distances between mesh points ∆x, ∆y, and ∆z along each direction are

33

constant.

(cid:90) xi+1/2

(cid:90) xi+1/2
(cid:90) yj+1/2

(cid:90) yj+1/2
(cid:90) zk+1/2

yj−1/2

xi−1/2

xi−1/2

yj−1/2

zk−1/2

F (x(cid:48), y(cid:48))dx(cid:48)dy(cid:48) ≈ F (xi, yj)∆x∆y

F (x(cid:48), y(cid:48), z(cid:48))dx(cid:48)dy(cid:48)dz(cid:48) ≈ F (xi, yj, zk)∆x∆y∆z

With the index notation, the numerical approximations to the one-, two- and three- dimen-

(3.24)

(3.25)

(3.26)

(3.27)

(3.28)

sional integrals are rewritten as

xi−1/2

(cid:90) xi+1/2
(cid:90) yj+1/2
(cid:90) zk+1/2

yj−1/2

(cid:90) xi+1/2
(cid:90) yj+1/2

xi−1/2

(cid:90) xi+1/2

f (x(cid:48))dx(cid:48) ≈ f (i)∆x

F (x(cid:48), y(cid:48))dx(cid:48)dy(cid:48) ≈ F (i, j)∆x∆y

xi−1/2

yj−1/2

zk−1/2

F (x(cid:48), y(cid:48), z(cid:48))dx(cid:48)dy(cid:48)dz(cid:48) ≈ F (i, j, k)∆x∆y∆z

3.2.3 Scharfetter-Gummel Discretization

Since the numerical solver of Poisson’s equation (3.7) can be achieved by following the

standard discretization discussed in section 3.2.1 and 3.2.2, one might expect to apply the

same discretization scheme to the rest of the governing equations, which are the continuity

equation (3.8) and the drift-diﬀusion equation (3.9). However, numerical instability limits

the eﬀectiveness of this approach.

Consider the one-dimensional case where the domain of interest is discretized into a set

of mesh points indexed by i = xi with equal spacing ∆x. The continuity equation (3.8) is

discretized by using the numerical formula (3.13), thus we need to know the value of the

hole current at the mid-point Jp(i + 1/2) evaluated from the drift-diﬀusion equation (3.9).

34

Following the standard numerical discretization discussed in section 3.2.1 and 3.2.2, we have

(cid:34)

Jp,x(i + 1/2) = qµpNvVt

pf (i + 1/2)

ψ(i) − ψ(i + 1)

∆x

− pf (i + 1) − pf (i)

∆x

(3.29)

(cid:35)

This is based on the assumption that the hole mobility is constant within the interval

[xi, xi+1]. The value of hole density at the mid-point is approximated by the arithmetic

average of the neighboring mesh points, such that pf (i + 1/2) ≈ (cid:2)pf (i) + pf (i + 1)(cid:3) /2.

Then, we have

(cid:34)

Jp,x(i + 1/2) = qµpNvVt

pf (i)

(cid:18) ψ(i) − ψ(i + 1)

+

(cid:19)

2∆x

+ pf (i + 1)

1
∆x

(cid:18) ψ(i) − ψ(i + 1)

2∆x

(cid:19)(cid:35)

− 1
∆x

(3.30)

To guarantee a smooth ﬂow of current between mesh points, the numerical calculation is

restricted by Reynolds number [53]. Therefore, for the stability of the calculation, we need to
have ψ(i+1)−ψ(i) (cid:54) 2 in the scaled units, which is equivalent to ψ(i+1)−ψ(i) (cid:54) 2kBT /q ≈

0.052V at room temperature (300K). Then the mesh interval needs to be suﬃciently small

for this condition to hold throughout the domain.

Scharfetter and Gummel have proposed the optimum solution to this problem [54]. With

their approach, two fundamental assumptions need to be made:

(i) hole current density Jp,x is constant between mesh points,

(ii) electric potential ψ is linear between mesh points, implying that electric ﬁeld Ex

is also constant between mesh points.

35

We can rewrite the hole current density at the mid-point as

(cid:18)

(cid:19)

Jp,x(i + 1/2) = qµpNvVt
⇒ dpf
dx

Ex(i + 1/2) pf (x) − dpf
dx
− Ex(i + 1/2) pf (x) = −Jp,x(i + 1/2)
qµpNvVt

(3.31)

where Ex(i + 1/2) = [ψ(i) − ψ(i + 1)] /∆x = −∆ψ(i + 1/2)/∆x. Clearly, equation (3.31) is

a ﬁrst-order diﬀerential equation that has general solution in the form

pf (x) = a + becx

(3.32)

Substituting this expression into equation (3.31), we have

[Ex(i + 1/2) − c] ecx =

Jp,x

qµpNvVt

− Ex(i + 1/2) a ⇒ 0.

(3.33)

The left-hand side of the equation is x-explicitly-dependent, while the right-hand side is not.

This equation is only valid when the values on both sides are zero, so the constants a and c

are then determined as

a =

Jp,x

qµpNvVtEx(i + 1/2)

and c = Ex(i + 1/2).

And, we now have

pf (x) =

Jp,x

qµpNvVtEx(i + 1/2)

+ b · eEx(i+1/2)x

(3.34)

(3.35)

36

Let’s consider the value of hole density at mesh points i and i + 1.

pf (i) =

Jp,x

qµpNvVtEx(i + 1/2)

+ b · eEx(i+1/2) xi

pf (i + 1) =

Jp,x

qµpNvVtEx(i + 1/2)

+ b · eEx(i+1/2) xi+1

Combining the two equations, and eliminating the constant b, we obtain,

Jp,x(i + 1/2) =

qµpNvVt

∆x

(cid:104)B(∆ψ(i + 1/2)) · pf (i) − B(−∆ψ(i + 1/2)) · pf (i + 1)

(3.36a)

(3.36b)

(cid:105)

,

(3.37)

where the Bernoulli function

B(x) =

x

ex − 1

with B(0) = 1.

(3.38)

Similarly, we can also derive

Jp,x(i − 1/2) =

qµpNvVt

∆x

(cid:104)B(∆ψ(i − 1/2)) · pf (i − 1) − B(−∆ψ(i − 1/2)) · pf (i)
(cid:105)

.

(3.39)

The discretization given by eq.(3.37) and eq.(3.39) is known as Scharfetter-Gummel dis-

cretization; and is an important result in device physics [53, 55, 56, 57, 58, 59, 60].

3.2.4 Numerical Calculation Scheme

Adopting the strategy of Koster et al. [52], a self-consistent scheme will be used to solve the

set of governing equations (3.7)-(3.9). A ﬂow diagram of the numerical calculation scheme

is shown in ﬁgure 3.3. To evaluate the current-voltage characteristic of the hole-only device,

37

we ﬁrst deﬁne an initial guess of the solutions, which are the electric potential ψo and free

hole carrier density po

f . We then employ the numerical technique to solve the discretized

Poisson’s equation, obtaining the new values of electric potential ψnew. Next, the discretized

continuity equation coupled with the discretized drift-diﬀusion equation is solved using the

same numerical technique and the updated values of the electric potential, yielding the new

values of free carrier density pnew

f

. The correction of the resulting free hole carrier density

determines if the simulation is converged, i.e., when the relative correction of free hole carrier

density is smaller than its tolerance δ. The tolerance used in this thesis is in the range of
10−12 − 10−15. The values of free hole carrier density are updated, and the process is

repeated until convergence is reached.

To improve the rate of convergence, we will use a weighted average of the two most recent

estimates to obtain the next estimate of the solution. In general, the updated numerical

solution F is expressed as

F = w · F n + (1 − w) · F n−1

(3.40)

where n is the iteration number and w is the weighting parameter (0 ≤ w ≤ 1). When w

is unity, we simply maintain the standard iterative scheme where the whole new value of

solution is passed to the next iteration of the calculation.

To solve this mathematical problem numerically, we developed a computer program writ-

ten in Fortran 90, which incorporates OpenMP to perform parallel computing. The calcula-

tions were carried out at the high performance computing center (HPCC) at Michigan State

University.

38

Figure 3.3: A ﬂow diagram used to calculate the current-voltage characteristic of hole-only
devices. The algorithm is based on the self-consistent scheme. We ﬁrst establish the device
setup with its morphology and parameter values, and then generate initial trial values for
electric potential ψo and free hole carrier po
f . Subsequently, Poisson’s equation is solved
numerically to ﬁnd the estimated values of potential ψn. Utilizing the new trial values
of potential ψ, the estimated values of free hole carrier density pn
f are evaluated from the
continuity equation and the drift-diﬀusion equation, leading to the new trial values pf . This
entire process is repeated until convergence is reached, when the relative correction of free
hole carrier density is less than the value of the tolerance δ.

39

3.3

1-D DD-SCLC Model: Homogeneous Planar Ge-

ometry

Figure 3.4: Schematic illustration of a planar hole-only device made of a thin-ﬁlm of organic
semiconductor. The arrows illustrate the hole current path responding to hole injection at
the top electrode. When the material is homogeneous, the vertical conductance is uniform.
A 1-D DD-SCLC model is often suﬃcient to describe the hole carrier transport in such a
device.

To model charge transport in an organic semiconductor device, it is common to assume

that the transport layer in a planar device is homogeneous, such that the dielectric con-

stant is uniform, the hole carrier mobility is isotropic, and trap sites are evenly distributed

throughout the device with constant trap density and trap energy. Therefore, only the verti-

cal conductivity along the transport direction from one electrode to the other is considered,

and the mathematical description is reduced to a one-dimensional problem. Shown in ﬁgure

3.4 is a schematic illustration of hole-only device that consists of an organic semiconductor

thin ﬁlm of thickness L fabricated in between conventional plane-parallel electrodes, and

the scaled governing equations speciﬁcally for this 1-D DD-SCLC model, including Poisson’s

equation, the continuity equation at steady state, and the drift-diﬀusion equation for holes,

are written as

40

d2ψ/dx2 = (qNv/Vt)(pf + ft · pf

rt)

dJp/dx = 0

Jp = −qµpNvVt(pf · dψ/dx + dpf /dx)

The boundary conditions at the electrodes (x = 0, L) are assigned as

ψ(x = 0) = Va/Vt

and

pf (x = 0) = 1

ψ(x = L) = 0

and

pf (x = L) = 1

(3.41)

(3.42)

(3.43)

(3.44a)

(3.44b)

It is worth emphasizing that this is the ﬁrst step to extend the simpliﬁed SCLC theory to

include the mechanism of charge diﬀusion.

3.3.1 Discretization

In ﬁgure 3.5, the one-dimensional device previously shown in ﬁgure 3.4 is uniformly dis-

cretized. The one-dimensional domain of length L is replaced by the total N+2 mesh

points indexed by i = 0, 1, 2, .., N + 1 with a constant mesh spacing h = L/(N + 1) or
xi = i · L/(N + 1).

Figure 3.5: Illustration of the spatial discretization in 1-dimensional Cartesian coordinate
system. The device thickness L is discretized into the grid of N + 2 mesh points. Each mesh
point is equally separated by ∆x = L/(N + 1).

41

3.3.2 Potential Solver

On a uniform mesh, Poisson’s equation (3.41) can be discretized by using the numerical

approximation of the second-order derivative (3.20)

−ψ(i − 1) + 2ψ(i) − ψ(i + 1) = h2(qNv/Vt)(cid:2)pf (i) + ft · pf (i)rt(cid:3)

(3.45)

where i = 1, 2, 3, .., N . We have the boundary conditions for the electric potentials as

ψ(0) = Va/Vt and ψ(N + 1) = 0

(3.46)

The mathematical treatment is straightforward. The N equations derived from equation

(3.45) together with the two boundary conditions (3.46) result in a linear system of N

equations with N unknown variables, which can be represented by



0

0

0

0
2 −1

2 −1
−1
0 −1
...

0
2 −1
. . .

0
. . .
. . .
0 −1

···
···
···

···
···
···

0

0

0
...

ψ(1)

ψ(2)

ψ(3)
...

2 −1

0
2 −1

0 −1

ψ(N − 2)
ψ(N − 1)









=

b(1) + ψ(0)

b(2)

b(3)
...

b(N − 2)
b(N − 1)



(3.47)

where b(i) = h2(qNv/Vt)(cid:2)pf (i) + ft · pf (i)rt(cid:3). Equation (3.47) can also be written in the

ψ(N )

b(N ) + ψ(N + 1)

0 −1

2

vector form

A (cid:126)ψ = (cid:126)b

42

(3.48)

where the vector (cid:126)ψ contains the unknown values of the electric potential, and the vector (cid:126)b
contains the constant b(i) and the boundary conditions ψ(0) and ψ(N +1). The N×N matrix
A contains the coeﬃcients of the variables, which are constant values and predominantly zero.

Therefore, we employ sparse matrix storing techniques to reduce the memory usage while

performing the numerical calculation.

There are two general schemes to solve the linear systems shown in equation (3.48)

numerically: (i) direct elimination method and (ii) iterative method [61, 62]. In this case,
the coeﬃcient matrix A is a symmetric tridiagonal matrix. The inversion of A can thus be

calculated by using the explicit inversion formula suggested by R. A. Usmani [63, 64]. The

solution may be found by multiplying both sides of equation (3.48) by the matrix inverse
A−1.

(cid:126)ψ = A−1(cid:126)b

(3.49)

3.3.3 Hole Density Solver

On a uniform mesh, the continuity equation (3.42) can be discretized by using the numerical

approximation of the ﬁrst-order derivative (3.19)

(cid:12)(cid:12)(cid:12)(cid:12)i

dJp
dx

≈ Jp(i + 1/2) − Jp(i − 1/2)

h

= 0.

(3.50)

Substituting the Scharfetter-Gummel discretization for the drift-diﬀusion equation (3.37)-

(3.39), we have

−B(∆ψ(i − 1/2)) · pf (i − 1) +

(cid:104)B(−∆ψ(i − 1/2)) + B(∆ψ(i + 1/2))

(cid:105) · pf (i)

−B(−∆ψ(i + 1/2)) · pf (i + 1) = 0 (3.51)

43

where i = 1, 2, 3, .., N and B is the Bernoulli function. The boundary conditions for the free

hole density are

pf (0) = 1 and pf (N + 1) = 1.

(3.52)

Following the similar mathematics, equation (3.51) may be written as

A (cid:126)pf = (cid:126)b

(3.53)

where the vector (cid:126)pf contains the unknown values of the free hole carrier density, the vector
(cid:126)b contains the boundary conditions ((cid:126)b = {B(∆ψ(i − 1/2)) · pf (0), 0, ..., 0,B(−∆ψ(i + 1/2)) ·
pf (N + 1)}) , and the matrix A contains the coeﬃcients of the variables, involving Bernoulli
functions. Again, the coeﬃcient matrix A is a N × N matrix with predominantly zero

elements, and is stored as a sparse matrix.

The matrix A is again a tridiagonal matrix, which can be inversed using the formula

proposed by Usmani [63, 64]. Subsequently, we ﬁnd the solution as

(cid:126)pf = A−1(cid:126)b

(3.54)

3.3.4 Veriﬁcation

We have developed a program written in Fortran 90 to carry out the potential solver and

hole density solver simultaneously by following the algorithm of the self-consistent scheme

shown in ﬁgure 3.3. General guidelines for design and veriﬁcation of the model include the

following considerations :

(i) minimize the error from numerical calculations

44

(ii) verify the theoretical modeling

First of all, the result of numerical solutions are approximate. The two major sources of

errors generally found in numerical calculations are round-oﬀ error and discretization error.

Regardless of the methodology employed, the calculation usually has the round-oﬀ error,

arising from the approximation of real numbers by the standard computer representation of

the ﬂoating point numbers with ﬁnite digits. In our model, most parameters and variables

are real. They are represented by double precision values in the program to minimize the

corresponding round-oﬀ error. Discretization error, on the other hand, comes from the

numerical approximations of the equations, which are usually proportional to a power of the

mesh interval.

The eﬀect of the discretization error can be seen from the numerical solution. To demon-

strate this eﬀect, we perform the numerical calculation to model hole carrier transport in the

hole-only devices of 80nm thickness. Following the discretization scheme shown in ﬁgure 3.5,

we gradually increase the mesh number, relating to the mesh interval through h = L/(N +1),

where N is the number of the unknown mesh points. The current density as simulated from

our 1-D DD-SCLC model with varying numbers of mesh points is shown in ﬁgure 3.6. It ap-

proaches the saturated regime with respect to the increasing number of mesh points, whereas

the discretization error is minimized. The calculation with a ﬁne discretization requires a

signiﬁcant amount of computational time and resources, especially in the higher dimensional

models. A good compromise between the computational eﬀort and the computational accu-

racy was found to be at 2nm mesh interval, which results in approximately a maximum of

7% diﬀerence compared to the converged numerical solution of the system with 0.1nm mesh

interval. Additionally, the 2nm mesh interval is an appropriate length scale for organic

semiconductor devices consisting of nanoscale morphologies. Lowering the mesh interval

45

ever further only improves the numerical solutions quantitatively owning to the decrease

of numerical errors, but does not vary the qualitative behaviors of devices. Therefore, the

standard 2nm mesh interval will be used in the subsequent calculation.

Next, we would like to verify the capability of our model to capture the physical behavior

of charge injection and transport in hole-only devices. According to the theory of SCLC, the

injecting electrode does not limit the hole current in the devices and consequently is assumed

to be ohmic, i.e., an inﬁnite source of charge carriers. Practically this implies a large charge

carrier density at the injecting electrode. In our DD-SCLC model, the hole carrier density

at the injecting contact is assumed to be equal to the eﬀective DOSs in the HOMO (Nv), as

expressed in eq.(3.6). The value of this parameter is determined from the regime where the

current density is slowly varying over the value of Nv as shown in ﬁgure 3.7. A very large

value of Nv, resulting in a steep gradient of charge density around the electrode, poses a

challenge to the calculation convergence. Empirically, to have an ohmic injecting electrode

yet ensure a well converged model, calculations will be performed with the value of Nv as
1.25 × 1021cm−3, which is more than suﬃcient for comparison with experiments.

To further verify the physical interpretation of our 1-D DD-SCLC model, we would like

to demonstrate the eﬀect of charge diﬀusion which is omitted from the simpliﬁed SCLC

theory (discussed in chapter 2). Again, the calculation is performed in the absence of traps
with a low hole mobility of µp = 10−4cm2V −1s−1. Figure 3.8 shows the simulated hole
current density (Jp) as a function of applied voltage (Va). By ﬁtting to a power law J ∝ V m

(solid line), the current-voltage characteristic exhibit two distinct regimes.

In the high

voltage regime (Va > 10) where drift current is dominant, the ﬁt of the current-voltage

relation yields the exponent m = 2, agreeing with the number predicted by the Mott-Gurney

equation (2.13). On the other hand, in the low voltage regime (Va < 10), the exponent is

46

Figure 3.6: Mesh size Veriﬁcation – Convergence of the simulated hole current density with
respect to the number of mesh points, when the applied voltage Va = 0.1, 1, 10V . Hole
current density is calculated from the 1-D DD-SCLC model of the trap-free (ft = 0) hole-only
devices of 80nm thickness (Nv = 1.25 × 10−21cm−3, µp = 10−4cm2V −1s−1). Discretizing
the system into 41 mesh points with the constant 2nm mesh interval provides a well converged
solution, and results in a maximum 7% error compared to the result from 801 mesh points
(0.1nm mesh interval).

47

Number of Mesh Points (N+1)60000700008000090000100000Current Density (A/m2)Va = 10 V0501001502005060708090Va = 0.1V100012001400160018002000Va =   1 VFigure 3.7: Nv Veriﬁcation – Convergence of the simulated hole current density with respect
to the eﬀective density of states in HOMO (Nv), when the applied voltage Va = 0.1, 1, 10V .
Hole current density is calculated from the 1-D DD-SCLC model of the trap-free (ft = 0)
hole-only devices of 80nm thickness (µp = 10−4cm2V −1s−1). The Nv value of 1.25 ×
1021cm−3 will be used in all following calculations.

48

The Effective Density of States Nv (cm-3)100100010000100000Current Density (A/m2)Va = 10 V101410151016101710181019102010211022102310241025110100Va = 0.1V101001000Va =   1 VFigure 3.8: Mott-Gurney convergence – The theoretical current-voltage characteristic cal-
culated from the 1-D DD-SCLC model of the trap-free (ft = 0) hole-only device obeys the
power law J ∝ V m, demonstrating the transition from the diﬀusion aﬀected regime of the
ﬁtting exponent m = 1.52 to the Mott-Gurney regime of the ﬁtting exponent m = 2 with
increasing applied voltage.

49

110100Applied Voltage (V)102103104105106107Current Density (A/m2)m = 2.00m = 1.52lowered due to the contribution of charge diﬀusion. The current-voltage characteristic in the

drift-dominated regime is clearly consistent to the prediction of the simpliﬁed SCLC model,

validating the 1-D DD-SCLC model and its numerical approach.

3.4

2-D DD-SCLC Model: Semi-Homogeneous Tip-

Plane Geometry

We have found that a 1-D DD-SCLC model is often suﬃcient to simulate the hole transport

in the case of planar organic semiconductor devices. However, the organic semiconductor ﬁlm

is often composed of both amorphous and crystalline phases. Many studies have reported

the exceptionally high conductance observed along the crystallized nanoﬁbers of the organic

semiconductor, resulting in an anisotropic mobility of hole carrier measured in thin ﬁlms

[34, 65, 66, 67]. Since the planar geometry demonstrates only the vertical conductance, one

way to observe the eﬀect of this lateral conductance is by replacing the top planar electrode by

a tip electrode. The so-called tip-plane geometry is analogous to the electrode conﬁguration

of c-AFM measurements commonly used to study organic semiconductors.

In order to simulate the charge transport in a hole injection system of tip-plane geometry,

the calculation is performed in the cylindrical coordinate system with cylindrical symmetry.

The tip electrode is then approximated as the circular contact of radius (Rtip) as shown in

ﬁgure 3.9. When this tip radius (Rtip) is equal to the device radius (R), the device setup

reverts to the planar geometry. The governing equations for the 2-D DD-SCLC model in the

50

Figure 3.9: Schematic illustration of a three-dimensional hole-only device of tip-plane ge-
ometry. The numerical calculation is performed using a cylindrical coordinate system with
the cylindrical symmetry, reducing the calculation domain to the two-dimensional system as
shown. Hole transport is then described by a 2-D DD-SCLC model.

cylindrical coordinate system are written as

(cid:73)

S=∂V

(ρ, z)(cid:126)E(ρ, z) · d (cid:126)A =

qNv
Vt

(cid:73)

(cid:0)pf (ρ, z) + S(ρ, z)ft · pf (ρ, z)rt(cid:1) dV

(cid:90)

V

(3.55)

(3.56)

(3.57)

S=∂V

Jp,ξ(ρ, z) = qµp,ξ(ρ, z)NvVt

(cid:126)Jp(ρ, z) · d (cid:126)A = 0

(cid:0)pf (ρ, z)Eξ(ρ, z) − ∇ξpf (ρ, z)(cid:1)

where (cid:126)Jp(ρ, z) = (cid:80)

ξ

Jp,ξ(ρ, z) ˆξ and ξ = {ρ, z}. Poisson’s equation (3.55) and the conti-

nuity equation (3.56) are presented in the integral form using the divergence theorem. In

contrast, the drift-diﬀusion equation (3.57) is shown in the diﬀerential form. The internal

inhomogeneity of the material is included in the model by implementing the spatially varying

parameters of dielectric constant (ρ, z) = r(ρ, z)o, hole mobility (cid:126)µp(ρ, z) =(cid:80)

µp,ξ(ρ, z) ˆξ,

and trap distribution S(ρ, z), where o is the permittivity of vacuum.

ξ

In the two-dimensional conﬁguration of ﬁgure 3.9, the organic semiconductor is grounded

51

on the bottom electrode, driven by the top electrode with the potential Va, and bounded

elsewhere by an insulator. Similar to the 1-D DD-SCLC model, the boundary conditions at

the tip-planar electrodes are deﬁned as

ψ(ρ, z) = Va/Vt and pf (ρ, z) = 1;

ψ(ρ, z) = 0 and pf (ρ, z) = 1;

where ρ ∈ [0, Rtip], z = L
where ρ ∈ [0, R], z = 0

(3.58a)

(3.58b)

Due to the absence of hole transport through the insulator, the normal component of hole

current vanishes at the interface. The insulating boundary conditions for potential and hole

density are assigned as

ˆn · ∇ψ = 0 and ˆn · ∇pf = 0

(3.59)

where ˆn is the unit vector normal to the insulating interface. Consider the physical meaning

of the insulating boundary conditions. The zero-gradient of electric potential and current

density can be interpreted as ‘no electric ﬂux’ and ‘no current ﬂux’ through the interface,

which is the Neumann boundary condition for a fully insulating boundary.

3.4.1 Discretization

The cylindrical domain is discretized as shown in ﬁgure 3.10. The inhomogeneity of the

material is also discretized. Each discrete volume, which later is reduced to a discrete area

according to the cylindrical symmetry, has a uniform dielectric constant, hole mobility and

trap parameters. Mesh points are then deﬁned on the edge of each discrete volume/area.

Therefore, the two-dimensional domain of radius R and thickness L shown in ﬁgure 3.9 is
replaced by discrete areas of equal size ∆ρ × ∆z and a set of mesh points (ρi, zk); i =

52

Figure 3.10: Schematic illustration of the spatial discretization in a cylindrical coordinate
system. With cylindrical symmetry, the two-dimensional device is divided into the sub-area
of size ∆ρ × ∆z with mesh points located at each corner. Each sub-area carries the single
values of dielectric constant, hole mobility, and trap parameters. While it is homogeneous
within sub-area, each sub-area can have diﬀerent parameter values, producing the inhomo-
geneity in the system.

Figure 3.11: Schematic illustration of Gaussian Surface of (a) axis cell for the centered mesh
at z-axis (i = 0), and (b) regular cell for the centered mesh elsewhere

{0, 1, 2, .., N} and k = {0, 1, 2, .., M}, where ρi = i · R/N and zk = k · L/M .

In order to solve the governing equations in the integral form, we deﬁne a Gaussian surface

(SG) bounding a ﬁnite volume surrounding each mesh point over the small angle ∆ϕ. Two

kinds of Gaussian surfaces are chosen as demonstrated in ﬁgure 3.11. The axis cell is the

wedge-like volume, corresponding to the mesh points at the z-axis (i = 0), which is adjacent

to three nearest neighbors listed in ﬁgure 3.11(a). On the other hand, the regular cell is a
parallelepipedic volume, corresponding to the mesh points elsewhere (i = {1, 2, 3, .., N}). A

53

regular mesh point is surrounded by four nearest neighbors listed in ﬁgure 3.11(b).

Because Poisson’s equation and the continuity equation in the integral form evaluate the

ﬂux of electric ﬁeld and hole current density at the Gaussian surface, three fundamental

assumptions can be made:

(i) hole current density Jp is constant between mesh points,

(ii) electric potential ψ is linear between mesh points, implying that the electric ﬁeld E is

also constant between mesh points.

(iii) the normal component of the electric ﬁeld at the Gaussian surface is continuous across

the interface of two discrete volumes/areas.

The ﬁrst two assumptions are meant to be consistent with those made in the Scharfetter-

Gummel discretization (see section 3.2.3). The last assumption is based on the commonly-

used interface boundary condition of the electric ﬁeld, stating that the component of the

electric ﬁeld that is tangential to the interface of two materials is continuous. Recall that

the internal inhomogeneity of the device is generated by varying the electrical and morpho-

logical properties of the discrete volumes/areas. In our numerical approach, the locations

of mesh points are speciﬁcally chosen to be at the corners of the discrete units as shown

in ﬁgure 3.10, and the correspondent Gaussian surfaces are always perpendicular to the in-

terfaces between the units as seen in ﬁgure 3.19, reasoning the assumption (iii). With the

combination of assumption (ii) and (iii), our numerical calculations selectively treat the in-

ternal heterogeneity of the devices locally through the average of the dielectric constant and

mobility over the Gaussian surface, which will be clearly seen in the following sections.

54

3.4.2 Potential Solver

Poisson’s equation (3.55) in the integral form is discretized by using the numerical approx-

imation of the surface and volume integral, as expressed in eq.(3.27) and eq.(3.28) respec-

tively, on the chosen Gaussian surface as shown in ﬁgure 3.11. The axis cell is speciﬁcally

selected for the discretization of mesh points at the z-axis (i = 0), while the regular cell is
for mesh points elsewhere (i = {1, 2, 3, ..., N}). The discretized Poisson’s equation based on

two Gaussian surface are

(ρ, z) ρdϕdz · Eρ(i + 1

2 , k) −

(cid:90)(cid:90)

(ρ, z) ρdρdϕ · Ez(i, k − 1
2)

SG(k− 1
2 )

(ρ, z) ρdρdϕ · Ez(i, k + 1
2)

At z-axis (ρ = 0 or i = 0)

(cid:34)(cid:73)

(cid:35)

(cid:126)E · d (cid:126)A

SG

(cid:90)(cid:90)
(cid:90)(cid:90)

SG(i+ 1
2 )

+
SG(k+ 1
2 )

(cid:40)

4¯ρ(i + 1

=

i,k

=

−

+

∆z

¯z(i, k − 1

2, k) · [ψ(i, k) − ψ(i + 1, k)]
(cid:18)∆ρ
(cid:19)2
(cid:19)2
(cid:18)∆ρ
(cid:18)∆ρ
(cid:19)2

∆ϕ∆z ·(cid:104)

¯z(i, k + 1

∆z

(cid:41)

2) · [ψ(i, k − 1) − ψ(i, k)]
2) · [ψ(i, k) − ψ(i, k + 1)]

pf (i, k) + ft(i, k) · pf (i, k)rt(i,k)(cid:105)

· 1
8

∆ϕ∆z

(3.60)

≡ qNv
Vto

· 1
2

2

55

Otherwise (ρ = i · ∆ρ where i = {1, 2, 3, .., N})

(cid:34)(cid:73)

(cid:35)

SG

i,k

=

(cid:126)E · d (cid:126)A

= −

(ρ, z) ρdϕdz · Eρ(i − 1

2 , k) +

(ρ, z) ρdρdϕ · Ez(i, k − 1

2) +

(cid:40)

SG(i− 1
2 )
−
SG(k− 1
2 )
−
1 − 1
2i

(cid:90)(cid:90)
(cid:90)(cid:90)
(cid:19)
(cid:18)
(cid:19)
(cid:18)
(cid:18)∆ρ
(cid:19)2 · ¯z(i, k − 1
(cid:19)2 · ¯z(i, k + 1
(cid:18)∆ρ
· i(∆ρ)2∆z∆φ ·(cid:104)

· ¯ρ(i − 1

· ¯ρ(i + 1

+

1 +

−

+

1
2i

∆z

∆z

(cid:90)(cid:90)
(cid:90)(cid:90)

SG(i+ 1
2 )

SG(k+ 1
2 )

(ρ, z) ρdϕdz · Eρ(i + 1

2, k)

(ρ, z) ρdρdϕ · Ez(i, k + 1
2)

2 , k) · [ψ(i − 1, k) − ψ(i, k)]
2, k) · [ψ(i, k) − ψ(i + 1, k)]
(cid:41)
2) · [ψ(i, k − 1) − ψ(i, k)]
2) · [ψ(i, k) − ψ(i, k + 1)]
pf (i, k) + ft(i, k) · pf (i, k)rt(i,k)(cid:105)

· i∆z∆φ

(3.61)

≡ qNv
Vto

·(cid:104)

where the local averages of the dielectric constant over each section of Gaussian surface

around the mesh point (i, k) are

(cid:105)

¯ρ(i ± 1

2 , k) =

¯z(i, k ± 1

2) =

1
2



(i ± 1

2 , k − 1

(i, k ± 1
2);
2 , k ± 1
(i − 1
2) ·

2) + (i ± 1
(cid:19)
(cid:18)

for i = 0,
1 − 1
4i

2 , k + 1
2)

for i = {0, 1, 2, .., N},(3.62a)

;

+ (i + 1

2 , k ± 1
2) ·

2

(cid:18)

1 +

1
4i

(cid:19)

;

(3.62b)

for i = {1, 2, 3, .., N},

and k = {0, 1, 2, ..., M}. When the domain is uniformly divided (∆ρ = ∆z = h), the

discretized equations (3.60)-(3.61) are simpliﬁed as

56

At z-axis (ρ = 0 or i = 0)

− 4¯ρ(i + 1

2 , k) · ψ(i + 1, k) − ¯z(i, k − 1
+ ¯(i, k) · ψ(i, k) =

2) · ψ(i, k − 1) − ¯z(i, k + 1
qNvh2
Vto

pf (i, k) + ft(i, k) · pf (i, k)rt(i,k)(cid:105)

2) · ψ(i, k + 1)

·(cid:104)

Otherwise (ρ = i · ∆ρ where i = {1, 2, 3, .., N})

(cid:18)

−

(cid:19)

· ¯ρ(i − 1

1 − 1
2i
− ¯z(i, k − 1

(cid:18)

· ¯ρ(i + 1

2 , k) · ψ(i + 1, k)

(cid:19)
pf (i, k) + ft(i, k) · pf (i, k)rt(i,k)(cid:105)
·(cid:104)

2) · ψ(i, k + 1) + ¯(i, k) · ψ(i, k)

1 +

1
2i

2 , k) · ψ(i − 1, k) −
2) · ψ(i, k − 1) − ¯z(i, k + 1
qNvh2 · i

=

Vto

(3.63)

(3.64)

where the local average of the dielectric constant at center point (i, k) are

4¯ρ(i + 1

1 − 1
2i

(cid:18)

¯(i, k) =

(cid:18)
(cid:19)
2 , k) + ¯z(i, k − 1
¯ρ(i − 1

2 , k) +

(cid:19)

2) + ¯z(i, k + 1
2);

for i = 0,

1 +

1
2i

¯ρ(i + 1

2 , k) + ¯z(i, k − 1

(3.65)
2) + ¯z(i, k + 1
2);

for i = {1, 2, 3, .., N},

and k = {0, 1, 2, ..., M}.

We have established the discretized Poisson’s equation (3.63) and (3.64), covering the
whole domain (∀i ∈ [0, N ], ∀k ∈ [0, M ]). In order to ﬁnd the numerical solutions of the

electric potential, all boundaries need to be assigned speciﬁc conditions, subject to the

device geometry shown in ﬁgure 3.9. The boundary conditions for the electric potential at

57

the electrodes are

Va/Vt

0

ψ(i, k) =

i ∈ [0, Rtip/h], k = M
i ∈ [0, R/h], k = 0

for

for

(3.66)

Additionally, boundaries away from the electrodes are insulating, and the gradient of the

potential normal to the interfaces are zero. In order to implement these condition, we deﬁne

an additional set of mesh points along these interfaces, and we then apply the insulating

boundary condition by setting,

ψ(i, k + 1)

ψ(i + 1, k)

ψ(i, k) =

i /∈ [0, Rtip/h], k = M (top boundary)
i = N, k = {0, 1, 2, ..., M} (side boundary)

for

for

(3.67)

These four boundary conditions (3.66)-(3.67) together with the discretized Poisson’s equa-

tion (3.63) - (3.64) result in a system of linear equation, which can be represented in the

matrix form (3.48); where vector (cid:126)b contains the values of hole carrier density and boundary
conditions, and Matrix A reﬂects the morphological eﬀect of the spatially-dependent dielec-
tric constant. Unlike the 1-D DD-SCLC model, the matrix A derived from tip-plane geometry

is not a tridiagonal matrix, and the inversion formula suggested by Usmani can no longer be

used. In our study, we used the biconjugate gradient stabilized method (BiCGSTAB) [68] to

ﬁnd solutions to this non-symmetric linear system iteratively.

58

3.4.3 Hole Density Solver

By pursuing the Schafetter-Gummel approach [54], hole current densities at the mesh mid-

points in cylindrical coordinates are expressed as

Jp,ρ(i ± 1

2 , k) = µp,ρ(i ± 1

2 , k) · qVtNv

∆ρ

(cid:40)
B [ψ(i ± 1, k) − ψ(i, k)] · pf (i, k)

(cid:41)
∓B [ψ(i, k) − ψ(i ± 1, k)] · pf (i ± 1, k)

≡ µp,ρ(i ± 1

qµp,z(i, k ± 1
∆z

Jp,z(i, k ± 1

2) =

(cid:40)
2 , k) · Jp,ρ(i ± 1
2 , k)
B [ψ(i, k ± 1) − ψ(i, k)] · pf (i, k)

2)VtNv

(cid:41)
∓B [ψ(i, k) − ψ(i, k ± 1)] · pf (i ± 1, k)

(3.68a)

≡ µp,z(i, k ± 1

2) · Jp,z(i, k ± 1
2)

(3.68b)

where B is the Bernoulli function. The angular component of hole current density (Jp,ϕ) is

zero, according to the cylindrical symmetry of the system.

Similar to the numerical method used for Poisson’s equation, the ﬂux of hole current

density is calculated on the chosen Gaussian surfaces. We then obtain

59

At z-axis (ρ = 0 or i = 0)

(cid:34)(cid:73)

(cid:126)Jp · d (cid:126)A

SG

=

i,k

(cid:90)(cid:90)
(cid:90)(cid:90)

SG(i+ 1
2 )

+
SG(k+ 1
2 )

(cid:40)

µp,ρ ρdϕdz · Jp,ρ(i + 1

2 , k) −

(cid:90)(cid:90)

µp,z ρdρdϕ · Jp,z(i, k − 1
2)

SG(k− 1
2 )

µp,z ρdρdϕ · Jp,z(i, k + 1
2)

2 , k)

4¯µp,ρ(i + 1

(cid:18)∆ρ
(cid:18)∆ρ

∆z

∆z

−

+

2, k) · Jp,ρ(i + 1
(cid:19)2
(cid:19)2

¯µp,z(i, k − 1

¯µp,z(i, k + 1

2) · Jp,z(i, k − 1
2)
2) · Jp,z(i, k + 1
2)

(cid:41)

=

≡ 0

· 1
8

∆ϕ∆z

(3.69)

Otherwise (ρ = i · ∆ρ where i = {1, 2, 3, .., N})

(cid:126)Jp · d (cid:126)A

= −

µp,ρ ρdϕdz · Jp,ρ(i − 1

2 , k) +

(cid:34)(cid:73)

SG

(cid:90)(cid:90)
(cid:90)(cid:90)

SG(i+ 1
2 )

SG(k+ 1
2 )

µp,ρ ρdϕdz · Jp,ρ(i + 1

2 , k)

µp,z ρdρdϕ · Jp,z(i, k + 1
2)

(cid:35)

(cid:35)

i,k

=

SG(i− 1
2 )
−
SG(k− 1
2 )
1 − 1
−
2i

(cid:40)

(cid:90)(cid:90)
(cid:90)(cid:90)
(cid:18)
(cid:18)
(cid:18)∆ρ
(cid:18)∆ρ

∆z

1 +

∆z

+

−

+

µp,z ρdρdϕ · Jp,z(i, k − 1

2) +

· ¯µp,ρ(i − 1

(cid:19)
(cid:19)
(cid:19)2 · ¯µp,z(i, k − 1
(cid:19)2 · ¯µp,z(i, k + 1

· ¯µp,ρ(i + 1

1
2i

2 , k)

2 , k)

2 , k) · Jp,ρ(i − 1
2 , k) · Jp,ρ(i + 1
2) · Jp,z(i, k − 1
2)
2) · Jp,z(i, k + 1
2)

(cid:41)

· i∆z∆φ

≡ 0

(3.70)

60

where the local averages of the normal component of hole mobility over each section of

Gaussian surface around the mesh point (i, k) are

(cid:105)

·(cid:104)

¯µp,ρ(i ± 1

2, k) =

¯µp,z(i, k ± 1

2) =

1
2



µp,ρ(i ± 1

2 , k − 1

2) + µp,ρ(i ± 1

;

2, k + 1
2)
for i = {0, 1, 2, .., N},

(3.71a)

µp, z(i, k ± 1
2);
2 , k ± 1
2) ·
µp,z(i − 1

(cid:18)

(cid:19)

for i = 0,
1 − 1
4i

(cid:18)

(cid:19)

(3.71b)

+ µp,z(i + 1

2

1 +

2 , k ± 1
2) ·

1
4i
for i = {1, 2, 3, .., N},

;

and k = {0, 1, 2, ..., M}. By choosing the uniform grid where ∆ρ = ∆z = h, the discretized

equations (3.69)-(3.70) are simpliﬁed as

At z-axis (ρ = 0 or i = 0)

4¯µρ(i + 1

2 , k) ·(cid:104)B(∆ρψ(i + 1, k)) · pf (i, k) − B(−∆ρψ(i + 1, k)) · pf (i + 1, k)
(cid:105)
2) ·(cid:104)B(∆zψ(i, k)) · pf (i, k − 1) − B(−∆zψ(i, k)) · pf (i, k)
(cid:105)
(cid:105)
2) ·(cid:104)B(∆zψ(i, k + 1)) · pf (i, k) − B(−∆zψ(i, k + 1)) · pf (i, k + 1)

− ¯µz(i, k − 1

+ ¯µz(i, k + 1

= 0 (3.72)

61

Otherwise (ρ = i · ∆ρ where i = 1, 2, 3, ...)

(cid:18)
(cid:19)

(cid:18)

+

1 +

−

1
2i

(cid:19)

1 − 1
2i

· ¯µp,ρ(i + 1

· ¯µp,ρ(i − 1

2 , k) ·(cid:104)B(∆ρψ(i, k)) · pf (i − 1, k) − B(−∆ρψ(i, k)) · pf (i, k)
2, k) ·(cid:104)B(∆ρψ(i + 1, k)) · pf (i, k) − B(−∆ρψ(i + 1, k)) · pf (i + 1, k)
2) ·(cid:104)B(∆zψ(i, k)) · pf (i, k − 1) − B(−∆zψ(i, k)) · pf (i, k)
2) ·(cid:104)B(∆zψ(i, k + 1)) · pf (i, k) − B(−∆zψ(i, k + 1)) · pf (i, k + 1)

(cid:105)
(cid:105)
(cid:105)
(cid:105)

−¯µp,z(i, k − 1

+¯µp,z(i, k + 1

= 0

(3.73)

where the spatial diﬀerence of electric potential is calculated as follow

∆ρψ(i, k) = ψ(i, k) − ψ(i − 1, k)
∆zψ(i, k) = ψ(i, k) − ψ(i, k − 1)

(3.74a)

(3.74b)

Boundary conditions are required to solve the discretized continuity equation (3.72)-(3.73)

numerically. At the electrodes, the boundary conditions are described as

1

1

pf (i, k) =

i ∈ [0, Rtip/h], k = M
i ∈ [0, R/h], k = 0

for

for

(3.75)

Following the same technique used when obtaining the insulating boundary conditions in

the potential solver, we have

pf (i, k + 1)

pf (i + 1, k)

pf (i, k) =

i /∈ [0, Rtip/h], k = M (top boundary)
i = N, k = {0, 1, 2, ..., M} (side boundary)

for

for

(3.76)

62

Once again, these four boundary conditions together with the discretized continuity equa-

tions (3.72)-(3.73) constitute a system of linear equations. In the corresponding matrix no-
tation (3.53), the vector (cid:126)b contains the boundary conditions, and the matrix A encompasses

the morphological eﬀect through the spatially-dependent hole mobility. Unlike the 1-D DD-
SCLC model, the matrix A is asymmetric, and the numerical solutions of free hole carrier

density (cid:126)pf is achieved by the iterative method called BiCGSTAB [68].

3.4.4 Veriﬁcation

Similar to the work done with the 1-D DD-SCLC model, we have developed a Fortran 90

program to ﬁnd the unknown electric potentials and free hole densities of 2-D DD-SCLC

models in cylindrical coordinate systems with cylindrical symmetry by following the algo-

rithm shown in ﬁgure 3.3. Due to the asymmetric properties emerged from the tip-plane

geometry, the direct method of inversion used in the 1-D DD-SCLC model is no longer appli-

cable. Therefore, we have developed a Fortran 90 module to execute a routine of BiCGSTAB,

which is an iterative method to solve the non-symmetric linear system. The programing code

is written speciﬁcally to handle sparse matrices. To minimize the error associated with nu-

merical calculation, all parameters and variables are represented by double precision values

in the program.

First of all, we would like to verify the BiCGSTAB module by comparing the results with

those acquired from the direct method. By retroﬁtting the module into the 1-D DD-SCLC

model, we simulate charge transport of trap-free hole-only devices in planar geometry. The

comparison of the current-voltage characteristics calculated from the two distinct methods,

inversion of tridiagonal matrices and BiCGSTAB, is illustrated in ﬁgure 3.12. The error

between the two methods is negligible.

63

Figure 3.12: BiCGSTAB VS Inversion – The theoretical current-voltage characteristic of
the trap-free (ft = 0) hole-only device solved iteratively by using BiCGSTAB module is
compared with that solved directly by inversion of a tridiagonal matrix. In both calculations,
we use a 1-D DD-SCLC model with L = 80nm, Nv = 1.25 × 1021cm−3,r = 3, µp =
10−4cm2V −1s−1.

64

Applied Voltage (V)101102103104105Current Density (A/m2)InversionBiCGSTAB0.11100.00000.00100.00200.00300.0040%DifferenceFigure 3.13: The comparison of the theoretical current-voltage characteristics of the trap-
free (ft = 0) hole-only device calculated from 1-D and 2-D DD-SCLC models.
In both
models, the transport layer (orange) is homogeneous with the constant parameters: r = 3,
µp = 10−4cm2V −1s−1, and Nv = 1.25 × 1021cm−3.

65

Next, we would like to verify the 2-D DD-SCLC model by reproducing the result of

homogeneous planar device. Consider the trap-free hole-only device of thickness 80nm and

radius 300nm. To capture the planar case, the tip radius is extended to the size of the device

radius. In ﬁgure 3.13, the current-voltage characteristic simulated from the 2-D DD-SCLC

model is compared with that simulated from the 1-D DD-SCLC model. Although the 2-D

DD-SCLC model has the possibility of the additional lateral transport compared to the 1-D

DD-SCLC model where only vertical transport occurs, their current-voltage characteristics

exhibit an inconsequential error less than 0.0035%. Therefore, only vertical conductance

matters in charge transport of these planar devices, and the 2-D DD-SCLC model is veriﬁed

for this case.

Replacing one of the planar electrodes with a tip-shaped electrode, one can expect a

diﬀerent amount of charge injection due to their distinct contact areas. Moreover, unlike in

the conventional planar device, the charge transport in a tip-plane device is subject to the

inﬂuence of the lateral current spreading induced by the tip geometry. As the 2-D DD-SCLC

model can only simulate the transport behavior in ﬁnite systems surrounded by insulating

boundaries, we need to ensure that the computational domain is suﬃciently large to make the

boundary eﬀect negligible. To demonstrate this eﬀect, we perform the numerical calculations

to model the hole-only device, where holes are injected from a tip electrode of 1nm radius

into a device of 80 nm thickness. The other electrode is planar of the same size as the device

radius, which will be varied in the study. It is important to note that the approximate 1nm

tip radius is the smallest tip allowed in the discretized domain of 2nm spacing, where it

can be assigned to a single mesh point. The broadest lateral spreading corresponds to the

smallest tip size.

Unlike the case of planar geometry, the hole-only device in tip-plane geometry no longer

66

maintain the uniform vertical current. Therefore, it is only informative to present the overall

hole transport in term of a total current, which is the surface integral of current density over

the device boundary. Due to the cylindrical symmetry of the 2-D DD-SCLC model, the total

current in the lateral direction is zero. Therefore, the total hole current injected at the tip

is simply the surface integral of the vertical current over the cross-sectional surface of the

device, which is

R(cid:90)

0

Ip = 2π

Jp,z(ρ, z) ρdρ.

(3.77)

This is the physical quantity that can be measured from the experiments with the tip-plane

electrode conﬁguration, e.g., the c-AFM measurements. It is also useful to emphasize that the

total current is independent of the vertical position of the device cross-section, i.e., the total

current is constant through the device. This manner is governed by the continuity equation
at the steady state. In the discretized domain (ρ ∈ [0, R] and ρ = i· ∆ρ; i = {0, 1, 2, .., N}),

total current is calculated from

Ip =

π
4

∆ρ2 Jp,z(0, k + 1

2) + 2π∆ρ2

N(cid:88)

i=1

iJp,z(i, k + 1
2)

(3.78)

where Jp,z(i, k + 1

2) is the discretized current density associated to Schafetter-Gummel dis-

cretization, as written in eq.(3.68b).

Figure 3.14 plots the hole current as a function of the radius to thickness ratio (R/L) of

the device, where the simulation is carried out using a 2-D DD-SCLC model with tip-plane

geometry. The transport layer is assumed to be homogeneous with isotropic hole mobility.

For a ﬁxed device thickness, the current approaches the saturated regime when device radius

is suﬃciently large so that the ﬁnite-size eﬀect can be safely ignored. A good balance point

67

Figure 3.14: Finite-size eﬀect–Convergence of the simulated hole current with respect to
device dimensional ratio of radius to thickness, when the applied voltage Va = 0.1, 1, 10V .
By using a 2-D DD-SCLC model, the trap-free (ft = 0) hole-only device is modeled in
tip-plane geometry (Rtip ∼ 1nm), assuming homogeneous transport layer (r = 3, Nv =
1.25 × 1021cm−3) with isotropic hole mobility (µp,ρ = µp,z = 10−4cm2V −1s−1). Using the
dimensional ratio R/L = 3.75 provides a well-converged solutions while the ﬁnite-size eﬀect
is negligible. All the following calculations of tip-plane devices with an isotropic mobility
will be performed with this dimensional ratio.

68

R/L550600650700750800Current  (pA)Va = 10 V02468100.0750.0760.077Va = 0.1V6.577.5Va =   1 VFigure 3.15: Finite-size eﬀect–Convergence of the simulated hole current with respect to
device dimensional ratio of radius to thickness, when the applied voltage Va = 0.1, 1, 10V .
By using a 2-D DD-SCLC model, the trap-free (ft = 0) holse-only device is modeled in
tip-plane geometry (Rtip ∼ 1nm), assuming homogeneous transport layer (r = 3, Nv =
1.25 × 1021cm−3) with anisotropic hole mobility (µp,ρ/µp,z = 25, µp,z = 10−4cm2V −1s−1).
Using the dimensional ratio R/L = 6.25 provides a well-converged solutions while the ﬁnite-
size eﬀect is negligible. All the following calculations of tip-plane devices with anisotropic
mobility will be performed with this dimensional ratio.

69

R/L200030004000500060007000Current  (pA)Va = 10 V0510150.450.50.550.6Va = 0.1V3040506070Va =   1 Vwith the feasible numerical domain and negligible ﬁnite-size eﬀect is achieved at the device

dimensional ratio of 3.75, which incurs less than a 1% diﬀerence compared to the converged

numerical solution of a system with dimensional ratio of 10.

In addition to the isotropic system, we have also considered the case of hole-only device

with anisotropic mobility. When the lateral mobility is larger than the vertical one, we

expect an increasing amount of current spreads out beneath the tip, which consequently

enhances the ﬁnite-size eﬀect. The evidence is shown in ﬁgure 3.15. Current predicted from

the simulation with anisotropic mobility (µρ/µz = 25) converges at higher dimensional ratio

compared to that of a homogeneous device with isotropic mobility. Empirically, to safely

omit the ﬁnite-size eﬀect at reasonable computational overhead, the simulation of tip-plane

device with large anisotropic mobility (µρ/µz = 25) shall be conducted at the dimensional

ratio of 6.25, which results in approximately a maximum 2% diﬀerence compared to the

converged numerical solutions from devices with ratio R/L = 12.50.

The phenomena of current spreading was emphasized in the recent work by Ginger and

his coworkers [32], in which they modiﬁed the well-known Mott-Gurney equation to include

this tip-induced eﬀect, and established a semi-empirical formula to describe current-voltage

characteristics in the SCLC regime measured from c-AFM. Because it is important to un-

derstand the c-AFM experimental results, we extend their theoretical model to include the

morphological eﬀect, i.e., exponentially distributed trap and anisotropic mobility. We also

alter the insulating boundary condition to increase the stability and improve the convergence

of the calculation, especially in the region near the tip electrode. Here, we would like to ver-

ify that the insulating boundary of device is well described by our choice of the insulating

boundary conditions used in the model, as expressed in eq.(3.59). To explicitly insulate the

device, the top boundary that is not the injecting electrode is replaced by a vacuum layer

70

Figure 3.16: Current as a function of vacuum layer thickness simulated using the given device
structure, compared to the current calculated from a 2-D DD-SCLC model with insulating
boundary conditions described by eq.(3.59), when Va = 0.5, 1, 10V . At high vacuum layer
thickness, the consistency of the numerical current between the two device geometries is
observed with approximately a maximum of 1% diﬀerence. For the transport layer (orange),
we use L = 80nm, r = 3, µp,ρ = µp,z = 10−4cm2V −1s−1, an5d Nv = 1.25 × 1021cm−3.

71

of vacuum permittivity (r = 1) and zero hole mobility ( (cid:126)µp = 0). The injecting electrode is

transformed from a circular contact of 1nm radius to a cylindrical rod of the same radius.

The conducting rod is then lengthened through the thickness of the vacuum layer, and con-

nected to the topmost planar electrode. We have conﬁrmed that the device radius is large

enough that we can safely keep the insulating boundary condition on the side of the devices.

The schematic representation of the new device structure is shown in ﬁgure 3.16.

Similarly, we carried out the numerical calculation to model the homogeneous hole-only

device of 80nm thickness and 300nm radius with varying thickness of the insulating vacuum

layer. Figure 3.16 plots the current versus the vacuum layer thickness when the applied

voltage Va = 0.5, 1, 10V . The currents simulated from the 2-D DD-SCLC model with the

insulating boundary condition, as expressed in eq.(3.59), are plotted (dashed red line) for

the comparison. It is shown that, with the increasing thickness of the vacuum layer, the

simulated currents approach the saturated regime, demonstrating a good consistency (< 1%

diﬀerence at the 100nm vacuum layer thickness) with the current calculated using only the

insulating boundary conditions. Note that the maximum diﬀerence between the two device

geometries occurs at the high applied voltage, which is Va = 10V in our calculations. The

explanation to this observation is related to the ﬁeld distortion caused by the planar electrode

at the topmost interface that becomes more severe with increasing applied voltage.

3.5 The Fully 3-D DD-SCLC Model in a Cartesian Co-

ordinate System

In the previous sections, we have introduced a 1-D DD-SCLC model to describe the vertical

charge transport in a homogeneous planar device, and a 2-D DD-SCLC model that allows us

72

to observe the additional lateral transport in a tip-plane device with cylindrical symmetry.

It is already well-known that electric conductance of an disordered organic semiconductor is

considerably aﬀected by its nanoscale morphology. Superior charge transport has been found

in semi-crystalline materials due to the presence of a nanoﬁbrillar network. To investigate

the eﬀect of such complex morphology, a fully three-dimensional model is obviously a good

choice. The model also permits the implementation of tip geometry and spatially-varying

traps.

Figure 3.17: Illustration of a three-dimensional hole-only device of tip-plane geometry. The
numerical calculation is performed using a three-dimensional coordinate system, allowing
full customization of trap distribution and nanoscale morphology in the transport layer.

Figure 3.17 illustrates the three-dimensional device structure modeled as a Lx × Ly × Lz

cuboid of tip-plane geometry. For simplicity, we approximate the tip electrode as a square
contact of Ltip side length. By enlarging the tip to cover the area Lx × Ly, the device

structure turns into a planar geometry. Similar to the 1-D and 2-D DD-SCLC models,

charge transport in the hole injection system can be described by three governing equations,

including Poisson’s equation, the continuity equation at steady state and the drift-diﬀusion

73

equation. We have

(cid:73)

S=∂V

((cid:126)r)(cid:126)E((cid:126)r) · d (cid:126)A =

qNv
Vt

(cid:73)

S=∂V

(cid:0)pf ((cid:126)r) + S((cid:126)r)ft · pf ((cid:126)r)rt(cid:1) dV

(cid:90)

V

(cid:126)Jp((cid:126)r) · d (cid:126)A = 0

(cid:0)pf ((cid:126)r)Eξ((cid:126)r) − ∂pf /∂ξ(cid:1)

(3.79)

(3.80)

(3.81)

Jp,ξ((cid:126)r) = qµp,ξ((cid:126)r)NvVt

Jp,ξ((cid:126)r) · ˆξ and ξ = {x, y, z}. The internal inhomogeneity, i.e., nanoscale

where (cid:126)Jp((cid:126)r) =(cid:80)
dielectric constant ((cid:126)r) = r((cid:126)r)o, hole mobility (cid:126)µp((cid:126)r) =(cid:80)

morphology and traps, is included in the model via the spatially varying parameters of
µp,ξ((cid:126)r) · ˆξ, and trap distribution
S((cid:126)r). Poisson’s equation and the continuity equation are presented in the integral form,

ξ

ξ

which is convenient to treat the inhomogeneity under the discretization.

Similar to the 2-D DD-SCLC model, the device is grounded on the bottom electrode

(z = 0), and driven at the top electrode (z = Lz) with the applied potential Va. The top

boundary is insulating. Therefore, the boundary condition at the top and bottom boundaries

are deﬁned as

ψ(x, y, Lz) = Va/Vt

and pf (x, y, Lz) = 1,

(cid:12)(cid:12)(cid:12)(cid:12)(x,y,Lz)

∂ψ
∂z

(cid:12)(cid:12)(cid:12)(cid:12)(x,y,Lz)

= 0,

= 0 and

∂pf
∂z

where {x, y} ∈ Dtip
where {x, y} /∈ Dtip

(3.82a)

(3.82b)

ψ(x, y, 0) = 0

and pf (x, y, 0) = 1,

where x ∈ [0, Lx], y ∈ [0, Ly]

(3.82c)

where Dtip is the two-dimensional domain of the tip contact area, e.g., in the case of square
tip shown in ﬁgure 3.17, we have Dtip ∈ [Ltip, Ltip]. It should be noted that the 3-D DD-

74

SCLC model does not limit to this particular tip geometry, but can be generalized for any

shape of tip contact that can be discretized by the numerical approach explained in the

following section. Furthermore, for the more complex condition that the conducting tip

is penetrated into the material, which commonly occurs in the c-AFM measurements, the
tip is thus modeled by deﬁning Dtip as the three-dimensional domain of the tip geometry

penetrated into the material domain. This is absolutely executable in the numerical tool

we developed for the 3-D DD-SCLC model. However, to keep the derivation simple, we will

assume that the penetration depth is negligible and the tip is approximated by the surface

contact area at the top boundary.

In contrast to the 2-D DD-SCLC model, we apply periodic boundary conditions on

the four sides of the device to mitigate ﬁnite-size eﬀects. However, when considering tip-

plane geometry, the device dimensions need to be large enough to allow complete current

spreading. Based on the analysis of the 2-D DD-SCLC model (shown in ﬁgure 3.14-3.15), the

proper dimensions for the device with isotropic mobility and extreme anisotropic mobility

(µp,x = µp,y = 25 µp,z) are Lx = Ly = 3.75 Lz and Lx = Ly = 6.25 Lz respectively.

3.5.1 Discretization

Choosing a uniform mesh separation equally in all directions (∆x = ∆y = ∆z = h), the

three-dimensional device domain is replaced by a number of cube-shaped cells of volume
h3 and a set of mesh points (xi, yj, zk) located on the edge of the cell; i = {0, 1, 2, .., Nx},
j = {0, 1, 2, ..., Ny}, and k = {0, 1, 2, .., Nz}, where xi = i · h, yj = j · h, zk = k · h,

and h = Lx/Nx = Ly/Ny = Lz/Nz, as shown in ﬁgure 3.18(a). A single mesh point is

surrounded by eight adjacent cells and six nearest neighbors listed in ﬁgure 3.18(b), and

each cell is assumed to carry invariable values of dielectric constant, hole mobility, and trap

75

parameters. The internal inhomogeneity of the material is then represented by the variation

of cell properties from one cell to the other.

(a)

(b)

(c)

Figure 3.18: Illustration of the spatial discretization in a three-dimensional Cartesian coor-
dinate system; (a) the device is divided into cube cells (h × h × h) with mesh points located
at each corner; (b) each mesh point is surrounded by six nearest neighbors, and is adjacent
to eight cells; (c) each cell containss constant values of dielectric constant, hole mobility, and
trap parameters. The inhomogeneity is generated by varying the values of these parameters
from one cell to another.

Figure 3.19: Gaussian surface SG as a cube of size h× h× h surrounding mesh point (i, j, k)

Considering Poisson’s equation and the continuity equation in the integral form, we
introduce a cube-shaped Gaussian surface (SG) bounding the ﬁnite volume h × h × h, as

shown in ﬁgure 3.19. Similar to the 2-D DD-SCLC model, the numerical surface and volume

integration of the physical variables around each mesh point will be performed with respect to

its Gaussian surface based on three assumptions: (i) hole current density is constant between

76

mesh points, (ii) electric ﬁeld is also constant between mesh points, (iii) the tangential

component of electric ﬁeld is continuous at the interface of neighboring cells. It is worth

emphasizing that by applying the assumption (ii) and (iii), the normal component of electric

ﬁeld at the Gaussian surface is constant and continuous due to the fact that the Gaussian

surface is chosen to be perpendicular to the interface of neighboring cells.

3.5.2 Potential Solver

By applying the numerical approximation of surface integral (3.27) and volume integral

(3.28) on the chosen Gaussian surface, Poisson’s equation (3.79) in the three-dimensional

Cartesian coordinate system can be discretized as

(cid:34)(cid:73)

(cid:35)

(cid:126)E · d (cid:126)A

= −

SG

i,j,k

((cid:126)r) dydz · Ex(i − 1

2, j, k) +

((cid:126)r) dydz · Ex(i + 1

2 , j, k)

(cid:90)(cid:90)
(cid:90)(cid:90)
(cid:90)(cid:90)

(cid:90)(cid:90)
(cid:90)(cid:90)
(cid:90)(cid:90)

SG(i+ 1
2 )

SG(j+ 1
2 )

((cid:126)r) dxdz · Ey(i, j + 1

2, k)

((cid:126)r) dxdy · Ez(i, j, k + 1
2)

2 , k) +

2) +

((cid:126)r) dxdy · Ez(i, j, k − 1

((cid:126)r) dxdz · Ey(i, j − 1

SG(i− 1
2 )
−
SG(j− 1
2 )
−
SG(k− 1
2 )
= − ¯x(i − 1
2 , j, k) h2 Ex(i − 1
2 , k) h2 Ey(i, j − 1
− ¯y(i, j − 1
2) h2 Ez(i, j, k − 1
− ¯z(i, j, k − 1

(cid:2)pf (i, j, k) + ftS(i, j, k) · pf (i, j, k)rt(cid:3) .

2) + ¯z(i, j, k + 1

≡ qNvh3
Vt

2 , j, k) + ¯x(i + 1

2 , k) + ¯y(i, j + 1

SG(k+ 1
2 )
2 .j, k) h2 Ex(i + 1
2 , k) h2 Ey(i, j + 1
2 , k)
2) h2 Ez(i, j, k + 1
2)

2 , j, k)

(3.83)

77

Follow the numerical approximations for the ﬁrst order derivative, and the discretized electric

ﬁelds are then given by

Ex(i ± 1
Ey(i, j ± 1
Ez(i, j, k ± 1

2 , j, k) = [ψ(i, j, k) ∓ ψ(i ± 1, j, k)] /h
2 , k) = [ψ(i, j, k) ∓ ψ(i, j ± 1, k)] /h
2) = [ψ(i, j, k) ∓ ψ(i, j, k ± 1)] /h

where the local averages of dielectric constant over each Gaussian surface (SG) are

(3.84a)

(3.84b)

(3.84c)

(3.85a)

(3.85b)

(3.85c)

(cid:105)
(cid:105)
(cid:105)

¯x(i ± 1

2 , j, k) =

¯y(i, j ± 1

2 , k) =

¯z(i, j, k ± 1

2) =

(cid:104)
(cid:104)
(cid:104)

(i ± 1
+(i ± 1
(i − 1
+(i − 1
(i − 1
+(i − 1

2 , j − 1

2 , k − 1

2) + (i ± 1

2 , j + 1

2 , k − 1
2)

2 , j − 1

2 , k + 1

2) + (i ± 1

2 , j + 1

2 , k + 1
2)

/4

2 , j ± 1

2 , k − 1

2) + (i + 1

2 , j ± 1

2 , k − 1
2)

2 , j ± 1

2 , k + 1

2) + (i + 1

2 , j ± 1

2 , k + 1
2)

/4

2 , j − 1

2 , k ± 1

2) + (i + 1

2 , j − 1

2 , j + 1

2 , k ± 1

2) + (i + 1

2 , j + 1

2 , k ± 1
2)
2 , k ± 1
2)

/4

The discretized Poisson’s equation (3.83) is then rewritten in terms of electric potential as

− ¯x(i − 1

2 , j, k) · ψ(i − 1, j, k) − ¯x(i + 1

2 , j, k) · ψ(i + 1, j, k)

− ¯y(i, j − 1

2 , k) · ψ(i, j − 1, k) − ¯y(i, j + 1

− ¯z(i, j, k − 1

2) · ψ(i, j, k − 1) − ¯z(i, j, k + 1

+ 6 ¯(i, j, k) · ψ(i, j, k) =

qNvh2

Vt

2 , k) · ψ(i, j + 1, k)
2) · ψ(i, j, k + 1)

(cid:2)pf (i, j, k) + ftS(i, j, k) · pf (i, j, k)rt(cid:3) (3.86)

78

where ¯(i, j, k) is the locally averaged dielectric constant at mesh point (i, j, k), which is

equivalent to the arithmetic average of the values from eight neighboring cells surrounding

the considering mesh point. Meanwhile, the ¯(i, j, k) can also be expressed in term of the

locally averaged values described in equation (3.85) as follows,

(cid:104)

¯(i, j, k) =

¯x(i − 1

2, j, k) + ¯x(i + 1

2 , j, k) + ¯y(i, j − 1

2 , k)

(cid:105)

+¯y(i, j + 1

2, k) + ¯z(i, j, k − 1

2) + ¯z(i, j, k + 1
2)

/6.

(3.87)

In order to achieve the unique solution, the boundary conditions are required. Adopting

the technique we used when deriving the discretized boundary conditions in the 2-D DD-

SCLC, we have the complete set of electrode- and insulating- conditions for both top (k =

Nz) and bottom (k = 0) boundary as follow.

Va/Vt

ψ(i, j, Nz + 1)

ψ(i, j, Nz) =

ψ(i, j, 0) = 0

for

for

for

(i, j) ∈ Dtip
(i, j) /∈ Dtip
i = {0, 1, 2, .., Nx}, j = {0, 1, 2, ..., Ny}

(3.88a)

(3.88b)

Four side boundaries of the system satisfy the periodic boundary condition, such that

ψ(0, j, k) = ψ(Nx, j, k)

ψ(i, 0, k) = ψ(i, Ny, k)

for

for

j = {0, 1, 2, .., Ny}, k = {0, 1, 2, ..., Nz}
i = {0, 1, 2, .., Nx}, k = {0, 1, 2, ..., Nz}

(3.89a)

(3.89b)

By substituting these boundary conditions into the discretized Poisson’s equation (3.86),

we have the system of linear equations in matrix and vector notation (3.48); where vector

79

(cid:126)b reﬂects the disordered properties of the material in terms of trapped hole density, while
matrix A carries the morphological eﬀect of the spatially varying dielectric constant. To ﬁnd

the numerical solution of potential ( (cid:126)ψ), we carried out the iterative calculation using the

BiCGSTAB module.

3.5.3 Hole Density Solver

By pursuing the Schafetter-Gummel approach [54], hole current densities at the mesh mid-

points in the three-dimensional Cartesian coordinate system are expressed as

Jp,x(i ± 1

2 , j, k) = µp,x(i ± 1

2 , j, k) · qVtNv

≡ µp,x(i ± 1
2 , k) = µp,y(i, j ± 1

Jp,y(i, j ± 1

≡ µp,y(i, j ± 1

qµp,z(i, j, k ± 1

h

Jp,z(i, k ± 1

2) =

(cid:40)
B [ψ(i ± 1, j, k) − ψ(i, j, k)] · pf (i, j, k)
(cid:41)
h
∓B [ψ(i, j, k) − ψ(i ± 1, j, k)] · pf (i ± 1, j, k)

(3.90a)

(cid:40)
2 , j, k) · Jp,x(i ± 1
2 , j, k)
B [ψ(i, j ± 1, k) − ψ(i, j, k)] · pf (i, j, k)
2 , k) · qVtNv
(cid:41)
h
∓B [ψ(i, j, k) − ψ(i, j ± 1, k)] · pf (i, j ± 1, k)
(cid:40)
B [ψ(i, j, k ± 1) − ψ(i, j, k)] · pf (i, j, k)

(cid:41)
∓B [ψ(i, j, k) − ψ(i, j, k ± 1)] · pf (i, j, k ± 1)

2 , k) · Jp,y(i, j ± 1

2)VtNv

2 , k)

(3.90b)

≡ µp,z(i, j, k ± 1

2) · Jp,z(i, j, k ± 1
2)

(3.90c)

80

where B is the Bernoulli function. Similar to the numerical method for Poisson’s equation,

the ﬂux of hole current density is calculated on the chosen Gaussian surfaces. We obtain

µp,x dydz · Jp,x(i − 1

2 , j, k) +

µp,x dydz · Jp,x(i + 1

2 , j, k)

(cid:34)(cid:73)

(cid:35)

(cid:126)Jp · d (cid:126)A

= −

SG

i,j,k

(cid:90)(cid:90)
(cid:90)(cid:90)
(cid:90)(cid:90)

SG(i− 1
2 )
−
SG(j− 1
2 )
−
SG(k− 1
2 )
− ¯µp,x(i − 1

(cid:34)

=

(cid:90)(cid:90)
(cid:90)(cid:90)
(cid:90)(cid:90)

SG(i+ 1
2 )

SG(j+ 1
2 )

µp,y dxdz · Jp,y(i, j − 1

2 , k) +

µp,z dxdy · Jp,z(i, j, k − 1

2) +

SG(k+ 1
2 )
2 , j, k) · Jp,x(i − 1, j, k) + ¯µp,x(i + 1

− ¯µp,y(i, j − 1
− ¯µp,z(i, j, k − 1

2 , k) · Jp,y(i, j − 1, k) + ¯µp,y(i, j + 1
2) · Jp,z(i, j, k − 1) + ¯µp,z(i, j, k + 1

µp,y dxdz · Jp,y(i, j + 1

2 , k)

µp,z dxdy · Jp,z(i, j, k + 1
2)

2 , j, k) · Jp,x(i + 1, j, k)
(cid:35)
2 , k) · Jp,y(i, j + 1, k)
2) · Jp,z(i, j, k + 1)

h2

≡ 0

(3.91)

where the locally-averaged hole mobilities over the Gaussian surface are written as,

¯µp,x(i ± 1

2 , j, k) =

¯µp,y(i, j ± 1

2, k) =

¯µp,z(i, j, k ± 1

2) =

(cid:104)
(cid:104)
(cid:104)

µp,x(i ± 1
+ µp,x(i ± 1
µp,y(i − 1
+ µp,y(i − 1
µp,z(i − 1
+ µp,z(i − 1

2 , j − 1

2 , k − 1

2) + µp,x(i ± 1

2 , j + 1

2 , k − 1
2)

2 , j − 1

2 , k + 1

2) + µp,x(i ± 1

2 , j + 1

2 , k + 1
2)

/4

(3.92a)

2 , j ± 1

2 , k − 1

2) + µp,y(i + 1

2 , j ± 1

2 , k − 1
2)

2 , j ± 1

2 , k + 1

2) + µp,y(i + 1

2 , j ± 1

2 , k + 1
2)

/4

2, j − 1

2 , k ± 1

2) + µp,z(i + 1

2 , j − 1

2, j + 1

2 , k ± 1

2) + µp,z(i + 1

2 , j + 1

2 , k ± 1
2)
2 , k ± 1
2)

/4.

(3.92b)

(3.92c)

(cid:105)
(cid:105)
(cid:105)

81

This is similar to the locally-averaged dielectric constants, as in eq.(3.85), but not as simpli-

ﬁed due to the fact that hole mobility is vector.

Substituting the discretized current density (3.90), we have the discretized continuity

equation for the 3-D DD-SCLC model.

−¯µp,x(i − 1

(cid:105)
2, j, k) ·(cid:104)B(∆xψ(i, j, k)) · pf (i − 1, j, k) − B(−∆xψ(i, j, k)) · pf (i, j, k)
2, j, k) ·(cid:104)B(∆xψ(i + 1, j, k)) · pf (i, j, k) − B(−∆xψ(i + 1, j, k)) · pf (i + 1, j, k)
(cid:105)
2, k) ·(cid:104)B(∆yψ(i, j, k)) · pf (i, j − 1, k) − B(−∆yψ(i, j, k)) · pf (i, j, k)
(cid:105)
(cid:105)
2 , k) ·(cid:104)B(∆yψ(i, j + 1, k)) · pf (i, j, k) − B(−∆yψ(i, j + 1, k)) · pf (i, j + 1, k)
2) ·(cid:104)B(∆zψ(i, j, k)) · pf (i, j, k − 1) − B(−∆zψ(i, j, k)) · pf (i, j, k)
(cid:105)
(cid:105)
2) ·(cid:104)B(∆zψ(i, j, k + 1)) · pf (i, j, k) − B(−∆zψ(i, j, k + 1)) · pf (i, j, k + 1)

−¯µp,z(i, j, k − 1

+¯µp,z(i, j, k + 1

+¯µp,x(i + 1

−¯µp,y(i, j − 1

+¯µp,y(i, j + 1

where the spatial diﬀerence of electric potential is calculated as follows

∆xψ(i, j, k) = ψ(i, j, k) − ψ(i − 1, j, k)
∆yψ(i, j, k) = ψ(i, j, k) − ψ(i, j − 1, k)
∆zψ(i, j, k) = ψ(i, j, k) − ψ(i, j, k − 1)

= 0 (3.93)

(3.94a)

(3.94b)

(3.94c)

Boundary conditions are needed to ﬁnd the unique solutions of free hole densities. For the

hole density solver of the 3-D DD-SCLC model, the complete set of electrode- and insulating-

82

conditions for the top (k = Nz) and bottom (k = 0) boundaries are discretized as follow

1

pf (i, j, Nz + 1)

pf (i, j, Nz) =

pf (i, j, 0) = 1

for

for

for

(i, j) ∈ Dtip
(i, j) /∈ Dtip
i = {0, 1, 2, .., Nx}, j = {0, 1, 2, ..., Ny}

(3.95a)

(3.95b)

where Dtip is the two-dimensional domain of the tip contact. Note that the top boundary

not corresponding to the tip contact is treated with the insulating boundary conditions,

following the numerical approach previously used in the 2-D DD-SCLC model. Lastly, Four

boundaries on the sides, i.e., S(i = 0), S(i = Nx), S(j = 0) and S(j = Ny), satisfy the

periodic boundary conditions, such that

pf (0, j, k) = pf (Nx, j, k)

pf (i, 0, k) = pf (i, Ny, k)

(3.96a)

(3.96b)

Following similar mathematics to that used for the potential solver, the combination of the

boundary conditions (3.95)-(3.96) and the continuity equation (3.93) results in an asymmetric

system of linear equations, which can be rewritten in the simple matrix form (3.53). While
vector (cid:126)b contains the boundary conditions, matrix A captures the morphological eﬀect in

terms of the spatially-varying hole mobility. The numerical solutions for (cid:126)pf are achieved by

the iterative method BiCGSTAB [68].

83

3.5.4 Veriﬁcation

We have developed a Fortran 90 program to simulate current-voltage characteristics of the

fully three-dimensional model in a Cartesian coordinate system by following the algorithm

shown in ﬁgure 3.3. Due to the similar mathematical structure, the BiCGSTAB module

previously used in the 2-D DD-SCLC model is implemented to handle both the potential

and hole density solvers. As a standard protocol, all parameters and variables are represented

by double precision values in the program, in order to minimize the error arisen from the

iterative method used in the numerical calculations. To verify the 3-D DD-SCLC model,

we carried out the numerical calculations of two limited cases that have been studied before

using the 1-D and 2-D DD-SCLC models.

First, we would like to verify the 3-D DD-SCLC model by reproducing the results of
the homogeneous planar devices. Consider the trap-free hole-only device of 600nm× 600nm

cross-sectional area and 80nm thickness. To resemble the planar geometry of the 1-D DD-

SCLC model, the tip contact area is extended to the size of device’s cross-section. Figure

3.20 shows the comparison of current-voltage characteristics of a trap-free hole-only device

simulated from 1-D and 3-D DD-SCLC models. It is clear that the numerical results from

the two diﬀerent models are identical.

Secondly, we would like to include the 2-D DD-SCLC model in this part of veriﬁcation

and consider the device simulations in the tip-plane geometry. As schematically shown in
ﬁgure 3.21, we portrait two possible tip designs to be constructed on the 2nm× 2nm square

grid of the 3-D DD-SCLC model, including (a) a circular contact of 6nm radius and (b)

a square contact of 10nm side length, which are chosen to be roughly consistent with a

circular tip used in the 2-D DD-SCLC model, as depicted by the red circle of 6nm radius.

84

Figure 3.20: The comparison of the theoretical current-voltage characteristics of trap-free
(ft = 0) hole-only devices calculated from 1-D and 3-D DD-SCLC models. In both models,
the transport layer (orange) is homogeneous with the constant parameters: r = 3, µp =
10−4cm2V −1s−1, and Nv = 1.25 × 1021cm−3.

85

Figure 3.21: The schematic representation of tip contacts in a 3-D DD-SCLC model, con-
structed on the 2nm × 2nm square grid as (a) an approximately circular geometry of 6nm
radius and (b) a square geometry of 10nm side-length. In a 2-D DD-SCLC model, the tip
contact is assigned as the circular area of 6nm radius (red dashed circle). Due to the diﬀer-
ence in the discretization, the tip models in the 3-D model are slightly smaller that in the
2-D model, such that A3−D(a)/A2−D = 0.85 and A3−D(b)/A2−D = 0.88, respectively.

We then carried out the numerical calculations to simulate charge transport in the trap-free

hole-only devices. In ﬁgure 3.22-3.23, the current-voltage characteristics simulated from the

3-D DD-SCLC model are compared with those simulated from the 2-D DD-SCLC model.

In the models of the inset, the calculation domain is considered as a circular cylinder of

300nm radius and 80nm thickness for the 2-D DD-SCLC model, and as a rectangular prism
of 300nm × 300nm cross-sectional area and 80nm thickness for the 3-D DD-SCLC model.

As illustrated in ﬁgure 3.22-3.23, we observe a good agreement of the current-voltage

characteristics of hole-only devices calculated from the two models.

In the case of the

round tip model, as shown in ﬁgure 3.22, the maximum of approximately 16% diﬀerence is

expected, and can be explained by the slightly diﬀerent contact area caused by the limitation

of discretization in the cylindrical and Cartesian coordinate systems. However, we can narrow

down these diﬀerences in the numerical solutions between the two models by using the square

tip contact, as shown in ﬁgure 3.23.

86

Figure 3.22: The comparison of the theoretical current-voltage characteristics of the trap-free
(ft = 0) hole-only devices calculated from 2-D and 3-D DD-SCLC models, in which the tip
electrode is approximated by a circular contact of 6nm radius, as illustrated in ﬁgure 3.21(a).
In both models, the transport layer (orange) is homogeneous with the constant parameters:
r = 3, µp = 10−4cm2V −1s−1, and Nv = 1.25 × 1021cm−3.

87

Figure 3.23: The comparison of the theoretical current-voltage characteristics of the trap-
free (ft = 0) hole-only devices calculated from 2-D and 3-D DD-SCLC models, in which
the tip electrode is approximated by a circular contact of 6nm radius and a square contact
of 10nm side length, respectively. (See ﬁgure 3.21(b)) In both models, the transport layer
(orange) is homogeneous with the constant parameters: r = 3, µp = 10−4cm2V −1s−1, and
Nv = 1.25 × 1021cm−3.

88

3.6 Conclusions

In this chapter the DD-SCLC model for hole injection systems was introduced. At the

center of the model are three governing equations including Poisson’s equation, the continuity

equation at steady state, and the drift-diﬀusion equation for holes. Adopting the numerical

technique presented by Kosteret al. [52], all governing equations are solved simultaneously

using a self-consistent scheme. We have developed the numerical tools written in Fortran

90 to carry out the numerical calculations for one-, two- and three- dimensional DD-SCLC

models, which provide the insight to charge transport of hole-only devices with diﬀerent

geometries.

• The one-dimensional (1-D) DD-SCLC model is the simplest one among the three mod-

els, yet very essential. The numerical calculations are performed in a one-dimensional

Cartesian coordinate system, and capable to simulate both drift and diﬀusion currents

in homogeneous planar devices, in which traps are assumed to be uniformly distributed,

with dielectric constant and hole mobility being constant. According to this descrip-

tion, the one-dimensional model improves the simpliﬁed SCLC model by just including
the eﬀect of charge diﬀusion. We have shown that a Log J − Log V plot simulated

from our model exhibits slope 2 in the drift-dominated regime (Va > 10V ). This is the

signature of Mott-Gurney equation, validating our DD-SCLC model.

• The two-dimensional (2-D) DD-SCLC model is constructed in the cylindrical coordi-

nate system with the cylindrical symmetry. The model is very eﬃcient to study SCLC

in tip-plane geometry with mobility anisotropy.

• The fully three-dimensional (3-D) DD-SCLC model is a big leap forward in the simu-

89

lation of SCLC. The numerical calculations are carried out using a three-dimensional

Cartesian coordinate system, enabling the treatment of inhomogeneous systems, in-

cluding the tip-plane electrode conﬁguration, spatially-varying trap distributions, and

nanoscale morphologies. However, the simulation is highly demanding in both com-

puting time and memory.

90

Chapter 4

The Drift-Diﬀusion SCLC Transport

Models - New Results

In this chapter, we present theoretical studies of hole-only devices, including charge diﬀu-

sion, traps, and nanoscale morphologies; in both planar and tip-plane (c-AFM) geometries.

Examples of the applications of the drift-diﬀusion (DD) SCLC model will be carried out us-

ing the numerical tools we successfully developed for the one-, two- and three- dimensional

DD-SCLC models described in the previous chapter. The parameters used in the device sim-

ulations are summarized in Table 4.1. Note that these values are the typical material and

device parameters, which should be useful to capture the general physics of organic semicon-

ductor devices. Our theoretical models do not limit to this particular system. In the next

chapter, these methods and insights will be applied to an analysis of c-AFM experimental

data provided by Dr. Pengpeng Zhang’s group at Michigan State University.

4.1 Diﬀusion Eﬀects

In chapter 2, we presented the simpliﬁed SCLC model for hole injection and transport

problems in the conventional planar geometry. The exact analytical solution including the

current-voltage characteristic, as known as the Mott-Gurney equation, and other physical

91

Table 4.1: Summary of parameters used in the device simulations

Parameter

Symbol Numerical value

Dielectric constant
The eﬀective density of states in the
HOMOa
The eﬀective hole mobility in the verti-
cal direction
Temperature
Device thickness

2-D DD-SCLC model :
Circular tip contact radius b
Circular device radiusc

εr
Nv

µo

T

L, Lz

Rtip
R

3.0
1.25 × 1021 cm−3
1.0 × 10−4 cm2/V s

300 K
80 nm

2 nm
300 nm for the isotropic transport
750 nm for the anisotropic transport

3-D DD-SCLC model :
Square tip contact side lengthd
Ltip
Square device cross-section side lengthe Lx, Ly

5 nm
600 nm for the isotropic transport
1500 nm for the anisotropic transport

aThe number of Nv is determined from the regime where the current density is saturated over the value

of Nv, as shown in ﬁgure 3.7.

bFor the 2-D model, the computing domain is uniformly discretized with the 2nm mesh interval. There-
fore, we choose the smallest tip size allowed by the discretization precess to maximize the concentration of
the electric ﬁeld around the tip, and pronouncing the eﬀect of lateral spreading.

cThe circular device radius is determined from the regime where the current density is saturated over the
value of R/L, as shown in ﬁgure 3.14-3.15 for the case of isotropic- and anisotropic- hole mobility respectively.
dFor the 3-D model, the computing domain is uniformly discretized with the 5nm mesh interval. Same

as b, we choose the smallest tip size allowed by the discretization process.

eSame as c.

92

variables (such as electric potential, hole carrier density, etc.) are summarized in table

2.1. In ﬁgure 4.2, we plot the analytical proﬁles of the electric potential and hole carrier

density in order to illustrate the characteristics of this hole-only SCLC transport model. For

simplicity, the variables are plotted in their natural scales, such that the potential is scaled

by the externally applied voltage Va, the hole density is scaled by its spatial average ¯p (from

table 2.1), and the internal device coordinate is scaled by the device thickness L.

Figure 4.1: Based on the simpliﬁed SCLC model of a trap-free hole-only device in the planar
geometry, the analytical proﬁles of electric potential (solid black line) and hole carrier density
(dashed red line) are demonstrated using non-dimensional scales.

Although charge diﬀusion is omitted from the simpliﬁed SCLC model, the exact solution

of hole charge carrier density exhibits a non-uniform distribution through the device thick-

ness. When a gradient of charge carrier density exists, a process of diﬀusion occurs such that

the charge carriers migrate from the regions of high density toward regions of low density.

Assuming that the hole diﬀusion coeﬃcient obeys the Einstein relation, the corresponding

93

hole diﬀusion current in the absence of traps is given by

(cid:16) x

(cid:17)−3/2

Jp,di =

3
8

µp

VtVa
L3

L

(4.1)

Both drift and diﬀusion current densities are dependent on the externally applied voltage
a and Jp,di ∝ Va, respectively. Since the drift current density
Va, such that Jp,dr ∝ V 2

increases quadratically with the applied voltage while the diﬀusion current density increases

only linearly, the drift current density will dominate in the high voltage regime where the

diﬀusion process is negligible.

Furthermore, per eq.(4.1), the diﬀusion current density is also non-uniform through the

device thickness. It is inﬁnite at the injecting electrode (x = 0), then dramatically declines,

and gradually saturates toward the extracting electrode (x = L). This behavior explains one

considerable discrepancy of the simpliﬁed SCLC model on the boundary condition at the

injecting electrode. The model assumes the injecting contact to be ohmic, whence the local
electric ﬁeld E = 0 and consequently the drift current vanishes. Therefore, the model clearly

requires an additional form of charge transport to account for hole injection at the electrode.

That is where the diﬀusion current comes into play, which is basically inﬁnite at the ohmic

contact of an inﬁnite hole reservoir. In order to evaluate how signiﬁcant the diﬀusion eﬀect

is compared to the drift eﬀect in the model, ﬁgure 4.2 shows the plot of the current density

ratio Jp,di/Jp,dr versus the dimensionless coordinate along the direction of device thickness

x/L. It is convenient to introduce the characteristic coordinate xc as the location where the

current density ratio is unity. The characteristic coordinate is given by

(cid:18)3Va

(cid:19)−2/3

xc = L

Vt

94

(4.2)

Figure 4.2: Plot of the analytical ratio of diﬀusion current density to drift current density of
a planar ﬂow, trap-free, hole-only device, Jp,di/Jp,dr, versus the dimensionless coordinate,
x/L. While the drift current density is the immediate analytic solution of the simpliﬁed
SCLC model as known as Mott-Gurney equation (2.13), the diﬀusion current density is
obtained from the further derivation of the analytic solutions of the model. The extremely
large value of diﬀusion current near the injecting electrode (x = 0) suggests the diﬀusion
assists hole injection at the injecting contact.

In the region near xc, diﬀusion and drift current densities are of the same order of magnitude.

The diﬀusion current density then become dominant toward the injecting electrode (x < xc).

However, the dominance of charge diﬀusion diminishes with the increasing of the applied

voltage. When the applied voltage is suﬃciently high, the hole injection system is then well

described by drift transport mechanism alone.

It is important to note that, by including the approximate diﬀusion current density

expressed in equation (4.1) to the total current density, the continuity equation at steady

state is indeed violated and total current is no longer conserved. This is expected since the

diﬀusive transport is not initially considered in the simpliﬁed SCLC model. The drift current

may dominate over the diﬀusion current if the externally applied voltage is suﬃciently high,

95

00.20.40.60.81x / L0.00010.0010.010.1110100100010000Jp,di  / Jp,drVa = 0.1VVa =    1VVa =  10Vbut it is not thoroughly understood whether the diﬀusion current can be ignored or not in

the intermediate range of applied voltage. Therefore, even where such analytical solutions

are available, it is still useful to contemplate the full SCLC model that incorporates both

drift and diﬀusion transport mechanisms, the so-called standard drift-diﬀusion model. Then

the eﬀect of diﬀusion process can be studied by solving the drift-diﬀusion equation coupled

with Poisson’s equation under the conservation law of current density described by continuity

equation at steady state.

Now that we have developed a computer program that solves the three governing equa-

tions numerically in a self-consistent manner, it is instructive to ﬁrst compare the theoretical

proﬁles of electric potential and hole carrier density simulated from numerical calculation to

that derived from the simpliﬁed SCLC model, previously demonstrated in ﬁgure 4.1. Using

the parameters given in table 4.1, the simulations of hole-only devices with two distinct

geometries; planar (plane-parallel) and tip-plane (c-AFM) geometries, will be investigated

in detail in the following subsections.

In order to focus solely on the inﬂuence of charge

diﬀusion, the charge transport layer is assumed to be homogeneous and trap free, similar to

that of the simpliﬁed SCLC model.

4.1.1 Planar Geometry

To begin with, let us examine the conventional case of planar hole-only device. The the-

oretical proﬁles of electric potential and hole carrier density as a function of the distance

from the injecting electrode with the varying applied voltage (Va = 0.1, 1, 10V ) are shown

in ﬁgure 4.3, which have all been generated by a 1-D DD-SCLC model of a hole-only de-

vice described in section 3.3. The numerical simulation includes both drift and diﬀusion

transport mechanisms, resulting in the non-monotonic distributions of electric potential and

96

hole carrier density. Together, the theoretical proﬁle can be categorized into three distinct

regimes separated by two local transitions of charge-transport behaviors. The ﬁrst transition

is determined by the position of a potential maximum where electric ﬁeld and drift current

vanish. The second transition is located at the position of a hole density minimum where

diﬀusion current becomes zero. Since the calculations are performed on the device at steady

state, the total current which is the sum of drift and diﬀusion currents is conserved. By

keeping this fact in mind, we will be able to understand the charge-transport mechanism in

each regime qualitatively.

In the transport regime (I), space charges are built up near the injecting electrode (x = 0),

resulting in a potential barrier. By overcoming the reverse drift current, the charge diﬀusion

therefore is the dominant transport mechanism driving holes away from the injecting contact

into the device. After passing the potential barrier and entering into regime (II), holes then

travel forward driven by both drift and diﬀusion mechanisms. Finally, in regime (III), the

diﬀusion current is reversed, counteracting the drift current, and leading to the reduction

of the drift current toward the collecting electrode at x = 80nm. This is due to the ohmic

nature of the extracting contact.

When quantitatively comparing the theoretical proﬁle simulated from our 1-D DD-SCLC

model (ﬁgure 4.3) to the analytical proﬁle derived from the simpliﬁed SCLC model (ﬁgure

4.1), we observe subtle diﬀerences in the applied voltage dependence of the distribution. In

the simpliﬁed SCLC model, the electric potential and hole carrier density have monotonic

distribution along the device thickness, while the 1-D DD-SCLC model yields a more com-

plicated result consisting of three distinct regimes. A similar observation has been described

in the past by A.A. Grinberg and S. Luryi, particularly for the case of double junction (p-

i-p) diode [69]. The idea can be extended to the case of hole-only devices by introducing

97

Figure 4.3: The theoretical proﬁle of electric potential and hole carrier density as a function
of the distance from the injecting electrode for planar geometry, generated using the 1-D DD-
SCLC model described in section 3.3, using parameters given in table 4.1. The proﬁles exhibit
three distinct transport regimes, plotted for three diﬀerent applied voltages Va = 0.1, 1, 10V ,
from top to bottom.

98

a virtual anode and a virtual cathode, located at the position of zero drift current (zero

electric ﬁeld) and zero diﬀusion current, respectively. The two virtual electrodes essentially

mark the separation of the transport regimes (I), (II), and (III) on the simulated proﬁle.

With the increasing of the applied voltage, the virtual electrodes move toward the actual

(ﬁxed) electrodes. This phenomenon is governed by the charge diﬀusion, which is a missing

transport mechanism in the simpliﬁed SCLC model.

4.1.2 Tip-Plane Geometry

Apart from the conventional planar geometry, the eﬀect of charge diﬀusion is further studied

by modifying the injecting electrode into a circular conducting tip, similar to that of c-AFM

measurement. Assuming that the transport layer is homogeneous with an isotropic hole

mobility, the device will have a cylindrical symmetry around the injecting electrode. This

makes it easy to simulate the tip-plane electrode conﬁguration using the 2-D DD-SCLC

model, as presented in section 3.4. Charge transport comprises of two components; the

lateral transport along the radial direction (ρ) and the vertical transport along the thickness

direction (z). The inﬂuence of charge diﬀusion on transport behavior of tip-plane electrode

conﬁguration can be observed along the thickness direction in the same manner as that of

plane-parallel electrode conﬁguration. By choosing the vertical path at the center of the

circular tip/electrode (ρ = 0), the lateral transport vanishes and only the vertical transport

remains. As the result, ﬁgure 4.4 shows the theoretical proﬁle of electric potential and hole

carrier density as a function of z with the varying applied voltage Va = 0.1, 1, 10V , where

z is the coordinate along the vertical axis from tip center pointing to planar electrode. By

applying the concept of virtual electrodes, similar to that of planar geometry, the theoretical

proﬁle of the tip-plane geometry exhibits three transport regimes:

99

Figure 4.4: The theoretical proﬁle of electric potential and hole carrier density at the center
of a circular tip (ρ = 0) as a function of the distance from the injecting electrode for tip-plane
geometry, generated by the 2-D DD-SCLC model described in section 3.4, using parameters
given in table 4.1. The proﬁles exhibits three distinct transport regimes, plotted for three
diﬀerent applied voltages Va = 0.1, 1, 10V .

100

• regime (I) where the charge diﬀusion is the dominant form of charge transport, located

between the actual anode (z = 0) and the virtual anode,

• regime(II) of the drift-diﬀusion assisted transport located between the virtual anode

and the virtual cathode,

• regime (III) where drift becomes the dominant form of charge transport, located be-

tween the virtual cathode and the actual cathode (z = L).

Similarly to what has been observed in the study of planar geometry, the virtual electrodes

of tip-plane geometry also move closer to the actual electrodes with the increasing of the

externally applied voltage.

To illustrate the eﬀect of the applied voltage on the movement of the virtual electrodes,

we plot the displacement of virtual anode and virtual cathode as a function of applied

voltage for both plane-parallel and tip-plane electrode conﬁgurations, as illustrated in ﬁgure

4.5. Due to the absence of intrinsic charge carriers and low carrier mobility in an organic

semiconductor, without the external potential bias applied to the electrodes (the so-called

short-circuit condition), there is no actual charge transport in the device. Although no

current can be measured, theoretically we can demonstrate the non-uniform distribution

of electric potential and charge carrier density. For the planar hole-only device, this is

analogous to the case of a parallel plate capacitor ﬁlled with a dielectric, and the simulation

yields a symmetric arrangement that both virtual electrodes are located at the same position

half-way between the two actual electrodes. By applying the positive bias at the actual

anode, holes are injected into the device, traveling through the transport layer via drift

and diﬀusion mechanisms, and then extracted at the actual cathode. That is when the

position of the virtual electrode starts shifting toward the actual one. In the situation of

101

Figure 4.5: Position of virtual anode (dashed line) and virtual cathode (solid line) as a
function of applied voltage for planar (black) and tip-plane (red) geometries, as seen in
ﬁgure 4.3 and ﬁgure 4.4 respectively. The inset shows the convergence of the virtual tip
position.

102

00.511.52Applied Voltage (V)020406080Position (nm)[Planar Geo.] V. Cathode [Tip-Plane Geo.] V. Cathode[Planar Geo.] V. Anode[Tip-Plane Geo.] V. Anode0246810Applied Voltage (V)020406080Position (nm)planar geometry, by increasing the applied voltage, the virtual electrodes move away from

each other symmetrically. However, in the tip-plane electrode conﬁguration, the symmetric

behavior no longer holds. At the short-circuit condition, the virtual electrodes of the tip-

plane geometry are located toward the injecting electrode, instead of the middle of the device

thickness. This is due to an asymmetric pair of tip and plane electrodes and a lateral current

spreading beneath the tip-like structure.

Although the concept of virtual electrodes is chosen for the explanation of charge trans-

port involving the additional charge diﬀusion, it is not the new to the simpliﬁed SCLC

model. But, per the chosen boundary conditions at the electrodes, the simpliﬁed SCLC

model assumes that the virtual electrodes are ﬁxed as the actual ones. The inset of ﬁgure

4.5 shows the convergence of the virtual electrodes’ position with the externally applied volt-

age. When the virtual electrodes are co-located with the actual electrode, the inﬂuence of

charge diﬀusion is negligible and the current-voltage characteristic of planar and tip-plane

geometries are then well described by Mott-Gurney equation (2.13) and the semi-empirical

formula (2.31) introduced by Ginger’s group, respectively.

4.2 Trap Eﬀects

Up to now, we have assumed that hole-only devices are free of traps. This situation typically

does not exist since an organic semiconductor is commonly known for its energetic disorder.

The topic of traps in organic semiconductors has been intensively studied by experiments on

unipolar injection devices in planar geometry when either electrons or holes are injected into

the devices. The corresponding current-voltage characteristics in the SCLC regime are then

analyzed using the theory of trap-limited SCLC, as written in eq.(2.27). This is a variation

103

of the Mott-Gurney equation that has been developed to include trap eﬀects. Therefore,

the derivation leading to eq.(2.27) is based on two assumptions: (i) the electric ﬁeld at the

injecting electrode vanishes and (ii) the charge diﬀusion is negligible. These assumptions are

only valid when the electric ﬁeld across the device is suﬃciently high. More importantly,

traps often lead to the condition that diﬀusion current can no longer be neglected, which

makes the drift-only transport model problematic. [45, 70, 71]

For a more accurate description of trap-limited SCLC, we developed the theoretical model

described in chapter 3 that incorporates both drift and diﬀusion transport mechanisms. In

the model, the energetically distributed traps are described by an exponential distribution

within the band gap, characterized by trap density through the parameter ft = Nt/Nv and

trap temperature through the parameter rt = T /Tt, as expressed in eq.(3.4). With varying

trap parameters (ft, rt), the corresponding current-voltage characteristics will be analyzed

in the SCLC regime. We carry out the numerical calculations in two diﬀerent cases; planar

and tip-electrode geometries, using parameters given in table 4.1. In this section, traps are

assumed to be uniformly distributed throughout the device.

4.2.1 Planar Geometry

The eﬀect of traps in the hole injection system can be studied by comparing the theoret-

ical current-voltage characteristics simulated from the planar device with the presence of

the exponentially distributed traps to that simulated from the trap-free planar device. As

previously demonstrated in ﬁgure 3.8, the current-voltage characteristic of a trap-free hole
injection device obeys the power law (J ∝ V m), showing the transition of charge transport

behavior from the diﬀusion-inﬂuenced regime at low applied voltage (m < 2) to the drift-

dominant regime at high applied voltage (m = 2). This characteristic slope of 2 is the key

104

aspect of trap-free SCLC in the drift-only transport model, as predicted by the well-known

Mott-Gurney equation (2.13). Therefore, it shall be instructive to present and analyze the

current-voltage characteristics of hole-only devices in the double logarithmic scale where the

slope is the exponent m. We will mostly show the current-voltage characteristics in log-log

representation from this point forward.

Figure 4.6 shows the theoretical current-voltage characteristic of a hole-only device with

a low density of deep traps (ft = 0.005, rt = 0.30) generated by the 1-D DD-SCLC model

(described in section 3.3). Unlike the trap-free current-voltage characteristic, one can dis-

tinguish three distinct regimes corresponding to Ohm’s law (m = 1), trap limited SCLC

(m > 2), and trap-ﬁlled SCLC (m = 2). In the absence of injection barriers, the number

of holes injected into the device increases with an increase in the externally applied voltage.

In the low voltage regime, the amount of injected holes is low, and become negligible com-

pared to that of the intrinsic thermally activated holes. The electrical conduction in this

regime is then due to the intrinsic charge carriers, and hence can be described by Ohm’s

law (m = 1). When increasing the applied voltage, the device accommodates an increasing

number of hole carriers. The rapid growth of current with the small increment of applied

voltage is caused by the ﬂow of injected holes through the device. In an organic semicon-

ductor device, this bulk transport is typically controlled by space charges. The intermediate

voltage regime where the exponent m > 2 is uniquely related to the exponential distribution

of traps within the band gap. When trap states are completely ﬁlled, the exponent then

reduces to 2 as predicted by the Mott-Gurney equation for trap-free SCLC. Note that the

theoretical current-voltage characteristic is plotted with the data of applied voltage up to

100V in order to observe all three key regimes. The working range of applied voltage for an

organic semiconductor device is normally below 10V , which is large enough to observe the

105

Figure 4.6: The trap inﬂuenced current-voltage characteristic of a hole-only device in planar
geometry obeys the power law J ∝ V m, exhibiting three regimes: Ohm’s law m = 1 (red),
trap-limited SCLC m > 2 (green), and trap-ﬁlled SCLC m = 2 (blue).

106

0.010.1110100Applied Voltage (V)10-210-1100101102103104105106107Current Density (A/m2)ft = 5e-3, rt = 0.30m = 1.00m = 2.89m = 2.0trap ﬁlling regime.

In order to determine the dependence of current-voltage characteristics on traps, we

performed a large number of device simulations over the wide range of trap parameters:
ft = 0−0.05 and rt = 0.20−0.40. Examples of the theoretical current-voltage characteristics

of hole-only devices in planar geometry are shown in ﬁgure 4.7 when (a) ft = 0.005 with

varying rt, and (b) rt = 0.30 with varying ft. These data are plotted in the range of the

intermediate applied voltage regime involving the trap ﬁlling process in SCL transport, and

can be described by the semi-empirical expression

J(ft, rt) = µo · kJ (ft, rt) · V m(ft,rt)

(4.3)

with the ﬁtting parameters kJ and m depending on trap parameters ft and rt. By using

this expression to ﬁt all simulated current-voltage characteristics, we obtain the database

of ﬁtting parameters kJ and m which lead to the understanding of how traps aﬀect the

current-voltage relations. This analysis is very challenging due to its nature of nonlinear

regression and the large amount of data involved. Therefore, we have developed a program

written in Python to carry out this extensive analysis. As a result, the ﬁtting parameters

are summarized in table A.1 of appendix A.

Although we cannot present the semi-empirical formula of the ﬁtting parameters, we

observe a monotonic dependence of the ﬁtting exponent m on traps as illustrated by the

contour plot in ﬁgure 4.8b. The exponent m varies from 1.89 to 5.00 corresponding to trap
parameter variations over the ranges of ft = 0 − 0.05 and rt = 0.20 − 0.40. The contour
plot demonstrates that there are many possible pairs of trap parameters {ft, rt} along the
contour line which can give the same slope of Log J − Log V curve in the SCLC regime.

107

(a)

(b)

Figure 4.7: The theoretical current-voltage characteristics of planar hole-only devices in the
SCLC regime obey the power law (I ∝ V m) when (a) ft = 0.005 with varying rt yielding
the exponent in the range m = 2.23 − 4.71, and (b) rt = 0.30 with varying ft yielding the
exponent in the range m = 1.89 − 3.72

108

110Applied Voltage (V)10-1100101102103104105Current Density (A/m2)rt = 0.40rt = 0.35rt = 0.30rt = 0.25rt = 0.20110Applied Voltage (V)101102103104105Current Density (A/m2)ft = 0ft = 5e-04ft = 1e-03ft = 5e-03ft = 0.010ft = 0.025ft = 0.050(a)

(b)

Figure 4.8: Dependence of the exponent m on the two trap parameters; trap density char-
acterized by parameter ft and trap temperature characterized by parameter rt for hole-only
devices in planar geometry based on (a) the simpliﬁed SCLC model as described by eq.(2.27),
(b) 1-D DD-SCLC model as seen in ﬁgure 4.7 and described by eq.(4.3).

109

However, in the strong trap limit where the trap temperature is high and the trap density

is large, the slope depends only on the trap temperature, while in the intermediate regime,

the slope is inﬂuenced by both the trap temperature and the trap density. This dependence

is quite diﬀerent from what is predicted by the simpliﬁed SCLC model described by eq(4.7),

in which the exponent m only depends on the parameter rt related to trap temperature, as

shown in ﬁgure 4.8a.

4.2.2 Tip-Plane Geometry

Taking into account the tip-plane electrode conﬁguration, we simulate the charge injection

and transport in hole-only device using the 2-D DD-SCLC model described in section 3.4. As

the ﬁrst step, the hole mobility is assumed to be isotropic. Figure 4.9 shows the theoretical

current-voltage characteristic of a hole-only device in tip-plane geometry with a low density

of deep traps (ft = 0.005, rt = 0.30), exhibiting three distinct regimes corresponding to

• Ohm’s law in the low voltage regime,

• trap-limited SCLC in the intermediate voltage regime, and

• trap-ﬁlled SCLC in the high voltage regime.

This is similar to what has been observed in the theoretical current-voltage characteristic of

the hole-only device in planar geometry (previously demonstrated in ﬁgure 4.6). The regime

of interest is the trap-limited SCLC, which typically ranges from 1V to 10V . More examples

of the theoretical current-voltage characteristics in this intermediate regime are shown in

ﬁgure 4.10 when (a) ft = 0.005 with varying rt and (b) rt = 0.30 with varying ft. This
again suggests that the slope of the Log I − Log V relations, i.e., the exponent m, depends

on trap parameters.

110

Figure 4.9: The trap inﬂuenced current-voltage characteristic of a hole-only device in the
tip-plane geometry obeys the power law J ∝ V m, exhibiting three regimes: Ohm’s law
m = 1 (red), trap limited SCLC m > 2 (green), and trap ﬁlled SCLC m = 2 (blue). Hole
mobility is assumed to be isotropic.

111

0.010.1110100Applied Voltage (V)10-710-610-510-410-310-210-1100101102103104105106Current (pA)ft = 5e-3, rt = 0.30m = 1.00m = 3.74m = 2.00While the system has cylindrical symmetry, the hole mobility used in the calculation can

be non-uniform. In this section, we will treat the mobility anisotropy in a simple manner,

where we assume that the mobility is uniform in the lateral (ρ-) and vertical (z-) direction,

and yet with diﬀerent values of µρ and µz, respectively. The situation often occurs in

the experiments of semi-crystalline organic semiconductors (See chapter 5). The degree

of mobility anisotropy is deﬁned by the ratio µρ/µz. As a preliminary study, we perform

the numerical calculation with varying anisotropic mobility ratio up to 25, which has been

observed in hole-transporting P3HT thin ﬁlm composed of semi-crystaline whiskers [67].

Figure 4.11 shows the theoretical current-voltage characteristics in the low deep trap regime

(ft = 0.005, rt = 0.30) with varying of the mobility anisotropic ratios (µρ/µz). Interestingly,

the exponent m in the power law does not show a dependence on the mobility anisotropy.

Following the analysis previously done on the trap study in the planar geometry, we car-

ried out extensive calculations of the dependence of the theoretical current-voltage character-

istics on the trap parameters and the mobility anisotropy. We found that the current-voltage

relations can be described by the semi-empirical expression,

I(µρ/µz, ft, rt) = µo · kI (µρ/µz, ft, rt) · V m(ft,rt)

(4.4)

where kI is the ﬁtting parameter and m is the power-law exponent. While the parameter

kI depends on both trap parameters and mobility anisotropy, the exponent m only depends

on trap parameters. We have applied our Python program to ﬁt the simulated current-

voltage characteristics in the SCLC regime and have extracted the ﬁtting values of kI and

m. The ﬁtting results are summarized in table A.2 of appendix A. By analyzing the slope
of Log I − Log V curves, we observe the dependence of the exponent m on trap parameters

112

(a)

(b)

Figure 4.10: The theoretical current-voltage characteristics of tip-plane devices in the SCLC
regime obey the power law (I ∝ V m) when (a) ft = 0.005 with varying rt yielding the
exponent in the range m = 2.89 − 5.52, and (b) rt = 0.30 with varying ft yielding the
exponent in the range m = 2.00 − 4.19. Hole mobility is assumed to be isotropic.

113

0.515Applied Voltage (V)10-410-310-210-1100101102Current (pA)rt = 0.40rt = 0.35rt = 0.30rt = 0.25rt = 0.200.515Applied Voltage (V)10-210-1100101102Current (pA)ft = 0ft = 5e-04ft = 1e-03ft = 5e-03ft = 0.010ft = 0.025ft = 0.050Figure 4.11: The theoretical current-voltage characteristics of tip-plane devices in the SCLC
regime for a range of mobility anisotropic ratios, µρ/µz, with ﬁxed trap parameters ft =
0.005, rt = 0.30.

114

Figure 4.12: Dependence of the exponent m on the exponentially distributed trap DOSs
characterized by parameters ft and rt for hole-only devices in tip-plane geometry in the SCLC
regime, resulting from extensive calculations for theoretical current-voltage characteristics
as demonstrated in ﬁgure 4.10 and described by the semi-empirical expression (4.4).

115

as illustrated by the contour plot in ﬁgure 4.12. The exponent m varies from 2.00 to 5.73
corresponding to trap parameter variations over the ranges of ft = 0 − 0.05 and rt =
0.20 − 0.40. The dependence of the exponent m on trap parameters resulting from device

simulations in the tip-plane geometry is qualitatively similar to that of the planar geometry

(shown in ﬁgure 4.8b).

4.3 Morphological Eﬀects: Semi-Crystaline P3HT Thin

Films

It has been shown that the overall performance of organic semiconductor devices have

been improved signiﬁcantly with the structural enhancement of their nanoscale morpholo-

gies [72, 73]. Semiconducting P3HT is among the most extensively studied polymers in

the ﬁeld, being the best-known representative that exhibits highly crystalline self-assembled

nanostructures [1, 65, 66, 67, 74, 75, 76, 77]. Therefore, in this section, we carried out a the-

oretical study of the morphological eﬀect in semi-crystalline P3HT ﬁlms by analyzing the

theoretical current-voltage characteristics simulated from tip-plane device geometry in anal-

ogy to c-AFM measurements.

Figure 4.13 shows a schematic illustration of crystallized P3HT nanoﬁber. In general,
the self-organized P3HT exhibits a lamellar ordering via π − π interactions that lengthen

the nanoﬁber. Based on the alignment of P3HT lamellae with respect to the substrate,

the nanoﬁber can be categorized into three orientation schemes: (i) ‘face-on’ orientation

with the lamellae being parallel to the substrate, (ii) ‘edge-on’ orientation with the lamellae

being perpendicular to the substrate, and (iii) ‘end-on’ orientation with both lamellae and

polymer-backbone being perpendicular to the substrate, which are schematically shown in

116

(a)

(b)

(c)

Figure 4.13: Schematic illustration of crystallized P3HT nanoﬁbers in three possible ori-
entation schemes, including (a) ‘face-on’, (b) ‘edge-on’ and (c) ‘end-on’ orientations, corre-
sponding to the parallel and perpendicular alignments of self-organized P3HT lamellae to
the substrate, respectively.

117

ﬁgure 4.13 (a),(b) and (c), respectively. While polymer crystallinity depends signiﬁcantly on

regioregularity [65,78,79], molecular weight [80,81,82,83] and fabrication process [84,85,86],
P3HT typically forms micrometer long nanoﬁbers along the π − π stacking direction with

a width of 20nm along the polythiophene backbone direction and a thickness of 5nm along

the alkyl side chain direction [1, 31, 36, 37, 38, 87]. Furthermore, it has been reported that

producing the ‘end-on’ ﬁbers is very challenging due to the diﬃculty of stacking up the long

backbone direction [88, 89]. Therefore, in this study, we will focus on the ‘edge-on’ and

‘face-on’ ﬁbers, which is commonly found in the semi-crystalline P3HT thin ﬁlms.

Experimentally, it has been shown that the crystallized P3HT nanoﬁber exhibits mobility
anisotropy [65,67,90,91], such that the mobility along the π−π stacking direction (ﬁber axis)

is very high as compared to the mobility along the alkyl side chain direction, which is com-

parable to that of amorphous P3HT. Clearly, the high mobility channel through nanoﬁbers

is preferable for charge transport. While it is straightforward to make an assumption that

the presence of nanoﬁbers could lead to the enhancement of the overall device mobility,

experiments [34, 67] have shown that the ‘edge-on’ ﬁbers do not signiﬁcantly increase the

mobility measured in the plane-parallel electrode conﬁguration (planar geometry) and may

slightly lower the mobility as compared to amorphous ﬁlms. This observation can be ex-

plained by the characteristics of charge transport in the conventional planar geometry. In

the experiments of this particular geometry, the current-voltage response is mainly subject

to the vertical charge transport, and consequently the exceptionally high lateral mobility

of ‘edge-on’ ﬁbers may be left undetected. To obtain the comprehensive understanding of

charge transport in both vertical and lateral directions, the current-voltage measurement in

tip-plane geometry, similar to c-AFM, is preferred.

To study the morphological eﬀect of crystallized P3HT ﬁbers, we carried out the nu-

118

merical calculations using the fully 3-D DD-SCLC model, described in section 3.5, that is

capable of incorporating the three-dimensional morphology of nanoﬁbers as well as handling

the device with tip-plane electrode conﬁguration. The calculations were performed on the
cubic grid, where the conducting tip is approximated by a 5nm × 5nm square contact. In

the following theoretical study discussed in this section, two assumptions have been made:

(i) the device is free of traps; and (ii) while the anisotropic mobility is assigned to the ﬁbers,

the background mobility of amorphous regions is isotropic. The parameters used in the

numerical calculations are summarized in table 4.1 unless stated otherwise.

4.3.1 An Isolated Fiber: ‘Face-On’ Fiber and ‘Edge-On’ Fiber

Although it is not common to selectively grow only one P3HT nanoﬁber in the ﬁlm, the

theoretical study of an isolated single ﬁber can give us the fundamental understanding of its

morphological eﬀect on the theoretical current-voltage characteristics. Also, we can investi-

gate the contribution of an individual ﬁber without the interference from the near-by ﬁbers,

which cannot be observed directly from experiments.

To examine the morphological eﬀect of an isolated crystalized P3HT nanoﬁber, we consid-

ered the case of a high mobility ﬁber embedded in an amorphous background. As schemat-

ically shown in ﬁgure 4.14, P3HT nanoﬁbers are depicted by green wires with the cross-
sectional areas of 20nm × 5nm. To maximize the eﬀect, the ﬁber is lengthened across the

entire device, which is about 80nm for ‘face-on’ ﬁber and 600nm for ‘edge-on’ ﬁber. This

is generally consistent with the P3HT nanoﬁbers visualized in experiments [1, 31, 36, 37, 38].

When compared to the amorphous domain, it has been shown that the ordered structures
like nanoﬁbers have at least improved the inter-molecular transport along the π− π stacking

direction [65, 67, 90, 91]. To proceed with the device simulations, the hole mobility in the

119

(a)

(b)

Figure 4.14: Schematic illustration of device structures with the tip-plane electrode con-
ﬁguration to study the morphological eﬀect of an isolated ﬁber with (a) ‘face-on’ and (b)
‘edge-on’ orientations.

nanoﬁber needs to be parameterized.

Evidently, the mobility anisotropy of semi-crystalline P3HT ﬁlms have been observed in
experiments [65, 67, 90, 91]. Hole mobilities along the direction of π − π stacking and alkyl

side chains have been investigated by using ﬁeld-eﬀect transistor experiment [65] and c-AFM
measurement [91]. These papers have reported that hole mobility along the π − π stacking

direction (µπ) is signiﬁcantly larger than that along the side chain direction (µS) by 2-3

orders of magnitudes. On the other hand, the hole mobility along polythiophene backbone

direction (µB) cannot simply be determined by the same techniques. To the best of our

knowledge, there is no report on the experimental measurement of the hole mobility along the

backbone direction, and only one theoretical analysis by Lan and Haung [90] suggested that

the hole mobility along the backbone direction (µB) is approximately 2 orders of magnitude
larger than that along the π − π stacking direction (µπ). This implies that holes travel in

nanoﬁbers along the polymer backbone with mobility of 4 orders of magnitude higher than

that along the polymer side chains. From the computational point of view, we have found

120

that this drastic transport anisotropy in P3HT nanoﬁbers impacts the numerical stability of

the calculations. To solve this issue, we need to evaluate how the high mobility anisotropy

of nanoﬁbers aﬀects the overall hole transport in semi-crystalline P3HT devices.

The straightforward way to examine the overall hole transport in devices is by evaluating

their current-voltage characteristic. For the present study, the conducting tip is located on

nanoﬁbers to maximize their impact on the macroscopic transport, as schematically illus-

trated in ﬁgure 4.14. We carried out extensive calculations of the dependence of theoretical

current-voltage characteristics on the magnitude of hole mobility along the polythiophene

backbone direction of a nanoﬁber (µB). While varying hole mobility along the backbone
direction, hole mobility along the π − π stacking direction and the alkyl side chain direction

were ﬁxed with the anisotropic mobility ratio µπ/µS = 100. The background mobility of

the amorphous region is chosen to be the same as that of the alkyl side chains, which is

consistent with recent experiments [1, 34, 67].

Figure 4.15 shows the comparison of theoretical current-voltage characteristics simulated

from hole-only devices consisting of an isolated P3HT nanoﬁber embedded in an amor-

phous domain, as schematically illustrated in ﬁgure 4.14, for the diﬀerent values of mobility

anisotropic ratio of the polymer backbone to the alkyl side chains (µB/µS). The theoretical

current-voltage characteristic of an amorphous device (no ﬁber) is plotted as the baseline.

The inset shows the same data plotted in double logarithmic scale, obeying the power law
I ∝ V m with the exponent of m = 1.89, which is close to the trap-free value, and con-

sequently incorporating the high transport domain doesn’t change the ﬁtting value of the

exponent.

From ﬁgure 4.15, two diﬀerent ﬁber orientations reveal two distinct dependences of the

theoretical current-voltage characteristic on the mobility anisotropic ratio µB/µS. In the

121

Figure 4.15: Based on the hole-only devices with an isolated P3HT ﬁber shown in ﬁgure
4.14, the theoretical current-voltage characteristics when the tip is in contact with the ﬁber
are plotted with varying the ﬁber mobility along the polythiophene backbone direction (µB).
The theoretical current-voltage characteristic of an amorphous device is plotted as the ref-
erence to demonstrate the current enhancement in the presence of a single nanoﬁber. The
inset shows the same data plotted in the double logarithmic scale.

122

Applied Voltage (V)0200040006000Current (pA)Face-on µB/µS=1Face-on µB/µS=10Face-on µB/µS=100Edge-on µB/µS=1Edge-on µB/µS=10Edge-on µB/µS=100No fiber02468100200400600% Increase0110Applied Voltage (V)110100100010000Current (pA)case of an isolated ‘edge-on’ ﬁber, the theoretical current-voltage characteristics remain un-

changed with the variation of ratio µB/µS, indicating that hole transport in the ‘edge-on’
ﬁber is dominated by the high mobility along the π − π stacking direction. On the other

hand, the current of a device with an isolated ‘face-on’ ﬁber decreases with increasing hole

mobility of the polymer backbone (µB). This setback in hole transport is associated to the

enhancement of the lateral current spreading due to an increase of in-plane mobility along the

backbone direction of the ‘face-on’ ﬁber. It should be noted that we veriﬁed that increasing

the mobility anisotropic ratio µB/µS to 104, as in the analysis by Lan and Haung [90], has a

very small eﬀect on the current-voltage characteristic (< 1% diﬀerence from the case of ratio

µB/µS = 100 ). To reduce the diﬃculty in the numerical calculations without impacting
the mobility anisotropy eﬀect, the hole mobility along the π − π stacking direction and the

polythiophene backbone direction is set to be 100 times larger than that along the alkyl side

chain direction. This is to be applied to all the numerical calculations on P3HT ﬁbers from

this point onward.

To further understand the morphological eﬀect of an isolated nanoﬁber contained in the

theoretical current-voltage characteristics shown in ﬁgure 4.15, we also plot the percentage

increase of current when the tip is placed on an isolated nanoﬁber compared to when the tip

is placed on an amorphous domain (no ﬁber), as a function of externally applied voltage. It

is clear that the overall hole transport in the device is enhanced when tip is in contact with

a high mobility nanoﬁber; approximately 200% and 500% increase of current are observed

for an isolated ‘edge-on’ and ‘face-on’ ﬁber, respectively. This observation demonstrates

a remarkable dependence of the overall hole transport on the orientation of a nanoﬁber

in this particular setup.

In the case of a ‘face-on’ orientation, the high mobility ﬁber is

perpendicular to the substrate, and consequently improves the vertical hole transport from

123

the conducting tip to the bottom electrode, resulting in a tremendous increase of current.

However, the explanation based on the improvement of vertical mobility cannot elucidate

the current enhancement observed when tip is in contact with the ’edge-on’ ﬁber, in which

the high mobility direction is in-plane. Moreover, as previously seen in the case of a ‘face-on’

ﬁber that the increasing of lateral hole mobility results in the decreasing in current, one

might wonder why we do not observe the same setback in the case of edge-on ﬁber when

compares to the case of amorphous P3HT device (no-ﬁber case). With this complication, to

consider just the theoretical current-voltage characteristics is not suﬃcient, and we need to

explore the detailed information on current ﬂow.

Figure 4.16 shows a selection of current maps, demonstrating three-dimensional current

ﬂow simulated from hole-only devices in the tip-plane geometry at the applied voltage Va =

1V . Figure 4.16(a-e) presents the vertical current at diﬀerent locations along the direction of

device thickness, while ﬁgure 4.16(f-j) shows the lateral current at the top surface (z = 80nm)

in contact with the conducting tip. Due to the large amount of current injecting vertically

and laterally from the conducting tip, the tip contact area is visualized as a red dot in

current maps of the top surface. In ﬁgure 4.16, three diﬀerent morphologies associated to

the semiconducting P3HT have been used in the device simulations, including (a,f) a purely

amorphous domain (no ﬁber), (b-c, g-h) an isolated ‘face-on’ ﬁber embedded in an amorphous

domain, and (d-e, i-j) an isolated ’edge-on’ ﬁber embedded in an amorphous domain, which

are schematically shown in ﬁgure 4.14. In the last two morphologies involving nanoﬁbers,

we compare the case when the conducting tip is located on ﬁber and the case when the tip

is near but not in contact to the ﬁber. The current maps clearly show distinct correlation

to their morphologies.

In ﬁgure 4.16(a,f), the amorphous device exhibits the behavior of an isotropic system,

124

showing the circularly symmetric distribution of the hole current. Figure 4.16 (a) shows the

vertical current spreading into the larger cross-sectional area toward the bottom electrode.

This 3-D cone of the vertical current is recognized as the signature of the tip-plane geometry

caused by the lateral current spreading beneath the conducting tip, as obviously seen in

ﬁgure 4.16(f).

In the case of an isolated ‘face-on’ ﬁber (ﬁgure 4.16(b-c, g-h)), we clearly observed the

superposition of two features arising from the nature of tip-plane electrode conﬁguration

and the characteristics of high mobility ﬁbers. First, it exhibits the circularly symmetric

current spreading from the conducting tip shown as the background; second, the number of

standout red dots in the vertical current maps indicate the intense vertical current at the

speciﬁc location, revealing the exceptionally high-transport channel through the ﬁber. When

the tip is located on the ﬁber, the ﬁber simply connect the upper tip to the bottom electrode.

As shown in ﬁgure 4.16(b), the alignment of red dots could be seen as the ‘tip extension’ to

the bottom electrode, which lead to the current enhancement previously observed in ﬁgure

4.15. Even when the tip is not in contact with the ﬁber (ﬁgure 4.16(c,h)), the lateral current

spreading brings holes to the nearby ﬁber where holes are able to travel faster to the bottom

electrode. In other words, the ﬁber assists hole transport by providing a ‘highway’ for current

ﬂow.

The last morphology to be discussed is the one with an isolated ‘edge-on’ ﬁber. In ﬁgure

4.16(d-e,i-j), the spherical symmetry of current spreading from the conducting tip is no longer

maintained, and the current maps clearly characterize the surface morphology of the ﬁber

that lies in-plane. Adopting the intuitive picture of the ‘tip extension’ described earlier, one

can relate the enhancement of the vertical current when the tip is in contact with the ﬁber

(as shown in ﬁgure 4.15) to the increase of the tip contact area, as clearly depicted by the red

125

wire in ﬁgure 4.16(d,i). Furthermore, in the condition that the ﬁber and the tip are no longer

in contact (ﬁgure 4.16(e,j)), as long as the ﬁber is aﬀected by the lateral current spreading

from the tip, it is obvious that the ﬁber promotes current ﬂow, behaving as a ‘highway’ for

hole transport. It is less obvious, though, that this mechanism is suﬃcient to signiﬁcantly

raise the vertical current of the device.

In order to avoid confusion, it is worth noting that the current maps shown in ﬁgure

4.16 are not the same as the current map of c-AFM measurements, despite the similarity

in the electrode conﬁguration. While the simulated current map illustrates the 3-D current

distribution, including both lateral and vertical currents, when the biased tip is positioned

at the point of interest on the surface; the c-AFM current map is the 2-D plot of the

vertical current measured at each location on the map by scanning the biased tip over the

surface of the sample. Although the simulated current map of a single point probe is limited

to the local morphological and electrical properties where the tip is located, these maps

essentially provide the 3-D visualization of current ﬂow which cannot be observed directly

from experiments. Particularly, by comparing the current maps simulated from diﬀerent

morphologies as illustrated in ﬁgure 4.16, we reveal how the ‘face-on’ and ‘edge-on’ ﬁbers,

despite the diﬀerence in ﬁber orientation, enhance the current ﬂow in the tip-plane geometry.

This mechanism is qualitatively consistent with the contrast seen in the c-AFM current maps

of ﬁbrous ﬁlms [1, 31, 36, 37, 38].

126

Figure 4.16: The vertical (a-e) and lateral (f-j) current of the semiconducting P3HT devices when the biasing voltage Va = 1V is
applied to the conducting tip (red dot) and the bottom electrode is grounded. Three diﬀerent types of morphologies have been
used in the device simulation: (a)&(f) an amorphous (no-ﬁber) domain; (b-c)&(g-h) an isolated ‘face-on’ ﬁber embedded in an
amorphous domain; and (d-e)&(i-j) an isolated ‘edge-on’ ﬁber embedded in an amorphous domain. In the last two morphologies,
the comparison of current ﬂows when tip is and is not in contact with the ﬁber is demonstrated.

127

Up to now, the current maps of a single-point probe have been helpful to qualitatively

elucidate the mechanism for c-AFM contrast in ﬁbrous ﬁlms. However, to verify the ability

of c-AFM current maps in visualizing each individual ﬁber, it is necessary to analyze the

contrast of the map quantitatively. Resembling the c-AFM measurement, we carried out

extensive calculations at diﬀerent positions along the chosen proﬁle trace. In this study, we

consider the same morphologies used in simulating the current map, which consists of an

isolated ﬁber embedded in an amorphous domain, as schematically shown in ﬁgure 4.14.

Figure 4.17 shows the dependence of the vertical current as a function of tip position on

a trace crossing the ‘face-on’ ﬁber along; (a) the polythiophene backbone direction and (b)

the alkyl side chain direction, at the applied voltages Va = 0.1, 1, 10V . Shown in the inset

is the surface morphology of the device in ﬁgure 4.14(a) with the proﬁle trace indicated by

the dashed line. Because of the symmetry of the morphology, we plotted the vertical current

collected from one half of the proﬁle trace, starting from the center of the ﬁber where tip

position is set to zero. The vertical currents in ﬁgure 4.17 have been normalized to the

vertical current of an amorphous (no ﬁber) device in order to make the comparison among

the three values of the applied voltage. Despite the diﬀerence in ﬁber dimensions, two obvious

similarities of the current distribution in the backbone direction and the side chain direction

are: (i) the vertical current is at its maximum value when tip position is at 0nm (the center of

ﬁber), (ii) the vertical current abruptly decreases when the tip is positioned outside the ﬁber

edge (FE). The ﬁber edge (FE) is the key parameter that pinpoints the interface between

the crystallized structure of the ﬁber and the amorphous domain, which is at the position

10nm and 2.5nm in the backbone direction and the side chain direction, respectively. To

understand the information contained in these data, the distribution of the vertical current

is parameterized in terms of the half-width at half maximum (HWHM) as marked and

128

(a)

(b)

Figure 4.17: Simulating the c-AFM measurements, the dependence of the vertical current
as a function of tip position on the proﬁle trace crossing the ‘face-on’ ﬁber along; (a) the
polythiophene backbone direction and (b) the alkyl side chain direction, at the applied
voltages Va = 0.1, 1, 10V . The proﬁle trace is indicated by the dashed line shown in the
inset.

129

labeled in the ﬁgure, yielding the averaged values HWHM of 14.0nm and 7.2nm in the

backbone direction and the side chain direction, respectably. To generalize this parameter,
we consider the value of HWHM − FE, which simply indicates the tip position away from

the ﬁber edge where the vertical current decreases to the half value of its maximum. As a
result, HWHM − FE is estimated to be 4.0nm and 4.7nm for the c-AFM simulations along

the proﬁle trace of the backbone direction and the side chain direction, respectively. These

values are comparable to the tip size, and conﬁrm the high lateral resolution of c-AFM

measurements on a ‘face-on’ ﬁber.

Figure 4.18(a) shows the dependence of the vertical current as a function of tip position

on the trace crossing the ‘edge-on’ ﬁber along the polythiophene backbone direction at the

applied voltages Va = 0.1, 1, 10V . In the inset, the proﬁle trace is indicated by the dashed

line cross the surface morphology of the device schematically shown in ﬁgure 4.14(b). Similar

to the case of ‘face-on’ ﬁber, the vertical current is maximized when the tip is positioned at

the center of ﬁber (0nm), then falls oﬀ rapidly when the tip moves beyond the ﬁber edge

(FE) at the position of 10nm. As marked and labeled in ﬁgure 4.18(a), the values of HWHM

evaluated at the three diﬀerent applied voltages are close, giving an average value of HWHM
as 16.0nm. Consequently, HWHM − FE is estimated to be 6.0nm, which is comparable to

the tip size. This conﬁrms the high lateral resolution of c-AFM measurements on a ‘edge-on’

ﬁber.

Figure 4.18(b) illustrates an analysis of the vertical sensitivity of c-AFM measurement on

an isolated ‘edge-on’ ﬁber at the applied voltage Va = 0.1, 1, 10V . To investigate the vertical

sensitivity of the measurement, we ﬁxed the tip position and modiﬁed the morphology by

embedding the ‘edge-on’ ﬁber away from the surface. In Figure 4.18(b), dependence of the

vertical current as a function of the vertical separation between tip and ﬁber is illustrated

130

(a)

(b)

Figure 4.18: Simulating c-AFM measurements, (a) the dependence of the vertical current
as a function of tip position on the proﬁle trace crossing the ‘edge-on’ ﬁber along the poly-
thiophene backbone direction and (b) the dependence of the vertical current as a function
of tip-electrode separation, at the applied voltage Va = 0.1, 1, 10V .

131

for the varying applied voltages. The vertical current decreases with the increase of the

tip-ﬁber separation. This is expected since the eﬀect of the ‘edge-on’ ﬁber subsides as the

ﬁber is located further away from the surface, as well as the tip. Similarly, it is useful to

parameterize the decay of the vertical current with tip-ﬁber separation in term of HWHM,

as given in the ﬁgure 4.18(b). However, in this study, the physical interpretation of HWHM

is the vertical separation between tip and ﬁber that corresponds to reducing the vertical

current by half of the maximum current measured when tip is in contact with ﬁber. The

averaged HWHM over three applied voltages is estimated to be 6.9nm, indicating the vertical

sensitivity of c-AFM measurement to a buried ‘edge-on’ ﬁber.

4.3.2 An ‘Edge-On’ Fiber Network

In this section, we applied the numerical tool of the 3-D DD-SCLC model to study the mor-

phological eﬀect of a nanoﬁbrillar network. In an exemplary way, we carried out the numer-

ical calculations to simulate the charge transport through the network of ‘edge-on’ oriented

P3HT ﬁbers in the tip-plane electrode conﬁguration, in analogy to c-AFM measurements.

In this study, we consider three types of randomness typically observed in semi-crystalline

ﬁber structure in P3HT ﬁlms [1, 31, 36, 37, 38], including random ﬁber spacing, random ﬁber

orientation, and random ﬁber length. Each model morphology is generated to meet four

conditions:

i) no in-plane ﬁber crossing,

ii) no out-of-plane ﬁber crossing,

iii) no violation of the periodic boundary condition,

iv) constant ﬁber fraction in each cross-sectional area.

132

(a)

(b)

(c)

Figure 4.19: The 600nm × 600nm cross-sections of three model morphologies consisting of
‘edge-on’ ﬁbers aligned in-plane with (a) random spacing (ﬁber fraction = 0.40), (b) random
spacing and random orientation (ﬁber fraction = 0.40) and (c) random spacing and random
length (ﬁber fraction = 0.32).

133

Figure 4.20: Theoretical current-voltage characteristics of hole-only devices in the tip-plane
geometry with three model morphologies (see ﬁgure 4.19) simulated when the conducting tip
is in contact with the ﬁber and when it is not, as labeled by ‘ON’ and ‘OFF’, respectively. The
inset shows the same data plotted in the double logarithmic scale. Despite the diﬀerences
in model morphologies, the total currents simulated from the devices with an ‘edge-on’
ﬁber network increase by approximately 50% compared to those simulated from the devices
with a single With the presence of an ‘edge-on’ ﬁber network, the total currents, despite
the diﬀerences in model morphologies and tip locations, increases by approximately 50%
compared to those simulated from the devices with a single ‘edge-on’ ﬁber.

134

Applied Voltage (V)010002000300040005000Current (pA)ON  -40% Random SpacingOFF-40% Random SpacingON  -40% Random Spacing & Orien.OFF-40% Random Spacing & Orien.ON  -32% Random Spacing & LengthOFF-32% Random Spacing & LengthON  -Single FiberOFF-Single fiber0246810050100% Increase0.1110Applied Voltage (V)110100100010000Current (pA)Figure 4.19 illustrates the 600nm × 600nm cross-sections of the x − y plane for three

model morphologies, consisting of the ‘edge-on’ ﬁbers (depicted by green wires) that are

aligned in plane with (a) random spacing, (b) random spacing and random orientation, (c)

random spacing and random length. These ﬁbers are embedded in an amorphous domain,

colored in brown.

Initially, the fraction of ﬁber in each cross-sectional area is chosen to

be 0.40. However, when we assigned the random ﬁber length to the morphology(c), the

ﬁber fraction in each cross-sectional area is reduced to 0.32 instead. These values of ﬁber

fraction are roughly consistent with the value evaluated from c-AFM current maps in the

experiment [1].

In analogy to c-AFM measurements, the local variation in hole transport of the model

morphologies is achieved by locating the tip at the spot of interest and generating a theoret-

ical current-voltage characteristic. Figure 4.20 shows the comparison of theoretical current-

voltage characteristics simulated from two distinct locations on the top surface (z = 80nm)

of the device, i.e., ON- and OFF- ﬁber, as marked on the surface morphologies in ﬁgure

4.19. The term ‘ON’ and ‘OFF’ are simply referred to the measurement conditions when the

c-AFM tip is in contact and not in contact with ﬁber, respectively. The theoretical current-

voltage characteristics of ON- and OFF- an isolated ‘edge-on’ ﬁber, as studied in section

4.3.1, are plotted as the baselines. Note that the data for OFF-an isolated ﬁber is equivalent

to the data generated from an amorphous device. The inset shows the same data plot in
double logarithmic scale, described by the power law I ∝ V 1.89, which is close to the char-

acteristic of a trap-free no-ﬁber device. By plotting the current-voltage characteristic this

way, we can make two qualitative observations: (i) the current ﬂow is signiﬁcantly enhanced

by the network of high-mobility ﬁbers; (ii) the diﬀerence in the theoretical current-voltage

characteristics simulated from the three model morphologies is insigniﬁcant.

135

To quantitatively evaluate the morphological eﬀect of a nanoﬁbrillar network shown in

ﬁgure 4.20, we also plotted the percent increase of current generated from three model

morphologies compared to the baseline (an isotropic ﬁber) as a function of applied volt-

age. Despite the diﬀerences in model morphology, ﬁber fraction, as well as the measuring

condition, we observed approximately 50% increase of current in all devices consisting of a

nanoﬁbrillar network when compared with devices consisting of an isolated nanoﬁber.

4.4 Conclusions

In this chapter, we addressed the theoretical studies of DD-SCLC model for the description

of hole transport in hole-only devices. Applying the numerical tools developed in chapter

3, we demonstrated how the hole transport in the devices is aﬀected by charge diﬀusion,

energetically distributed traps and nanoscale morphologies.

First, we investigated the eﬀect of charge diﬀusion on the SCLC in the conventional planar

geometry, which is similar to the electrode conﬁguration of the simpliﬁed SCLC theory, by

using the 1-D DD-SCLC model.

In contrast to those of the simpliﬁed SCLC model, the

theoretical proﬁles of electric potential and hole carrier density simulated from the DD-

SCLC model shows the strong dependence on the externally applied voltage, exhibiting

three transport regimes: (i) diﬀusion-dominant transport regime, (i) drift-diﬀusion assisted

transport regime, and (iii) drift-dominant transport regime, along the device thickness from

the injecting electrode (an actual anode) to the extracting electrode (an actual cathode). The

three transport regimes are separated by two major transitions of the local charge transport

occurring at the position of the potential maximum, where the local electric ﬁeld and drift

current vanish, and the position of the hole density minimum, where the diﬀusion current

136

becomes zero. These conditions at the turning points are consistent with the boundary

conditions at the electrodes of the simpliﬁed SCLC theory, marking the location of the

virtual anode and a virtual cathode, respectively. We also found that the positions of the

virtual electrodes move toward those of the actual electrodes with increasing applied voltage.

The dynamics of the virtual electrodes and the corresponding transport regimes demonstrate

the crucial eﬀect of charge diﬀusion and clarify the transport mechanism near the electrodes.

Furthermore, we carried out the numerical calculation of the devices in tip-plane geometry

by using the 2-D DD-SCLC model. Similar to the planar geometry, the vertical transport

along the vertical path at the center of the conducting tip exhibits three transport regimes

and two virtual electrodes. However, due to the diﬀerent size of injecting and extracting

contacts, i.e., the actual anode and actual cathode, an asymmetric behavior of the virtual

electrodes is observed particularly at the low voltage regime.

Next, we examined the eﬀect of traps on the SCLC. In both simpliﬁed SCLC model and

DD-SCLC model, trap states are assumed to be an exponential distribution of the form

(3.4), characterizing by two parameters: (i) trap density (Nt) and (ii) trap temperature (Tt).
By assuming the power relation J(or I) ∝ V m, the theoretical current-voltage characteristic

simulated from the DD-SCLC model exhibits three distinct regimes corresponding to Ohm’s

law (m = 1), trap-limited SCLC (m > 2) and trap-ﬁlled SCLC (m = 2), ordered from low

to high applied voltage. Trap analysis is focused in the trap-limited SCLC regime where
trap states are partly ﬁlled, which is typically in the range of applied voltage 1 − 10 V. We

have found that the ﬁtting exponent m depends on both trap density and trap temperature.

However, in the strong trapping limit where the trap temperature is high and the number

of traps is large, the ﬁtting exponent depends only on the trap temperature, as previously

described by the simpliﬁed SCLC theory in equation (2.27). Over the range of the scaled

137

trap parameters ft(Nt) = 00.05 and rt(Tt) = 0.200.40, the planar geometry yields the ﬁtting
exponent in range m = 1.89 − 5.00; and the tip-plane geometry yields the ﬁtting exponent
in range m = 2.00 − 5.73. Particularly in the tip-plane geometry, the mobility anisotropy

does not aﬀect the exponent m. This makes the exponent m a suitable ﬁtting parameter for

trap characterization in the SCLC regime.

Lastly, we looked into the eﬀect of the nanoscale morphologies on SCLC in the c-AFM

geometry. In this study, the semi-crystalline P3HT thin ﬁlms, including two common ﬁber

orientations: (i) ‘face-on’ and (ii) ‘edge-on’ ﬁbers, are used as the representative systems to

demonstrate the eﬀect. The 3-D DD-SCLC model successfully models the nanoﬁbrillar mor-

phologies of semi-crystalline materials, which are composed of isolated nanoﬁbers embedded

in an amorphous domain, demonstrating a strongly ﬁlamentary structure of the current ﬂow

along the percolation path of the high mobility ﬁbers, leading to a signiﬁcant current en-

hancement when the conducting tip is in contact with a ﬁber. The device simulations take

the mobility anisotropy of the ﬁbers into account, showing the strong dependence of the hole

transport enhancement on the ﬁber orientation, such that the current simulated when the

tip is located on top of the ‘face-on’ ﬁber is approximately two times larger than that of

the ‘edge-on’ ﬁber. This quantitatively conﬁrms the current enhancement observed in the

c-AFM measurements.

Moreover, we found that the lateral and vertical electrical responses of an isolated ﬁber

are very local and over a distance comparable to the tip size, conﬁrming the high resolution

of c-AFM measurements. A clear example of this is shown in the study of the more complex

morphologies of the ‘edge-on’ ﬁber network embedded in an amorphous domain. We pre-

sented three model morphologies with variations in the ﬁber spacing, ﬁber orientation and

ﬁber length. Despite the diﬀerences in their nanoscale morphologies, the theoretical current-

138

voltage characteristics simulated from all three morphologies are consistent with each other,

demonstrating the local dependence on tip location rather than the overall morphologies.

Additionally, we found that the total current simulated from the case of a ﬁber network

increases by approximately 50% compared to the case of an isolated ﬁber, owing to the

increasing percolation path along the ﬁber network.

139

Chapter 5

The Drift-Diﬀusion SCLC Transport

Models - Application to an

Experimental System

The hole-only device we model in this chapter is based on the experiments conducted by

Jiebing Sun of Dr. Pengpeng Zhang’s group at Michigan State University [1]. In the ex-

periment, the eﬀect of a thermal annealing process on P3HT thin ﬁlm morphology was

investigated using c-AFM measurement.

It was reported that the crystallinity of P3HT

was signiﬁcantly improved by thermal annealing, based on the comparison of c-AFM im-

ages of topography and current maps collected from annealed and non-annealed samples.

This is accompanied by the current enhancement observed from the local current-voltage

measurement when c-AFM tip is located precisely on top of a crystallized P3HT ﬁber. Fur-

thermore, thermal annealing also improved the conductivity of P3HT thin ﬁlms by reducing

traps. This is conﬁrmed by the analysis of the current-voltage characteristics of annealed

and non-annealed samples.

In order to explain the experimental results, we apply our theoretical models, as described

in chapter 3, as well as the trap and morphology analyses presented in chapter 4. First,

following the systematic investigation of trap eﬀects on hole-only devices presented in section

140

4.2, we analyze the c-AFM current-voltage data of annealed and non-annealed samples and

also determine the eﬀective hole mobility in the vertical direction. Secondly, the fully 3-D

DD-SCLC model was employed to study the transport mechanism leading to c-AFM contrast

of annealed samples. The morphology of the ﬁlm is approximated by the three-dimensional

model morphology, as outlined in section 4.3.

5.1 Experiments and Results

According to the publication by Jiebing, et al. [1], the hole-only devices were fabricated

by the spin-coating process. Transparent indium-tin-oxide (ITO) was used as the bottom

electrode of the device. The solution of poly(3,4-ethylenedioxythiophene)/poly (styrenesul-

fonate) (PEDOT:PSS) was ﬁrst spin-coated on top of the ITO substrate, followed by the

spin-coating of the P3HT solution on top of the PEDOT:PSS ﬁlm. The annealed samples
were later completed by heating at 160◦ C for 20 min. A platinum (Pt) coated c-AFM tip

was selected as the top electrode. The tip contact diameter was estimated to be 12nm.

In the current-voltage measurement, the Pt tip was grounded and the biasing voltage was

applied to the ITO electrode.

The authors reported a very intriguing feature in P3HT thin ﬁlms after the thermal an-

nealing process. Unlike those collected from the non-annealed sample, the c-AFM images of

topography and current map collected from the annealed sample clearly show the dense ﬁbers

of crystallized P3HT, suggesting that these ﬁbers were developed through the annealing pro-

cess. It was found that these ﬁbers have width of 20 nm and length of several micrometers.

The combined techniques of X-ray diﬀraction (XRD) and ultraviolet-visible (UV-Vis) spec-

troscopy were used to verify that P3HT formed a crystalline nanoﬁbrillar structure with

141

.

Figure 5.1: Reported in ref. [1], (a) an image of the topography of the annealed sample for
ﬁxed-spot current-voltage (I-V ) measurement; (b) I-V spectra with the tip located on/oﬀ
the ﬁbers at locations labeled in (a). The inset of (b) shows log-log plots of the I-V spectra
(Log|I| vs. Log|V |) from the negative bias regime. (c) Typical I-V curves in the negative
sample bias for annealed samples (blue) and non-annealed samples (green), and the ﬁtting to
power laws IαV m (red). For annealed samples, the I-V spectra have three distinct segments:
Ohms law ( m = 1), space-charge limited transport with shallow traps (m = 2− 2.25), and a
higher slope regime (m > 5) with increase of the applied bias. In non-annealed samples, trap
eﬀects are pronounced in the SCLC regime as the exponent is signiﬁcantly higher (m = 4−5).

142

‘edge-on’ orientation as the primary mode.

With the presence of the ‘edge-on’ ﬁber network, the local variation of electrical conduc-

tivity has been determined by c-AFM measurement. Figure 5.1(a) shows the topographic

image of the annealed sample marking position A as on a crystallized ﬁber and position B

as on an amorphous region (oﬀ ﬁber). In ﬁgure 5.1(b), the corresponding current-voltage

responses measured at position A and B exhibit a distinct behavior under positive and neg-

ative bias. Under the negative device bias, holes are injected into the device, transport

through the P3HT layer, and result in the negative current. The current diminishes as the

negatively applied voltage is reduced, and then vanishes when device undergoes the positive

bias, where charge transport is prohibited. Since the positively-biased voltage regime has

no physical importance, we will focus on the negative device bias case only. The inset of

ﬁgure 5.1(b) presents the Log-Log representation of the absolute values of current-voltage

responses in the negatively-biased SCLC regime, showing that c-AFM current is much larger

when the measuring tip is located on ﬁber, even though the ﬁber is of ‘edge-on’ orientation

which exhibits high in-plane mobility. This observation is well-explained by the theoretical

study discussed in section 4.3.

More information about the electrical transition of P3HT thin ﬁlm through the ther-

mal annealing process was also investigated by comparing the charge transport of annealed

samples with that of non-annealed samples. To characterize the statistical variation in trans-

port properties across the surface of annealed and non-annealed samples, a number of local

current-voltage measurements at diﬀerent locations on each sample were performed. Figure

5.1(c) shows the typical current-voltage responses of annealed and non-annealed samples,

respectively. Quantitatively, at the same applied voltage, the current of annealed samples

is much greater than that of non-annealed samples. Assuming the current-voltage relation

143

as I ∝ V m, the current-voltage characteristics of hole-only devices exhibit three diﬀerent

regimes based on the ﬁtting exponent m: (a) Ohm’s law regime, (b) trap limited SCLC

regime, and (c) anomalous high slope regime, ordered from low to high applied voltage.

Particularly in the intermediate voltage regime where traps aﬀect hole transport, the ﬁtting
of annealed samples yields the exponent m close to 2 (m = 2− 2.25). Compared to annealed

samples, non-annealed samples result in a much larger value of the exponent m typically in
the range m = 4−5. It is clear from the discussion of the trap eﬀects on current-voltage char-

acteristic in section 4.2 that the much steeper dependence of current on voltage is evidence

of the higher density of deep traps. The comparison of the c-AFM current-voltage responses

measured from annealed and non-annealed samples suggests that thermal annealing also

enhances hole transport in the P3HT ﬁlm by reducing traps.

5.2 Theoretical Simulations and Analyses

Based on the fabrication method reported in the publication, the simple schematic repre-

sentation of the device structure and the corresponding energy level alignment is illustrated

in Figure 5.2. The energy alignment diagram of the interface Pt/P3HT/PEDOT:PSS/ITO,

as shown in ﬁgure 5.2(b), is based on the values of work functions and HOMO/LUMO level

of P3HT given in the literature [92, 93]. It is also important to note that the actual HOMO

levels of P3HT and PEDOT:PSS do not exactly match, and a hole injection barrier at this

interface was reported [92]. By using ultraviolet photoelectron spectroscopy (UPS), the hole

injection barrier was measured to be 0.2 eV, and a much larger interface dipole of -1.0 eV

was observed [92]. These two values suggest that fermi-level pinning occurs at this interface,

and the potential barrier is negligible.

144

Figure 5.2: The schematic representation of the device structure (a) and the corresponding
energy alignment (b) of the hole-only device. This device behaves similarly to a semicon-
ductor diode. When the device is positive biased, no charge carrier transport presents. Con-
versely, under the negative bias, holes are injected from Pt probe into the system, traveling
through the layer of P3HT and PEDOT:PSS, and collected at the ITO contact; while elec-
trons injected from ITO cannot transport into the system due to the layer of PEDOT:PSS.
Therefore, this is the hole only injection system.

145

The device design brings an essential diﬀerence between positive- and negative- bias.

When the device is positively biased, both electron and hole transports are prohibited. In

contrast, under the negative bias, current can be observed. Holes are injected from Pt tip

into the device, traveling through the layer of P3HT, and later collected at the ITO electrode.

Electrons, on the other hand, are injected from ITO electrode, but cannot transport into the

P3HT layer due to the PEDOT:PSS electron blocking layer. Therefore, this device allows

us to investigate the single carrier transport of holes in P3HT when applying the negative

bias. These behaviors of the device under the positive and negative biases were conﬁrmed

by the current-voltage data taken by c-AFM measurement, as previously shown in ﬁgure

5.1(b), which is the characteristic of hole-only devices.

The numerical calculations of the DD-SCLC model, as previously explained in chapter

3, were performed to study hole transport in the device. However, one minor adaptation

has been made on the electrode boundary conditions to particularly match the behavior

of the fabricated P3HT devices under the bias from c-AFM tip, assuming the perfect hole

extraction at the bottom electrode. The modiﬁed boundary conditions are written as

ψ = 0 and pf = Nv

at the c-AFM tip

(5.1a)

ψ = −V

and pf = 0

at the bottom (planar) electrode

(5.1b)

This is similar to the boundary conditions used in the simulation of hole-only device by Reid

et al. [32]. Unless stated otherwise, the parameters used in the device simulations presented

in this chapter are summarized in table 5.1.

146

Table 5.1: Summary of parameters used in the device simulations to analyze c-AFM data

Parameter

Symbol Numerical value

Dielectric constant
The eﬀective density of states in the HOMO
Tip contact diameter
Temperature
Device thickness

εr
Nv
D
T
L

3.0
1021 cm−3
12 nm
300 K
80 nm

5.2.1 Trap and Mobility Anisotropic Eﬀects

Evidently, the current-voltage characteristics of hole-only devices obey power laws. In the

low and intermediate applied voltage regimes, the c-AFM current-voltage data are in good

agreement with the theoretical current-voltage relations simulated using the 2-D DD-SCLC

model with tip-plane geometry, and hole transport behaviors can be described by Ohm’s law

and the theory of SCLC in the presence of traps, respectively (See section 4.2.2 for detailed

discussion). However, in the high applied voltage regime, one can observe a much steeper

increase of current with applied voltage than predicted by the model. This discrepancy is

possibly attributed to the mobility enhancement due to the electric ﬁeld, so called ‘Poole-

Frenkel’ mobility [32], which is not integrated in the model.

In this section, the systematic method to study trap eﬀects in hole-only devices presented

in section 4.2 is used to analyze the c-AFM current-voltage responses in the intermediate

applied voltage regime, which is typically between 0.5V and 5V , and also to determine the

eﬀective hole mobility in the vertical direction of P3HT samples. First, extensive device sim-

ulations over the wide range of traps and a considerable range of hole mobility anisotropy

were performed using the parameters summarized in table 5.1. To accommodate the spread-

ing eﬀect, we consider the calculation domain as a circular cylinder of 750nm radius and

147

80nm thickness. The c-AFM tip is approximated by a circular contact of 6nm radius. In

the model, traps are described by the energy states distributed exponentially within the

band gap, which are characterized by parameters ft = Nt/Nv and rt = T /Tt where Nt is

trap density and TT is trap temperature. The mobility anisotropy is reﬂected by the ratio

of the in-plane (lateral) mobility to the out-of-plane (vertical) mobility (µρ/µz). The mo-

bility anisotropic ratio ranges from 1 (isotropic case) to 25, which is the strongest mobility

anisotropy observed in P3HT thin ﬁlms composed of semi-crystalline whiskers [67].

Examples of theoretical current-voltage characteristics of hole-only devices in tip-plane

electrode conﬁguration are illustrated in ﬁgure 5.3. Figure 5.3(a) shows the dependence

of theoretical current-voltage relations on mobility anisotropy in the low deep trap regime

(ft = 0.001, rt = 0.30); ﬁgure 5.3(b-c) show the current-voltage relations simulated with an

isotropic mobility when (b) ft = 0.005 with varying rt and (c) rt = 0.35 with varying ft.

These data are plotted in the intermediate applied voltage regime involving the trap ﬁlling

process in space-charge limited transport, which can be described by the semi-empirical

expression (4.4). Afterwards, the ﬁtting of all simulated current-voltage relations yield the

ﬁtting parameter kI and the ﬁtting exponent m, summarized in table A.3 of appendix A.

Figure 5.4 shows the contour plot of the exponent m as a function of trap parameters ft and

rt. The exponent m varies from 1.96 to 5.52 corresponding to trap parameter variations of
ft = 0− 0.05 and rt = 0.20− 0.40. Next, we will present how to determine the eﬀective hole

mobility in the vertical direction through the analysis of the experimental data.

To determine the eﬀective hole mobility in the vertical direction represented by µo in

the semi-empirical expression (4.4), we started by ﬁtting the experimental current-voltage

relations in the intermediate applied voltage regime with a power law, then obtained the
values of the coeﬃcient µo · kI and exponent m. Since the exponent m is ﬁxed for each

148

Figure 5.3: Based on the setup of c-AFM measurements [1], the theoretical current-voltage
characteristics in the SCLC regime are well-described by the semi-empirical expression (4.4)
when (a) trap parameters is ﬁxed (ft = 0.001 rt = 0.30) with varying mobility anisotropic
ratios µρ/µz; (b) mobility is isotropic and ft = 0.005 with the varying rt yielding the
exponent in range m = 2.31 − 5.19; and (c) mobility is isotropic and rt = 0.35 with varying
ft yielding the exponent in range m = 1.96 − 3.61.

149

Figure 5.4: Dependence of the exponent m on the parameters ft and rt for hole-only de-
vices in the SCLC regime, resulting from extensive calculations of theoretical current-voltage
characteristics demonstrated in ﬁgure 5.3 [1].

150

current-voltage curve, the values of ft and rt can only be taken along a contour line of the

plot in ﬁgure 5.4 that has the same value of m. We then compiled a table of kI values with

respect to ft, rt and µρ/µz, and deduced the mobility µo from the ﬁtting coeﬃcient.

Let us begin with the analysis on c-AFM current-voltage data of the annealed samples. In

the experiment, the exponent m found from current-voltage responses at diﬀerent locations

in the ﬁlms exhibits a signiﬁcantly larger variation as compared to the local variation of

current-voltage responses when the c-AFM tip is located on and oﬀ a ﬁber. This is a clear

evidence that the transport still is dominated by local traps in the amorphous regimes. Using

an average mobility anisotropy is then a reasonable approximation of ﬁbrous morphology

in characterizing transport of the semi-crystalline P3HT ﬁlm. Figure 5.5 illustrates how

the theoretical current-voltage characteristics apply to the experimental results of annealed

samples. Figure 5.5(a-b) show a representative example of c-AFM current-voltage response

measured from annealed P3HT ﬁlms, i.e., both data (black dots) are identical. The dashed

red lines are ﬁts using theoretical current-voltage characteristics simulated from two distinct

pairs of trap parameters that result in the same value of m: (a) ft = 0.0005 and rt = 0.35,

(b) ft = 0.001 and rt = 0.40. Clearly, agreement with the experimental data is excellent.

Additionally, table 5.2 lists four representative pairs of trap parameters ft and rt that provide

a good ﬁt to the current-voltage response. As the result, the eﬀective mobility demonstrates

a weak dependence on trap parameters ft and rt, while it is strongly dependent on the

mobility anisotropy, exhibiting approximately an order of magnitude diﬀerence from the

lowest to the highest values of anisotropic ratio. Therefore, the vertical mobility of each

current-voltage curve is approximated by averaging the values ﬁtted over the range of valid

trap parameters.

By carrying out the extensive analysis of a total of 48 c-AFM current-voltage responses

151

Figure 5.5: The representative current-voltage (I-V ) response from the annealed sample
(solid circle) is ﬁtted to the theoretical I-V relation (solid red line) in the low trap regime
with two pairs of trap parameters that result in the same slope: (a) ft = 0.0005 and rt = 0.35,
(b) ft = 0.001 and rt = 0.40. Shown in the inset are the same plot in linear scale. (c) A
histogram of the normalized vertical hole mobility extracted from 48 I-V datasets on the
annealed P3HT ﬁlm, with the comparison of an isotropic model and an anisotropic model
with µr/µz = 25. The mobility is normalized to the maximum value of each case.

152

of annealed samples, we found that the exponent m is ﬁtted in the narrow range of 2− 2.25,

suggesting that annealed samples exhibit the low deep trap regime with trap parameters
ft = 0 − 0.001 and rt = 0.30 − 0.40. Figure 5.5(c) shows the histogram of the vertical mo-

bilities found from the 48 locations on annealed ﬁlms when assuming isotropic mobility and

anisotropic mobility of µρ/µz = 25. Regardless of the mobility model, the eﬀective mobilities,

which are normalized to the maximum, demonstrates the same distribution. This scalable

property of mobility anisotropy is a signature of low deep trap regime, resulting from the lin-

ear dependence of the ﬁtting parameter kI on the mobility anisotropic ratio. Unfortunately,

we cannot simplify the semi-empirical expression (4.4) further since this linear relation be-

tween kI and µρ/µz still depends on trap parameters. Varying the mobility anisotropic ratio

yields the eﬀective mobility in the vertical direction, as summarized in table 5.3. Assuming
isotropic mobility, we found the eﬀective mobility of 2.01 ± 1.46 × 10−4cm2V −1s−1. This

is consistent with the values of hole mobility measured experimentally [32, 35, 94], validat-

ing the methodology used in this study. However, it is inevitable to raise the question on

the magnitude of mobility anisotropy. The eﬀective mobility extracted from the anisotropic
model, especially when µρ/µz (cid:62) 10, is signiﬁcantly lower than what has been reported. We

speculate that semi-crystalline P3HT ﬁlms produced by thermal annealing may not exhibit

anisotropic conductivity as strong as that composed of whisker-type nanoﬁbers. To prove

this statement, further measurements of hole mobility in the device geometry are needed.

Compared to annealed samples, the extensive c-AFM current-voltage measurements at

diﬀerent locations on non-annealed samples yield a signiﬁcantly steeper exponent m in the
range of 4−5, indicating a higher density of deep traps. From the contour plot of the exponent

m shown in ﬁgure 5.4, one can ﬁnd the corresponding trap parameters in the range of ft =
0.0025 − 0.05 and rt = 0.20 − 0.25. In addition, the authors reported no clear visualization

153

Figure 5.6: The representative current-voltage response from the non-annealed sample (solid
circle) is ﬁtted to the theoretical current-voltage relation (solid red line) in the low trap
regime with two pairs of trap parameters that result in the same slope: (a) ft = 0.0075 and
rt = 0.25, (b) ft = 0.01 and rt = 0.25. Shown in the inset are the same plot in linear scale.

154

Table 5.2: The eﬀective mobilities in the vertical direction found from ﬁtting a representative
examples of c-AFM current-voltage responses measured from an annealed P3HT sample for
a range of anisotropic ratios.

Table 5.3: An analysis of point-to-point variation in transport properties across the ﬁlm is
presented for c-AFM data taken at 48 locations on the surface of an annealed sample. The
average, standard deviation, minimum and maximum of the eﬀective mobility in the vertical
direction found from ﬁtting the experimental data for a range of anisotropic ratio.

Table 5.4: The eﬀective mobilities in the vertical direction found from ﬁtting a representative
example of c-AFM current-voltage responses measured from a non-annealed P3HT sample

155

of crystallized P3HT nanoﬁber observed in the c-AFM images of non-annealed samples.

Without the morphological origin of anisotropic hole transport, we assume hole mobility in

non-annealed ﬁlms to be isotropic. In ﬁgure 5.6, a representative example of current-voltage

response was ﬁtted with two theoretical current-voltage characteristic (dashed red lines)

individually simulated from two pairs of trap parameters that result in the same value of

m: (a) ft = 0.0075 and rt = 0.25, (b) ft = 0.01 and rt = 0.25. Clearly, agreement with

the experimental data is excellent. Consequently, table 5.4 lists three representative pairs

of trap parameters ft and rt that achieve a good ﬁt to the current-voltage response. Unlike

the annealed sample, the vertical mobility is strongly dependent on trap parameters ft and

rt. It was in fact revealed that to base the estimate of vertical mobility on the ﬁtting of

experimental current-voltage relations led to unrealistic results. Characterization of trap

parameters using thermally stimulated current [95, 96], transient spectroscopy [97, 98, 99,

100,101], or impedance spectroscopy techniques [44,102,103,104,105] will be required. If we

assume that the hole mobility of non-annealed samples is similar to that of annealed samples.

the representative example of c-AFM current-voltage response measured from non-annealed

samples can be described using trap parameters ft = 0.006 and rt = 0.25, corresponding

to signiﬁcantly more concentration of traps from the deeper states of the DOS distribution

when compared to an annealed sample. Evidently, the thermal annealing process enhances

the current ﬂow in the ﬁlm by reducing traps.

5.2.2 Traps and Morphological Eﬀects

The assumption of the cylindrical symmetry used in the analysis of c-AFM data is helpful to

elucidate the anisotropic aspect of the carrier transport without introducing the complexity

of fully three-dimensional calculations. However, to examine the ability of c-AFM to resolve

156

the semi-crystalline ﬁber structures observed in annealed P3HT sample, it is necessary to

consider the fully 3-D DD-SCLC model presented in section 3.5. To illustrate the eﬀect, we

carried out device simulations, incorporating the nanoscale morphology of a high mobility

ﬁber network embedded in a low mobility background.

In the numerical calculations, we adopted a three-dimensional device domain of 600nm×

600nm in cross-section and 80nm in thickness, and modeled the c-AFM tip as a square
contact of the size 10nm × 10nm. To illustrate the morphological eﬀect observed in c-

AFM measurements, we considered the case of a high mobility ﬁber network embedded

in a low mobility background. In the study of ‘edge-on’ P3HT ﬁber network presented in

section 4.3, we observed a consistent theoretical current-voltage characteristic across three

model morphologies with a variety of random ﬁber spacing, random ﬁber orientation and

random ﬁber length. In this study, we thus choose the simplest morphology with a single

randomness of the in-plane spacing between ﬁbers. Similar to assumptions made in the

preliminary study, both in-plane and out-of-plane ﬁber crossing are forbidden while periodic

boundary conditions are maintained on the four sides of the device that are not connected
to the electrodes. Figure 5.7 shows the 600nm × 600nm cross-section of the representative

morphology at the vertical position z. As depicted by green wires, each crystallized P3HT

ﬁber has the size of 20nm width and 5nm thickness based on the literature [1,31,36,37,38,87].

The background region in brown is then classiﬁed as amorphous P3HT. The fraction of ﬁber

in the cross-section area is chosen to be 0.33, which is roughly consistent with the ﬁber

density evaluated from c-AFM images of the annealed samples. One of the key aspects

observed in ‘edge-on’ ﬁbers is the mobility anisotropy. In the device simulations, the lateral
mobility along the π − π stacking direction ([0 1 0]) and in the polythiophene backbone

direction ([0 0 1]) are taken to be 100 times larger than the vertical mobility along the alkyl

157

Figure 5.7: The 600nm × 600nm cross-sections of the representative morphology at the
vertical position z in the device thickness direction demonstrating the crystallized P3HT
ﬁbers of 20nm width (green wires) aligned with random spacing in the amorphous region
(brown area). The ﬁber fraction of cross-section surface area is 0.33. The bottom electrode is
placed at z = 0nm, and the locations of tip contact are marked on the z = 80nm cross-section
for the on-ﬁber and oﬀ-ﬁber c-AFM simulations.

158

Figure 5.8: Adapted from [1], the on-ﬁber and oﬀ-ﬁber current-voltage responses of the
annealed P3HT sample are well described by the theoretical current-voltage characteristics
simulated from the representative 3-D morphology of an ‘edge-on’ ﬁber network shown in the
inset on the bottom-right corner. (see ﬁgure 5.7 for cross-section morphology) To described
the anisotropic transport in ﬁbers (green), the in-plane mobility is 100 times larger than the
vertical mobility and the background mobility of amorphous region (brown). The background
mobility is ﬁtted to be 7.78 × 10−5cm2V −1s−1 with trap parameters ft = 1.00 × 10−3 and
rt = 0.30. In the inset on the top-left corner, the vertical current ﬂow when the c-AFM tip
is located on the ﬁber is visualized.

159

side chain direction ([1 0 0]) [65, 67, 90, 91]. On the other hand, the hole transport in the

amorphous region is assumed to be isotropic with the mobility that is the same value as that

of alkyl side chains [34, 67].

To simulate the c-AFM measurements, we performed the device calculations when the

tip is located on the ﬁber and on the amorphous region (oﬀ the ﬁber). The tip locations

are indicated on the cross-section image at z = 80nm, as shown in ﬁgure 5.7. Figure 5.8

illustrates how the theoretical current-voltage characteristics apply to the experimental on-

ﬁber and oﬀ-ﬁber current-voltage responses. The dashed red lines are ﬁts using on-ﬁber

and oﬀ-ﬁber current-voltage characteristics simulated using the trap parameters ft = 0.001
and rt = 0.30, yielding the isotropic mobility of amorphous phase 7.78 × 10−5cm2V −1s−1.

Clearly, agreement with the experimental data is excellent. Evidently, a strong enhancement

of current ﬂow can be observed when tip is located on ﬁber, elucidating the contrast visualized

in the c-AFM current map, as discussed in chapter 4.

5.3 Conclusions

In this chapter, we utilized the 2-D and 3-D DD-SCLC models to analyze the current-voltage

responses of c-AFM measurements on hole-only devices fabricated from semiconducting

P3HT. Experimentally, the comparison of typical current-voltage responses measured from

annealed and non-annealed ﬁlms has shown the correlation between the enhancement of c-
AFM current and the thermal annealing process. By assuming the power relation I ∝ V m,
the annealed samples yield an exponent m close to 2 (m = 2 − 2.25) while the non-annealed
samples exhibit a much steeper exponent in the range of m = 4 − 5. The application of

the numerical tool based on the 2-D DD-SCLC model conﬁrms that the thermal anneal-

160

ing process has signiﬁcantly enhanced the hole transport in P3HT ﬁlm by alleviating the

eﬀect of traps. Furthermore, statistical analysis of 48 annealed samples suggest that the
device corresponds to the low trap regime with trap parameter ft = 1.08 ± 0.28 × 10−3 and
rt = 0.35 ± 0.04. By assuming that hole mobility in such devices is isotropic, the eﬀective
mobility in the vertical direction is estimated to be 2.01 ± 1.46 × 10−4cm2V −1s−1. This is

consistent with experiments [32, 35, 94], conﬁrming the validity of the numerical model and

methodology used in the analysis. However, if the device mobility is anisotropic, the eﬀec-

tive mobility in the vertical direction can be much lower. The comparison of current-voltage

data from c-AFM and device measurements is needed to determine the degree of mobility

anisotropy as well as the vertical mobility of holes in the annealed samples.

Additionally, the thermal annealing process improve the crystallinity of P3HT ﬁlms.

The c-AFM image of annealed ﬁlms shows a large spatial variation of the vertical current,

indicating a nanoscale morphology comprised of nanoﬁbers and the lower conductivity back-

ground. Further characterization using X-ray diﬀraction and UV-Vis spectroscopy indicates

the development of self-organized P3HT nanoﬁbers with the ‘edge-on’ orientation as the

primary mode. By incorporating the three-dimensional model morphology of a nanoﬁbril-

lar network, the simulation based on the 3-D DD-SCLC model in tip-plane geometry has

demonstrated the current enhancement when the conducting tip is in contact with the ﬁber,

which visualizes the c-AFM contrast of ﬁbrous ﬁlms.

161

Chapter 6

Conclusions

The simpliﬁed SCLC theory is a well-established model to explain the unipolar transport in

disordered organic semiconductors. However, a shortcoming lies within the assumption that

charge diﬀusion can be neglected. This can lead to the misinterpretation of experimental

data when the eﬀect of charge diﬀusion in the devices is no longer negligible. In this thesis,

we developed numerical approaches that eﬃciently simulate the hole-only SCLC model with

the full description of hole drift and diﬀusion transport mechanisms, i.e., the DD-SCLC

model. In the case of fully 3-D calculations, the numerical model is able to treat inhomo-

geneous systems including spatially varying trap distributions, nanoscale morphologies, and

the tip-plane (c-AFM) geometry. The application of the model to the analysis of c-AFM

experimental data is presented.

Starting in chapter 2, we provided an overview to the fundamental concepts of injection

and transport mechanisms of charge carriers in disordered organic semiconductor devices,

and further delved into the simpliﬁed SCLC theory which describes the drift-dominated

unipolar SCLC. While the earlier work by M.A. Lampert and P. Mark [43] was focused on

electron-only SCLC, we demonstrated the analytical derivation of the simpliﬁed SCLC theory

for holes. As the result, the analytical expressions to describe the relation of hole current to

externally applied voltage for the trap-free devices, a.k.a. the Mott-Gurney equation, and for

the devices with the exponentially distributed trap density are summarized in Table 2.1. For

162

the application to c-AFM measurements, we provided a brief review of the semi-empirical

formula by O.G. Reid et al. [32] which modiﬁed the Mott-Gurney equation to include the

eﬀect of tip geometry. While organic semiconductors are typically energetically disordered,

to the best of our knowledge, the theoretical model to describe trap-limited SCLC in the

c-AFM or tip-plane geometry has never been introduced.

In chapter 3, we introduced the drift-diﬀusion (DD) SCLC model, as in its name, describ-

ing hole-only SCL transport by both drift and diﬀusion mechanisms. The inhomogeneous

systems, including spatially varying trap distributions, nanoscale morphologies, as well as

tip-plane (c-AFM) geometry are also incorporated. As a matter of fact, traps are commonly

found in disordered organic semiconductors, and thus we modelled the trap DOSs for holes

as the tail states of HOMO DOS toward the energy gap, which is well approximated by the

exponential distribution [43, 44]. The DD- SCLC model is governed by three equations, in-

cluding Poisson’s equation, the drift-diﬀusion equation for holes, and the continuity equation

at steady state. To stabilize the numerical calculations, we apply the Scharfetter-Gummel

discretization [54] to the drift-diﬀusion equation, then numerically solve the three governing

equation simultaneously, similar to the numerical scheme used by L.J.A. Koster et al. [52].

We successfully developed the numerical tools for the DD-SCLC model in one-, two-,

and three- dimensional systems. Among the three models, the 1-D DD-SCLC model is the

simplest one and improves on the simpliﬁed SCLC model by including the mechanism of

charge diﬀusion. While the 1-D DD-SCLC model is numerically much more eﬃcient than

the fully 3-D SCLC model, the trade-oﬀ is the possibility to include the spatially distributed

traps, the nanoscale morphologies, and the c-AFM geometry. Unlike the 1-D and 3-D DD-

SCLC model that are discretized in Cartesian coordinate system, the 2-D DD-SCLC model

is set up in a system with cylindrical symmetry. Thus, the 2-D DD-SCLC model is very

163

eﬃcient to simulate the SCLC in the c-AFM geometry with mobility anisotropy in the lateral

and vertical transport. All 1-D, 2-D, and 3-D DD-SCLC models are veriﬁed to be consistent

with each other, and validated to be consistent with the simpliﬁed SCLC theory in the

drift-dominated regime.

The simpliﬁed SCLC theory is only the ﬁrst step to understand the SCL transport in

hole-only device, yet far from complete. We saw in chapter 4 that, while the theory has

made a major assumption that charge diﬀusion is negligible in the devices, the theoretical

hole carrier density is not uniformly distributed, leading to the extremely large amount

of diﬀusion current near the injecting electrode. This discrepancy becomes more severe

with decreasing applied voltage, demonstrating the importance of treating charge diﬀusion

in device modeling. What is its eﬀect on the SCLC? To answer these question, we have

achieved the prerequisite task by successfully developing the numerical tools for the DD-

SCLC model, and the device simulations, as addressed in chapter 4, have revealed a number of

important factors that aﬀect the SCLC, including charge diﬀusion, traps, as well as nanoscale

morphologies.

First, we investigated the eﬀect of charge diﬀusion on the trap-free SCLC. Unlike those

predicted by the simpliﬁed SCLC theory for the hole-only device in planar geometry, the the-

oretical proﬁles of electric potential and hole density simulated by using the 1-D DD-SCLC

model have a non-monotonic distribution along the device thickness, exhibiting three distinct

transport regimes: (i) diﬀusion-dominant transport regime, (ii) drift-diﬀusion assisted trans-

port regime, and (iii) drift-dominant transport regime, parting from the injecting electrode

to the extracting electrode. Two separations of the three transport regimes are determined

by (I) the position of the potential maximum, where the local electric ﬁeld and the corre-

spondent drift current vanish; and (II) the position of the hole density minimum, where the

164

diﬀusion current become zero. These separations are consistent with the boundary condi-

tions at the electrode of the simpliﬁed SCLC theory, and thus we reintroduced these turning

points as the virtual electrodes since they are no longer the actual ones in the DD-SCLC

model. The extensive device simulations over the range of applied voltage have shown that

the positions of the virtual electrodes move toward the actual electrodes with increasing

applied voltage. Moreover, we have found that the diﬀusion current plays an important role

near the actual injecting electrode, as in transport regime (i), by overcoming the reverse

drift current and driving holes away from the contact. Furthermore, the device simulations

in tip-plane geometry have demonstrated similar transport regimes, but with an asymmetric

behavior of the virtual electrodes due to the diﬀerent size of the actual injecting and extract-

ing electrodes. The dynamics of the virtual electrodes and the concurrent transport regimes

demonstrate the crucial eﬀect of charge diﬀusion on the SCLC, and also clarify the transport

mechanisms near the electrodes that can never be explained by the simpliﬁed SCLC theory.

Traps aﬀect the SCLC signiﬁcantly. We found that the theoretical current-voltage char-
acteristics obey the power relation J(or I) ∝ V m, exhibiting three distinct regimes, corre-

sponding to Ohm’s law (m = 1), trap-limited SCLC (m > 2), and trap-ﬁlled SCLC (m = 2),

ordering from low to high applied voltage. To analyze the eﬀect of traps, we focused on the

trap-limited SCLC regime where trap states are partly ﬁlled, typically found in the range
of applied voltage 1 − 10V . Extensive device simulations have shown that the ﬁtting expo-

nent m depends on two parameters of the ﬁeld-independent exponential distribution of trap

DOSs, including total trap density and trap energy level. However, in the strong trapping

limit where the trap temperature is high and the number of traps is large, the ﬁtting expo-

nent depends only on the trap temperature. This strong trap case is previously predicted

by the simpliﬁed SCLC theory, as described in equation (2.27). Particularly in the tip-plane

165

geometry, the mobility anisotropy does not aﬀect the exponent m. This makes the exponent

m a suitable ﬁtting parameter to characterize traps in the SCLC regime.

Lastly, we examined the eﬀect of the nanoscale morphologies on the SCLC in the tip-

plane geometry, in analogy to the c-AFM measurements. The semi-crystalline P3HT thin

ﬁlms, which consist of a mixture of crystalized nanoﬁber and amorphous domain, is chosen as

the reference system to demonstrate the eﬀect. We found that the percolation nature of hole

transport in the disordered organic semiconductor yields a strongly ﬁlamentary structure of

the current along the percolation paths of the high mobility ﬁbers, leading to a signiﬁcant

current enhancement when the conducting tip is in contact with the ﬁbers. Moreover, it

has been shown that the lateral and vertical electrical responses of the nanoﬁber are very

local and in a probed distance comparable to the tip size. These observations quantitatively

conﬁrm the high resolution of c-AFM measurements, supporting c-AFM as a very useful tool

to study the local electrical properties and local morphologies simultaneously.

On the aspect of the application to the experiments, in chapter 5, we utilized the DD-

SCLC model to elucidate the c-AFM measurements on the P3HT based hole-only devices.

The experiments addressed in this chapter were carried out by Professor Zhang’s group at

Michigan State University [1], aiming to investigate the eﬀect of thermal annealing process on

the electrical properties and nanoscale morphologies of the P3HT thin ﬁlms. With the device

simulations of the 2-D DD-SCLC model, we have shown that the thermal annealing process

has signiﬁcantly improved hole transport in P3HT thin ﬁlms by alleviating the eﬀect of

traps, yielding the eﬀective mobility of an annealed device that is consistent with previously

published experiments. Furthermore, the 3-D DD-SCLC model enables the simulation of hole

transport in the semi-crystalline P3HT thin ﬁlm, consisting of a mixture of crystalized P3HT

nanoﬁbers and amorphous P3HT domains. By incorporating the 3-D model morphology of

166

the nanoﬁbrillar network, the device simulation in the tip-plane geometry has evidently

shown the current enhancement when the conducting tip is in contact with the ﬁber, which

demonstrates the c-AFM contrast of ﬁbrous ﬁlms.

While the current state of DD-SCLC model is limited to the numerical solutions, the

model polishes our understanding of hole transport with the full description of drift and

diﬀusion, and elucidates mechanisms missing from the simpliﬁed SCLC theory. The DD-

SCLC model can be signiﬁcantly improved by considering the ﬁeld-dependent hole carrier

density as a boundary condition at the electrodes, based on the expression by J.C. Scott and

G.G. Malliaras [106]. Furthermore, the accuracy of the DD-SCLC model in the high voltage

regime can be enhanced by including the ﬁeld- and carrier-density- dependent mobilities,

the so-called ‘Poole-Frankel’ model [39, 107] and ‘Pasveer’ model [39, 108], respectively, as

well as using the general form of Einstein relation to describe the charge diﬀusion coeﬃcient.

We anticipate that these improvements will lead to more realistic models, shedding light

on the investigation of charge transport in disordered organic semiconductors. To have a

clearer understanding on the subject would be very beneﬁcial, not only for developing future

organic photovoltaics (OPVs), but also other organic semiconductor devices such as organic

light emitting diodes (OLEDs) and organic ﬁeld-eﬀect transistors (OFETs).

167

APPENDIX

168

Appendix A

Fitting Parameters of Theoretical

Current-Voltage Characteristics

169

Table A.1: Summary of parameters m and kJ estimated by ﬁtting theoretical current-voltage
characteristics of planar devices in the SCLC regime (see section 4.2.1 for details of device
simulations) to the semi-empirical expression (4.3) for a range of trap parameters ft and rt.
Note that R-square of regression is high (R2 > 0.9999).

170

Table A.2: Summary of parameters m and kJ estimated by ﬁtting theoretical current-voltage
characteristics in the tip-plane geometry in the SCLC regime (see section 4.2.2 for details of
device simulations) to the semi-empirical expression (4.4) for a range of trap parameters ft
and rt and a range of mobility anisotropy µρ/µz. Note that R2 > 0.9999 for all regressions.
Table continued on the next page.

171

Table A.2: (cont’d)

172

Table A.3: Summary of parameters m and kJ estimated by ﬁtting theoretical current-voltage
characteristics in the tip-plane geometry in the SCLC regime (see section 5.2.1 for details of
device simulations) to the semi-empirical expression (4.4) for a range of trap parameters ft
and rt and a range of mobility anisotropy µρ/µz. Note that R2 > 0.9999 for all regressions.
Table continued on the next page.

173

Table A.3: (cont’d)

174

BIBLIOGRAPHY

175

BIBLIOGRAPHY

[1] J. Sun, K. Pimcharoen, S. R. Wagner, P. M. Duxbury, and P. Zhang. Nanoscale imag-
ing of dense ﬁber morphology and local electrical response in conductive regioregular
poly(3-hexylthiophene). Organic Electronics, 15(2):441–448, 2014.

[2] M. A. Green. Third generation photovoltaics: Ultra-high conversion eﬃciency at low

cost. Progress in Photovoltaics: Research and Applications, 9(2):123–135, 2001.

[3] J. Nelson. The Physics of Solar Cells. London: Imperial College Press, 2003.

[4] A. C. Mayer, S. R. Scully, B. E. Hardin, M. W. Rowell, and M. D. McGehee. Polymer-

based solar cells. Materials Today, 10(11):28–33, 2007.

[5] M.C. Scharber and N.S. Sariciftci. Eﬃciency of bulk-heterojunction organic solar cells.

Progress in Polymer Science, 38(12):1929–1940, 2013.

[6] A. J. Heeger. 25th anniversary article: Bulk heterojunction solar cells: Understanding

the mechanism of operation. Advanced Materials, 26(1):10–28, 2014.

[7] C. Liu, C. Yi, K. Wang, Y. Yang, R. S. Bhatta, M. Tsige, S. Xiao, and X. Gong.
Single-junction polymer solar cells with over 10% eﬃciency by a novel two-dimensional
donoracceptor conjugated copolymer. ACS Applied Materials & Interfaces, 7(8):4928–
4935, 2015.

[8] C. J. Traverse, R. Pandey, M. C. Barr, and R. R. Lunt. Emergence of highly transparent

photovoltaics for distributed applications. Nature Energy, 2(11):849–860, 2017.

[9] A. Bree, D. J. Carswell, and L. E. Lyons. Photo- and semi-conductance of organic
crystals. Part I. Photo-eﬀects in tetracene and anthracene. Journal of the Chemical
Society, 0:1728–1733, 1955.

[10] L. E. Lyons. Photo- and semi-conductance in organic crystals. Part V. Ionized states

in molecular crystals. Journal of the Chemical Society, 0:5001–5007, 1957.

[11] D. Kearns and M. Calvin. Photovoltaic eﬀect and photoconductivity in laminated

organic systems. The Journal of Chemical Physics, 29(4):950–951, oct 1958.

176

[12] H. Kallmann and M. Pope. Electrolytic contacts for photoconductivity measurements.

Review of Scientiﬁc Instruments, 30(1):44–46, 1959.

[13] A. K. Ghosh, D. L. Morel, T. Feng, R. F. Shaw, and C. A. Rowe. Photovoltaic and
rectiﬁcation properties of Al/Mg phthalocyanine/Ag Schottkybarrier cells. Journal of
Applied Physics, 45(1):230–236, 1974.

[14] B. R. Weinberger, M. Akhtar, and S. C. Gau. Polyacetylene photovoltaic and devices.

Synthetic Metals, 4:187–197, 1982.

[15] G. A. Chamberlain. Organic solar-cells - a review. Solar Cells, 8(1):47–83, 1983.

[16] D. Whrle and D. Meissner. Organic solar cells. Advanced Materials, 3(3):129–138,

1991.

[17] C. W. Tang. Twolayer organic photovoltaic cell. Applied Physics Letters, 48(2):183–

185, 1986.

[18] G. Yu, J. Gao, J. C. Hummelen, F. Wudl, and A. J. Heeger. Polymer photovoltaic
cells: Enhanced eﬃciencies via a network of internal donor-acceptor heterojunctions.
Science, 270(5243):1789–1791, 1995.

[19] G. Zhang, K. Zhang, Q. Yin, X.-F. Jiang, Z. Wang, J. Xin, W. Ma, H. Yan, F. Huang,
and Y. Cao. High-performance ternary organic solar cell enabled by a thick active
layer containing a liquid crystalline small molecule donor. Journal of the American
Chemical Society, 139(6):2387–2395, 2017.

[20] W. Zhao, S. Li, H. Yao, S. Zhang, Y. Zhang, B. Yang, and J. Hou. Molecular opti-
mization enables over 13% eﬃciency in organic solar cells. Journal of the American
Chemical Society, 139(21):7148–7151, 2017.

[21] L. S. C. Pingree, M. C. Hersam, M. M. Kern, B. J. Scott, and T. J. Marks. Spatially-
resolved electroluminescence of operating organic light-emitting diodes using conduc-
tive atomic force microscopy. Applied Physics Letters, 85(2):344–346, 2004.

[22] B. Geﬀroy, P. Le Roy, and C. Prat. Organic light-emitting diode (OLED) technology:
materials, devices and display technologies. Polymer International, 55(6):572–582,
2006.

[23] M. Kuik, G.-J. A. H. Wetzelaer, H. T. Nicolai, N. I. Craciun, D. M. De Leeuw, and
P. W. M. Blom. 25th anniversary article: Charge transport and recombination in
polymer light-emitting diodes. Advanced Materials, 26(4):512–531, 2014.

177

[24] H. Sirringhaus, N. Tessler, and R. H. Friend. Integrated optoelectronic devices based

on conjugated polymers. Science, 280(5370):1741–1744, 1998.

[25] H. Sirringhaus. 25th anniversary article: Organic ﬁeld-eﬀect transistors: The path

beyond amorphous silicon. Advanced Materials, 26(9):1319–1335, 2014.

[26] V. Coropceanu, J. Cornil, D. A. da Silva, Y. Olivier, R. Silbey, and J. L. Bredas.
Charge transport in organic semiconductors. Chemical Reviews, 107(4):926–952, 2007.

[27] C. Goh, R. J. Kline, M. D. McGehee, E. N. Kadnikova, and J. M. J. Frchet. Molecular-
weight-dependent mobilities in regioregular poly(3-hexyl-thiophene) diodes. Applied
Physics Letters, 86(12):122110, 2005.

[28] Z. Chiguvare and V. Dyakonov. Trap-limited hole mobility in semiconducting poly(3-

hexylthiophene). Physical Review B, 70(23):235207, 2004.

[29] O. Douhret, L. Lutsen, A. Swinnen, M. Breselge, K. Vandewal, L. Goris, and J. Manca.
Nanoscale electrical characterization of organic photovoltaic blends by conductive
atomic force microscopy. Applied Physics Letters, 89(3):032107, 2006.

[30] A. Alexeev, J. Loos, and M. M. Koetse. Nanoscale electrical characterization of semi-
conducting polymer blends by conductive atomic force microscopy (C-AFM). Ultra-
microscopy, 106(3):191–199, 2006.

[31] M. Dante, J. Peet, and T. Q. Nguyen. Nanoscale charge transport and internal struc-
ture of bulk heterojunction conjugated polymer/fullerene solar cells by scanning probe
microscopy. Journal of Physical Chemistry C, 112(18):7241–7249, 2008.

[32] O. G. Reid, K. Munechika, and D. S. Ginger. Space charge limited current measure-
ments on conjugated polymer ﬁlms using conductive atomic force microscopy. Nano
Letters, 8(6):1602–1609, 2008.

[33] G. A. MacDonald, R. A. Veneman, D. Placencia, and N. R. Armstrong. Electrical
property heterogeneity at transparent conductive oxide/organic semiconductor inter-
faces: Mapping contact ohmicity using conducting-tip atomic force microscopy. ACS
Nano, 6(11):9623–9636, 2012.

[34] H. Yang, E. Glynos, B. Huang, and P. F. Green. Out-of-plane carrier transport in
conjugated polymer thin ﬁlms: Role of morphology. Journal of Physical Chemistry C,
117(19):9590–9597, 2013.

178

[35] C. Tanase, E. J. Meijer, P. W. M. Blom, and D. M. de Leeuw. Uniﬁcation of the
hole transport in polymeric ﬁeld-eﬀect transistors and light-emitting diodes. Physical
Review Letters, 91(21):216601, 2003.

[36] X. N. Yang, J. Loos, S. C. Veenstra, W. J. H. Verhees, M. M. Wienk, J. M. Kroon,
M. A. J. Michels, and R. A. J. Janssen. Nanoscale morphology of high-performance
polymer solar cells. Nano Letters, 5(4):579–583, 2005.

[37] G. Li, V. Shrotriya, Y. Yao, and Y. Yang. Investigation of annealing eﬀects and ﬁlm
thickness dependence of polymer solar cells based on poly(3-hexylthiophene). Journal
of Applied Physics, 98(4):043704, 2005.

[38] W. L. Ma, C. Y. Yang, X. Gong, K. Lee, and A. J. Heeger. Thermally stable, eﬃcient
polymer solar cells with nanoscale control of the interpenetrating network morphology.
Advanced Functional Materials, 15(10):1617–1622, 2005.

[39] P. W. M. Blom, M. J. M. de Jong, and M. G. van Munster. Electric-ﬁeld and temper-
ature dependence of the hole mobility in poly(p-phenylene vinylene). Physical Review
B, 55:R656–R659, 1997.

[40] H.-N. Lin, H.-L. Lin, S.-S. Wang, L.-S. Yu, G.-Y. Perng, S.-A. Chen, and S.-H. Chen.
Nanoscale charge transport in an electroluminescent polymer investigated by conduct-
ing atomic force microscopy. Applied Physics Letters, 81(14):2572–2574, 2002.

[41] C. Ionescu-Zanetti, A. Mechler, S.A. Carter, and R. Lal. Semiconductive polymer
blends: Correlating structure with transport properties at the nanoscale. Advanced
Materials, 16(5):385–389, 2004.

[42] R. Yang, A. Garcia, D. Korystov, A. Mikhailovsky, G. C. Bazan, and T. Q. Nguyen.
Control of interchain contacts, solid-state ﬂuorescence quantum yield, and charge
transport of cationic conjugated polyelectrolytes by choice of anion. Journal of the
American Chemical Society, 128(51):16532–16539, 2006.

[43] M. A. Lampert and P. Mark. Current injection in solids. Academic Press, New York,

1970.

[44] A. J. Campbell, D. D. C. Bradley, and D. G. Lidzey. Space-charge limited conduc-
tion with traps in poly(phenylene vinylene) light emitting diodes. Journal of Applied
Physics, 82(12):6326–6342, 1997.

[45] J. Dacuna and A. Salleo. Modeling space-charge-limited currents in organic semicon-

ductors: Extracting trap density and mobility. Physical Review B, 84(19):9, 2011.

179

[46] O. Douhret, A. Swinnen, S. Bertho, I. Haeldermans, J. D’Haen, M. D’Olieslaeger,
D. Vanderzande, and J. V. Manca. Highresolution morphological and electrical char-
acterisation of organic bulk heterojunction solar cells by scanning probe microscopy.
Progress in Photovoltaics: Research and Applications, 15(8):713–726, 2007.

[47] M. Pope and C. E. Swenberg. Electronic processes in organic crystals and polymers.

Oxford University Press, New York, 1999.

[48] A. J. Mozer and N. S. Sariciftci. Negative electric ﬁeld dependence of charge car-
rier drift mobility in conjugated, semiconducting polymers. Chemical Physics Letters,
389(46):438–442, 2004.

[49] D. Wood, I. Hancox, T. S. Jones, and N. R. Wilson. Quantitative nanoscale mapping
with temperature dependence of the mechanical and electrical properties of poly(3-
hexythiophene) by conductive atomic force microscopy. The Journal of Physical Chem-
istry C, 119(21):11459–11467, 2015.

[50] Y. Roichman and N. Tessler. Generalized Einstein relation for disordered semiconduc-
tors - Implications for device performance. Applied Physics Letters, 80(11):1948–1950,
2002.

[51] G. A. H. Wetzelaer, L. J. A. Koster, and P. W. M. Blom. Validity of the Einstein
relation in disordered organic semiconductors. Physical Review Letters, 107(6):4, 2011.

[52] L. J. A. Koster, E. C. P. Smits, V. D. Mihailetchi, and P. W. M. Blom. Device model
for the operation of polymer/fullerene bulk heterojunction solar cells. Physical Review
B, 72(8):085205, 2005.

[53] N. Kopteva, E. O’Riordan, and L. Vulkov. Computational methods for boundary and
interior layers. International Journal of Numerical Analysis and Modeling, 7(3):I–II,
2010.

[54] D. L. Scharfetter and H. K. Gummel. Large-signal analysis of a silicon Read diode

oscillator. IEEE Transactions on Electron Devices, 16(1):64–77, 1969.

[55] S. Selberherr. Analysis and Simulation of Semiconductor Devices. Springer Vienna,

1984.

[56] F. Brezzi, L. D. Marini, S. Micheletti, P. Pietra, R. Sacco, and S. Wang. Discretiza-
tion of semiconductor device problems (i). Special Volume: Numerical Methods in
Electromagnetics, Vol 13, 13:317–441, 2005.

180

[57] E. A. B. Cole. Mathematical and Numerical Modelling of Heterostructure Semiconduc-

tor Devices: From Theory to Programming. Springer-Verlag London, 2009.

[58] J. van Dijk, G. M. W. Kroesen, and A. Bogaerts. Plasma modelling and numerical

simulation. Journal of Physics D: Applied Physics, 42(19):190301, 2009.

[59] C. Jacoboni. Theory of Electron Transport in Semiconductors: Pathway from Elemen-
tary Physics to Nonequilibrium Green Functions, volume 165. Springer-Verlag Berlin
Heidelberg, 2010.

[60] T. Kirchartz and J. Nelson. Device modelling of organic bulk heterojunction solar

cells. Multiscale Modelling of Organic and Hybrid Photovoltaics, 352:279–324, 2014.

[61] . Bjrck. Numerical Methods in Matrix Computations. Springer, Cham, 2015.

[62] P. J. Olver and C. Shakiban. Applied Linear Algebra. Springer, Cham, 2018.

[63] R. A. Usmani. Inversion of Jacobi’s tridiagonal matrix. Computers & Mathematics

with Applications, 27(8):59–66, 1994.

[64] T. Yamamoto. Inversion formulas for tridiagonal matrices with applications to bound-
ary value problems. Numerical Functional Analysis and Optimization, 22(3-4):357–385,
2001.

[65] H. Sirringhaus, P. J. Brown, R. H. Friend, M. M. Nielsen, K. Bechgaard, B. M. W.
Langeveld-Voss, A. J. H. Spiering, R. A. J. Janssen, E. W. Meijer, P. Herwig, and
D. M. de Leeuw. Two-dimensional charge transport in self-organized, high-mobility
conjugated polymers. Nature, 401(6754):685–688, 1999.

[66] J. A. Merlo and C. D. Frisbie. Field eﬀect transport and trapping in regioregular poly-

thiophene nanoﬁbers. Journal of Physical Chemistry B, 108(50):19169–19179, 2004.

[67] W. D. Oosterbaan, J. C. Bolsee, A. Gadisa, V. Vrindts, S. Bertho, J. D’Haen, T. J.
Cleij, L. Lutsen, C. R. McNeill, L. Thomsen, J. V. Manca, and D. Vanderzande.
Alkyl-chain-length-independent hole mobility via morphological control with poly(3-
alkylthiophene) nanoﬁbers. Advanced Functional Materials, 20(5):792–802, 2010.

[68] H. A. van der Vorst. Bi-CGSTAB: A fast and smoothly converging variant of Bi-
CG for the solution of nonsymmetric linear systems. SIAM Journal on Scientiﬁc and
Statistical Computing, 13(2):631–644, 1992.

181

[69] A. A. Grinberg and S. Luryi. Space-charge-limited current and capacitance in double-

junction diodes. Journal of Applied Physics, 61(3):1181–1189, 1987.

[70] H. T. Nicolai, M. M. Mandoc, and P. W. M. Blom. Electron traps in semiconducting
polymers: Exponential versus gaussian trap distribution. Physical Review B, 83(19):5,
2011.

[71] H. T. Nicolai, M. Kuik, G. A. H. Wetzelaer, B. de Boer, C. Campbell, C. Risko,
J. L. Bredas, and P. W. M. Blom. Uniﬁcation of trap-limited electron transport in
semiconducting polymers. Nature Materials, 11(10):882–887, 2012.

[72] A. Salleo, R. J. Kline, D. M. DeLongchamp, and M. L. Chabinyc. Microstructural
characterization and charge transport in thin ﬁlms of conjugated polymers. Advanced
Materials, 22(34):3812–3838, 2010.

[73] F. Liu, Y. Gu, X. B. Shen, S. Ferdous, H. W. Wang, and T. P. Russell. Characterization
of the morphology of solution-processed bulk heterojunction organic photovoltaics.
Progress in Polymer Science, 38(12):1990–2052, 2013.

[74] D. A. Chen, A. Nakahara, D. G. Wei, D. Nordlund, and T. P. Russell. P3HT/PCBM
bulk heterojunction organic photovoltaics: Correlating eﬃciency and morphology.
Nano Letters, 11(2):561–567, 2011.

[75] B. A. Collins, J. R. Tumbleston, and H. Ade. Miscibility, crystallinity, and phase
development in P3HT/PCBM solar cells: Toward an enlightened understanding of
device morphology and stability. Journal of Physical Chemistry Letters, 2(24):3135–
3145, 2011.

[76] C. R. McNeill. Morphology of all-polymer solar cells. Energy & Environmental Science,

5(2):5653–5667, 2012.

[77] D. M. DeLongchamp, R. J. Kline, and A. Herzing. Nanoscale structure measurements
for polymer-fullerene photovoltaics. Energy & Environmental Science, 5(3):5980–5993,
2012.

[78] Y. Kim, S. Cook, S. M. Tuladhar, S. A. Choulis, J. Nelson, J. R. Durrant, D. D. C.
Bradley, M. Giles, I. McCulloch, C. S. Ha, and M. Ree. A strong regioregularity
eﬀect in self-organizing conjugated polymer ﬁlms and high-eﬃciency polythiophene:
fullerene solar cells. Nature Materials, 5(3):197–203, 2006.

182

[79] K. Sivula, C. K. Luscombe, B. C. Thompson, and J. M. J. Frchet. Enhancing the
thermal stability of polythiophene:fullerene solar cells by decreasing eﬀective polymer
regioregularity. Journal of the American Chemical Society, 128(43):13988–13989, 2006.

[80] Z. Y. Wu, A. Petzold, T. Henze, T. Thurn-Albrecht, R. H. Lohwasser, M. Sommer, and
M. Thelakkat. Temperature and molecular weight dependent hierarchical equilibrium
structures in semiconducting poly(3-hexylthiophene). Macromolecules, 43(10):4646–
4653, 2010.

[81] S. Joshi, S. Grigorian, and U. Pietsch. X-ray structural and crystallinity studies of
low and high molecular weight poly(3-hexylthiophene). Physica Status Solidi (a) -
Applications and Materials Science, 205(3):488–496, 2008.

[82] A. M. Ballantyne, L. Chen, J. Dane, T. Hammant, F. M. Braun, M. Heeney, W. Duﬀy,
I. McCulloch, D. D. C. Bradley, and J. Nelson. The eﬀect of poly(3-hexylthiophene)
molecular weight on charge transport and the performance of polymer:fullerene solar
cells. Advanced Functional Materials, 18(16):2373–2380, 2008.

[83] A. Zen, M. Saphiannikova, D. Neher, J. Grenzer, S. Grigorian, U. Pietsch,
U. Asawapirom, S. Janietz, U. Scherf, I. Lieberwirth, and G. Wegner. Eﬀect of
molecular weight on the structure and crystallinity of poly(3-hexylthiophene). Macro-
molecules, 39(6):2162–2171, 2006.

[84] J. A. Lim, F. Liu, S. Ferdous, M. Muthukumar, and A. L. Briseno. Polymer semicon-

ductor crystals. Materials Today, 13(5):14–24, 2010.

[85] H. H. Yang, S. W. LeFevre, C. Y. Ryu, and Z. N. Bao. Solubility-driven thin ﬁlm struc-
tures of regioregular poly(3-hexylthiophene) using volatile solvents. Applied Physics
Letters, 90(17):3, 2007.

[86] D. H. Kim, Y. Jang, Y. D. Park, and K. Cho. Layered molecular ordering of
self-organized poly(3-hexylthiophene) thin ﬁlms on hydrophobized surfaces. Macro-
molecules, 39(17):5843–5847, 2006.

[87] J. Liu, M. Arif, J. Zou, S. I. Khondaker, and L. Zhai. Controlling poly(3-
hexylthiophene) crystal dimension: nanowhiskers and nanoribbons. Macromolecules,
42(24):9390–9393, 2009.

[88] J. Ma, K. Hashimoto, T. Koganezawa, and K. Tajima. End-on orientation of semi-
conducting polymers in thin ﬁlms induced by surface segregation of ﬂuoroalkyl chains.
Journal of the American Chemical Society, 135(26):9644–9647, 2013.

183

[89] I. Osaka and K. Takimiya. Backbone orientation in semiconducting polymers. Polymer,

59:A1–A15, 2015.

[90] Y. K. Lan and C. I. Huang. Charge mobility and transport behavior in the ordered
and disordered states of the regioregular poly(3-hexylthiophene). Journal of Physical
Chemistry B, 113(44):14555–14564, 2009.

[91] J. C. Bolsee, W. D. Oosterbaan, L. Lutsen, D. Vanderzande, and J. Manca. C-AFM
on conjugated polymer nanoﬁbers: Capable of assessing one ﬁber mobility. Organic
Electronics, 12(12):2084–2089, 2011.

[92] F. Zhang, A. Vollmer, J. Zhang, Z. Xu, J. P. Rabe, and N. Koch. Energy level alignment
and morphology of interfaces between molecular and polymeric organic semiconduc-
tors. Organic Electronics, 8(5):606–614, 2007.

[93] Z. Xu, L.-M. Chen, M.-H. Chen, G. Li, and Y. Yang. Energy level alignment of
[6,6]-phenyl C61 butyric acid methyl ester bulk heterojunc-

poly(3-hexylthiophene):
tion. Applied Physics Letters, 95(1):013301, 2009.

[94] B. Huang, E. Glynos, B. Frieberg, H. Yang, and P. F. Green. Eﬀect of thickness-
dependent microstructure on the out-of-plane hole mobility in poly(3-hexylthiophene)
ﬁlms. ACS Applied Materials & Interfaces, 4(10):5204–5210, 2012.

[95] K. Kawano and C. Adachi. Evaluating carrier accumulation in degraded bulk het-
erojunction organic solar cells by a thermally stimulated current technique. Advanced
Functional Materials, 19(24):3934–3940, 2009.

[96] W. Weise, T. Keith, N. von Malm, and H. von Seggern. Trap concentration dependence
of thermally stimulated currents in small molecule organic materials. Physical Review
B, 72(4):235207, 2005.

[97] R. A. Street. Localized state distribution and its eﬀect on recombination in organic

solar cells. Physical Review B, 84(7):075208, 2011.

[98] J. H. Youn, Y. I. Lee, H. T. Moon, M. S. Ryu, J. Kim, and J. Jang. Trap energy
level of P3HT: PCBM-71 bulk heterojunction solar cells with PICTS (photo-induced
current transient spectroscopy). Current Applied Physics, 10(3):S525–S527, 2010.

[99] A. J. Campbell, D. D. C. Bradley, E. Werner, and W. Brtting. Transient capacitance
measurements of the transport and trap states distributions in a conjugated polymer.
Organic Electronics, 1(1):21–26, 2000.

184

[100] T. P. Nguyen. Defects in organic electronic devices. Physica Status Solidi (a) - Appli-

cations and Materials Science, 205(1):162–166, 2008.

[101] C. Renaud and T.-P. Nguyen.

Identiﬁcation of the nature of trapping centers in
polyspirobiﬂuorene based diodes by using electrical characterization. Journal of Ap-
plied Physics, 107(12):124505, 2010.

[102] K. Fukuda, T. Sekitani, and T. Someya. Eﬀects of annealing on electronic and struc-
tural characteristics of pentacene thin-ﬁlm transistors on polyimide gate dielectrics.
Applied Physics Letters, 95(2):023302, 2009.

[103] A. Sharma, P. Kumar, B. Singh, S. R. Chaudhuri, and S. Ghosh. Capacitance-voltage
characteristics of organic schottky diode with and without deep traps. Applied Physics
Letters, 99(2):023301, 2011.

[104] I. H. Campbell, D. L. Smith, and J. P. Ferraris. Electrical-impedance measurements

of polymer light-emitting-diodes. Applied Physics Letters, 66(22):3030–3032, 1995.

[105] G. Garcia-Belmonte, P. P. Boix, J. Bisquert, M. Sessolo, and H. J. Bolink. Simultaneous
determination of carrier lifetime and electron density-of-states in P3HT:PCBM organic
solar cells under illumination by impedance spectroscopy. Solar Energy Materials and
Solar Cells, 94(2):366–375, 2010.

[106] J. C. Scott and G. G. Malliaras. Charge injection and recombination at the metalor-

ganic interface. Chemical Physics Letters, 299(2):115–119, 1999.

[107] P. N. Murgatroyd. Theory of space-charge-limited current enhanced by frenkel eﬀect.

Journal of Physics D: Applied Physics, 3(2):151, 1970.

[108] W. F. Pasveer, J. Cottaar, C. Tanase, R. Coehoorn, P. A. Bobbert, P. W. M. Blom,
D. M. de Leeuw, and M. A. J. Michels. Uniﬁed description of charge-carrier mobilities
in disordered semiconducting polymers. Physical Review Letters, 94:206601, 2005.

185