To Marilyn, the strongest person I know.

iii

ACKNOWLEDGMENTS

This thesis is only possible thanks to the combined effort of a great number of people whom,
over the years, have helped guide and support me. I would specifically like to acknowledge
the support of Megan, Ania, Carol, Jiwu, Jenni, Giuseppe, Gift, Connor, Shannon, Jim,
Chelsea, Paul, the Jons, and Josh. Thanks to my family for supporting me through thick
and thin. Thanks to Marilyn for ceaselessly believing in me. And of course thanks to my
adviser, Phil Duxbury, who’s really to blame for all of this.

iv

TABLE OF CONTENTS

List of Tables . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . viii

List of Figures . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

ix

Chapter 1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
1.1 Motivation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
1.2 Outline . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

1
1
5

Chapter 2 The Basics of OPV . . . .
2.1 Solar Cell Devices . . . . . . . . . .
2.2 Organic Photovoltaics . . . . . . .
2.3 P3HT/PCBM OPV Device Design

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

8
8
17
21

Chapter 3 Monte Carlo Device Model for OPV . . . . . . . . . .
3.1 The Scales of Modeling and Literature Review . . . . . . . . . . .
3.2 Basics of the Dynamic Monte Carlo OPV Model . . . . . . . . . .
3.2.1 First Reaction Method . . . . . . . . . . . . . . . . . . . .
3.2.2 Exciton Events . . . . . . . . . . . . . . . . . . . . . . . .
3.2.2.1 DMC exciton generation . . . . . . . . . . . . . .
3.2.2.2 DMC exciton hop rates . . . . . . . . . . . . . .
3.2.2.3 DMC exciton dissociation rate . . . . . . . . . .
3.2.3 Charge Events . . . . . . . . . . . . . . . . . . . . . . . . .
3.2.3.1 DMC charge hop rates . . . . . . . . . . . . . . .
3.2.3.2 DMC charge collection and dark charge injection
3.2.3.3 DMC charge recombination rate . . . . . . . . . .
3.3 Marcus Hop Rates . . . . . . . . . . . . . . . . . . . . . . . . . .
3.4 Simple Bilayer Morphology Example . . . . . . . . . . . . . . . .
3.4.1 The effects of energetic disorder on a bilayer system . . . .
3.4.2 The effects of sample depth on a bilayer system . . . . . .
3.5 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.

25
25
28
29
30
30
32
33
34
37
38
39
40
42
43
46
47

Chapter 4 A Novel Method to Simulate
4.1 Principles of the Fourier Transform . .
4.2 Basics of scattering . . . . . . . . . . .
4.2.1 Single Scattering Event . . . . .
4.2.2 Multiple Atom Scattering . . .
4.2.3 Scattering strength . . . . . . .

.
.
.
.
.
.

.
.
.
.
.
.

.
.
.
.
.
.

.
.
.
.
.
.

.
.
.
.
.
.

.
.
.
.
.
.

50
51
54
57
58
60

v

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

Small
. . . .
. . . .
. . . .
. . . .
. . . .

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

Angle
. . . .
. . . .
. . . .
. . . .
. . . .

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

Scattering
. . . . . . .
. . . . . . .
. . . . . . .
. . . . . . .
. . . . . . .

LIST OF TABLES

Table 3.1

Parameters used in the bilayer example, except when otherwise specified. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

45

Table 6.1

Parameters extracted from SANS based on polydisperse sphere model 133

Table 6.2

SLD values for the four different fractions of dispersed PCBM morphologies seen in figure 5.17, used in simulated scattering in figure
5.17. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 143

viii

LIST OF FIGURES

Figure 2.1

Figure 2.2

Figure 2.3

Figure 2.4

Figure 2.5

Figure 2.6

Schematic of a simple PN junction in a solar cell at short-circuit
condition, as evidenced by the energy levels of the anode and cathode being equalized. Absorbed photon excites an electron from the
valence band into the conduction band. The electron and hole can
then follow the chemical potential gradient to the cathode and anode
respectively. The Fermi energy, Ef , exists between the valence and
conduction bands. For interpretation of the references to color in this
and all other figures, the reader is referred to the electronic version
of this dissertation. . . . . . . . . . . . . . . . . . . . . . . . . . . .

9

Demonstration of applied voltage to a simple solar cell band diagram.
As the open circuit voltage, Voc is passed, the runaway diode behavior
emerges. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

12

Equivalent circuit for a solar cell, showing photon absorption generating a short-circuit current density, JSC , which has losses represented
as a parallel and series resistances, RP and RS respectively. Dark
current, Jdark , is represented as a diode in the reverse direction of
the photogenerated current. . . . . . . . . . . . . . . . . . . . . . . .

14

Example of a current density (J) versus applied voltage (Vapplied )
graph, demonstrating the diode like behavior of a solar cell. Curves
are shown for both a device in the dark (Jdark ) and under illumination (Jlight ). Also shown on the graph are the open circuit voltage,
VOC , short-circuit current, JSC , thick dashed line to better understand the fill factor, F F . Notable behavior points occur at short
circuit, a , near maximum power, b , near open circuit, c , and finally at runaway diode-like behavior, d . . . . . . . . . . . . . . . .

15

Example of a molecular PV device, showing photoexcitation occurring in the donor material, then charge separation and transport to
the contacts. Note the lack of any band bending, but rather hopping
between HOMO levels (hole transport) and LUMO levels (electron
transport). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

19

Qualitative illustration of how a charge hops on a length scale related
to thermal vibration of a local energy landscape. . . . . . . . . . . .

20

ix

Figure 2.7

Figure 3.1

Figure 3.2

Figure 3.3

Figure 4.1

Schematic of a P3HT/PCBM based OPV device, with the physical
layers shown in the top figure and the corresponding energy levels
on the bottom. In this diagram, light passes in from the right hand
side and is absorbed in the active layer. Electrons travel to the back
metal contact (typically aluminum) and holes, passing through the
PEDOT:PSS, are collected at the top contact (ITO). . . . . . . . .

24

Parabola geometry used to derive Marcus rate activation energy. The
net energy change from the reaction is the Gibb’s free energy of reaction, ∆Go . λ is energy of reorganization. Gibb’s free energy of
activation is found by solving for the crossover point of the parabolas. Hij is the interaction energy at the intersection of the two states.
In the two states shown, the state on the right has a lower overall
configurational energy by ∆Go , but the activation barrier must be
overcome for charges to transport to this lower configuration. . . .

41

Example of how energetic roughness would effect a bilayer. Data
presented as collection efficiency versus applied voltage (Vapplied ),
demonstrating the well known diode-like behavior. Note that no
dark current is being simulated, so no runaway current in the opposite direction is possible. Figure 3.2a is the charge collection efficiency, χc-collect , while figure 3.2b is the exciton dissociation efficiency, χx-diss . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

44

DMC calculations showing the effect of varying total bilayer thickness
on charge collection efficiency, χc-collect (top), and exciton dissociation efficiency, χx-diss (bottom). Presented are the results for systems
of thickness L = 20, 40, 60, 80, and 200, with total morphology size
35 × 35 × L, and with periodic boundary conditions in the x and y
dimensions. No dark charge is simulated here. All other parameters
used in the simulation are shown in table 3.1. . . . . . . . . . . . . .

48

A demonstration of the Fourier relationship. Figure 4.1a shows an
arbitrary function composed of superimposed Gaussian distributions
of varying heights. Figure 4.1b demonstrates the real component
of the Fourier transform of the function in figure 4.1a. Figure 4.1c
demonstrates the auto-correlation of the function in figure 4.1a. . . .

53

x

Figure 4.2

Convolutions of Gaussian functions of varying widths with a series of
delta function peaks. The function describing the delta functions is
f (x) = δ(x + 10) + 2δ(x + 5) + 3δ(x) + 2δ(x − 5) + δ(x − 10). The
Gaussian is described as g(x) = exp −x2 /2σ 2 . Widths of Gaussian
distributions used in these three examples are σ 2 = .1 in figure 4.2a,
σ 2 = .5 in figure 4.2b, and σ 2 = 3 in figure 4.2c. . . . . . . . . . . .

54

Figure 4.3

Single scatter event example . . . . . . . . . . . . . . . . . . . . . .

57

Figure 4.4

Comparison of methods used to simulate small angle scattering, highlighting the overall mathematical similarities of the two methods. . .

65

˚
Comparison of scattering from a single R = 100 A sphere calculated
using the analytical form factor and the DFM algorithm. Analytical
form factor provided by the IRENA SAS Macros in Igor, where sphere
to solvent contrast is taken 1-to-0. DFM calculation performed over
a randomly sampled set of 107 pairs, or .05% of the total model pairs.
The lattice spacing here is taken as al = 1 ˚. . . . . . . . . . . . . .
A

74

˚
Simulated scattering from cylinders of length, L = 100 A, and ra˚, calculated using the analytical form factor for a
dius, R = 25, 50 A
cylinder and compared to the scattering calculated with the DFM
algorithm. Analytical form factors provided by IRENA SAS macros
in Igor, where cylinder to solvent contrast is taken as 1-to-0. DFM
calculation performed over a randomly sampled set of 107 pairs, or
.005% of the total system pairs. . . . . . . . . . . . . . . . . . . . .

75

Simulated scattering from an ideal chain polymer model with 106
monomers with resultant radius of gyration, rg ≈ 398 ˚. The DFM
A
calculation is performed by sampling the given number of random
pairs on 1000 random polymer configurations. Note that for the size
of polymer simulated here, the total calculation would require a sum
of nearly 5 × 1011 total pairs per polymer configuration. Even with
only 100 pairs sampled per polymer, the analytical solution is exactly
recovered up to high-q values. . . . . . . . . . . . . . . . . . . . . .

77

Figure 4.5

Figure 4.6

Figure 4.7

xi

Figure 4.8

Figure 4.9

Figure 4.10

Figure 4.11

Figure 5.1

DFM calculation of scattering from a system of 697 spheres of radius
3
r = 20 ˚ placed randomly inside a box of size 8003 ˚ which were seA
A
lected to have a minimum distance between sphere centers of L = 100
˚. Demonstrated are the effect of calculating the scattering on this
A
morphology while points were randomly removed from the structure
using a “fuzzying” routine via a Gaussian probability with width σ.
This procedure reduces the inherent hard-edged box effects seen as
high frequency oscillations at the expense of total signal contrast and
sample size. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

79

A demonstration of how simulated scattering can be used to study the
effects of annealing a complex, interpenetrating structure such as the
Ising model shown here. Also demonstrated is the use of applying a
smearing function to somewhat under-sampled data. The raw data is
presented as dashed lines and the σq = 10% smeared data is presented
as solid lines. Inset figures show cross-sections of these two-phase
morphologies as they are annealed from a = 4 to a = 16. . . . . . . .

83

˚
Spheres of constant radii, r = 10 A, randomly placed inside an L =
˚ box with a minimum center-to-center distance of D = 20, 40,
400 A
and 100 ˚. For comparison, the form factor for a single sphere of
A
radius r = 10 ˚ is also shown. The DFM calculation is performed by
A
9 pairs, and smearing resultants by 5%. . . . . . . . . .
sampling 10

86

Scattering from a series of 370 spheres randomly placed in a L =
800 ˚ box at a minimum center-to-center distance of D = 100 ˚.
A
A
The sphere radii are selected from a Schulz distribution shown in
the inset to the figure. For each distribution, the average radius is
r = 20 ˚, but the deviation used in each case is σ = 5, 10, and
A
15, which generated resultant polydispersity of P DI = 1.36, 3.27,
and 8.15 respectively. For comparison, also shown is the scattering
from spheres placed under similar conditions with constant radius,
r = 20 ˚, corresponding to P DI = 1.0. DFM calculation performed
A
by sampling 109 pairs, and smearing results 5%. The inset shows the
Schulz polydisperse sphere distributions used to select radii. . . . . .

88

Figure 5.1a demonstrates the work of Steve Forrest’s group using a
simple Ising model to roughen the interface between a set of bilayers[1].
Figure 5.1b demonstrates the work of Alison Walker’s group using
an Ising model to study the effect of BHJ feature size on device
performance[2]. . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

91

xii

Figure 5.2

Schematic of reflectometry measurement, with incident wavelength,
λi reflecting at an angle of θi off of the surface, with a transmitted wave having associated wavelength λt at angle θt . Note that
this figure demonstrates the reflection and transmission angles associated with only the first layer of a multilayer system. Further
reflections/transmissions are implied but not shown. . . . . . . . . .

95

Figure 5.3

Two qualitative 1-d examples of an incident wave, ψi , which originates in air (βair = 0) interacting with a sample. Figure 5.3a demonstrates a incident wave ψi incident on a sample with index of refraction n, with subsequent transmitted wave ψt and reflected wave ψr .
Figure 5.3b demonstrates the incident wave on a multilayered sample, with N different layers atop a substrate layer S. Each layer has
a propagating component, ψ − to the left and a reflected component
ψ + to the right. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 100

Figure 5.4

A demonstration of the simulated reflectometry data from a simple
step function, with various amounts of “roughing” of the edges to
smooth the transition between layers. . . . . . . . . . . . . . . . . . 107

Figure 5.5

An example of various step functions with combinations of sharp
and smooth transitions between steps, demonstrating how identical
reflectometry data can be generated from structures with different
SLD density profiles. . . . . . . . . . . . . . . . . . . . . . . . . . . 107

Figure 5.6

An example of simulated reflectometry data from layers of alternating SLD values, in order to demonstrate how SLD contrast effects
reflectometry data. . . . . . . . . . . . . . . . . . . . . . . . . . . . 108

Figure 5.7

A demonstration of the effect of varying the background layer of an
otherwise fixed SLD density profile. . . . . . . . . . . . . . . . . . . 108

xiii

Figure 5.8

Various images and models of polymer/fullerene BHJ morphologies.
Figure 5.8a shows bright field TEM images of an polymer/fullerene
BHJ OPV film composed of PTB7:PC61BM/DCB + DIO, with artificial black lines added to highlight differences between the polymerrich and nanoparticle-rich domains[3]. Figure 5.8b is a cartoon illustrating the nanomorphology proposed in BHJ OPV devices under the
“Rivers and Streams“ model[4]. Figure 5.8c is an image of the x-y
plane in a polymer(PF10TBT)/nanoparticle(PCBM) film generated
from electron tomography data, with the polymer regions highlighted
in light and PCBM regions in dark[5]. Figure 5.8d shows a sample
cross-section of reconstructed data taken on a P3HT/ZnO OPV device obtained through electron tomography with ZnO being yellow
and P3HT being transparent[6]. Figures 5.8e and 5.8f show sulfur
maps and carbon maps, respectively, taken through energy filtered
electron tomography (EFTEM) for 1:1 weighted P3HT/PCBM films
which were annealed for 30 minutes at zero defocus and with dimensions approximately 400 nm across. Light regions in the sulfur map
correspond to higher concentrations of P3HT-domains[7]. . . . . . . 112

Figure 5.9

Demonstration of Ising annealement over time on an example morphology annealed with parameters J = 2.0, kB T = 2.0, Jdens = 15.0,
and Jpop = 4.1, fitting the annealed reflectivity profile shown in figure
5.10. As the morphology was being annealed, exciton transport calculations were performed for each integer value of feature size reached
(a = 1, 2, 3, . . .) using a simplified form of the DMC model as is described in chapter 6. The results of these calculations are presented
in the inset as exciton dissociation efficiency versus characteristic feature size (a). The calculation was performed using decay of lengths,
Lx = [.01, .1, 1, 10, 100], which are shown. Four specific morphologies
are highlighted, with their corresponding results also highlighted in
the inset. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 113

Figure 5.10

PCBM density profile fits derived from neutron reflectivity data and
an associated cartoon theorizing the morphology[8]. . . . . . . . . . 114

Figure 5.11

Model morphologies (L = 400 × 400 × 400) fit to the unannealed
(as-cast) density profile in figure 5.10. Cross sections of the x-y plane
are shown for 4 different feature sizes, a = 4 (5.11a), a = 8 (5.11b),
a = 12 (5.11c), and a = 16 (5.11d). Brown represents PCBM sites
and white represents P3HT sites. . . . . . . . . . . . . . . . . . . . . 122

xiv

Figure 5.12

Model morphologies fit to the annealed density profile seen in figure
5.10. Cross sections of the x-y plane are shown for 4 different feature
sizes(figure), a = 4 (5.12a), a = 8 (5.12b), a = 12 (5.12c), and a = 16
(5.12d). Brown represents PCBM sites and white represents P3HT
sites. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 123

Figure 5.13

A comparison of electron tomography[5] performed on PPV:PCBM
BHJ solar cells to simple morphologies generated using the constrained
Ising model method. First, figure 5.13a demonstrates a zoomed in
400 nm × 400 nm section of electron tomography taken of a BHJ
morphology[5]. Next, in figure 5.13b, is an artificially colored twophase view of this image. Finally, figure 5.13c is a sample side-cut
of an example morphology generated with our model with comparable scale, containing feature size a = 8.1. The interpenetrating
morphology of figures 5.8b and 5.8c are evident. . . . . . . . . . . . 124

Figure 5.14

3-d views of the as-cast and annealed model morphologies, at feature
size a = 8, presented such that the abrupt change in PCBM density
seen near the contacts which dramatically alters the pathways to the
contacts is visible. Brown represents PCBM sites and white represent
P3HT sites. Note that as each of these cases has the same feature
size, the morphological differences are consequences of the utilized
density profiles. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 125

Figure 5.15

Assumed aggregate PCBM density profile from the three-phase model.126

Figure 5.16

2-d cross-sections demonstrating the sequestered PCBM morphologies, with figures 5.16a, 5.16b, 5.16c, and 5.16d representing x =
0, x = .1, x = .15, and x = .2 reduced volume fractions respectively.
Brown sites represent aggregate PCBM, white represents the mixed
PCBM/P3HT phase. . . . . . . . . . . . . . . . . . . . . . . . . . . 127

Figure 5.17

Superpositions of 15 cross-sections of the reduced aggregate PCBM
morphologies, where only aggregate PCBM is displayed as green for
sites connecting to the electron collecting electrode and white for
sites which are disconnected. The four cross sections shown are for
the x = 0, x = .1, x = .15, and x = 2. dispersed PCBM cases. . . . . 128

Figure 5.18

Figure 5.18a shows the P3HT density profiles used in model generation here, demonstrating both the crystallized fraction (c-P3HT) and
the total fraction. Figure 5.18b shows a top down view of c-P3HT
fiber network. Superposition of 20 total layers shown, with each fiber
being 4 layers deep. . . . . . . . . . . . . . . . . . . . . . . . . . . . 129

xv

Figure 5.19

Views of the unannealed c-P3HT fiber-network based morphologies.
Figure 5.19a shows a 2-d cross section of the initial configuration of a
model morphology, which has had a c-P3HT fiber network inlaid but
has not yet undergone annealing, such that the PCBM (brown) and
amorphous P3HT (white) are still well mixed. Figure 5.19b is a 3-d
view of the same morphology. . . . . . . . . . . . . . . . . . . . . . 129

Figure 5.20

Views of the annealed (a=16) c-P3HT fiber-network based morphologies. Figure 5.20a shows a 2-d cross section of the final configuration
of a model morphology, which has undergone annealing of the amorphous P3HT and PCBM (but left the c-P3HT network frozen) to
characteristic feature size, a=16, such that the phases are now well
segregated. Figure 5.20b is a 3-d view of the same morphology. . . . 130

Figure 6.1

Demonstration of how a lattice constant, blattice , can be extracted
from the Bragg-like nearest neighbor peak generated by the DFM
method. Once final scaling is performed such that the calculated
scattering profile is best it to experimental data, these lattice peaks
(as highlighted here) will correspond to 2π/blattice . Shown are the
fits for the annealed and as cast two-phase morphology seen in figure
6.5. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 134

Figure 6.2

Small angle neutron scattering (SANS) data of the P3HT/PCBM
BHJ thin film mimics, from the thesis of Jon Kiel[9]. Figure 6.2a
shows the original SANS data for the as cast and annealed systems,
as well as the polydisperse sphere fits. Figure 6.2b is a cartoon by
Jon Kiel, suggesting how the agglomerate PCBM size increases might
be interpreted. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 151

Figure 6.3

Demonstration of Schulz Distributions of the polydisperse fits for the
as cast experimental SANS data(R = 1.0, σ = .75) and the annealed
experimental SANS data (R = 5.9, σ = 2.4). . . . . . . . . . . . . . 152

Figure 6.4

Demonstration of how the plateau shifts to greater intensity and the
bend shifts to lower q as the constrained Ising model anneals a BHJ
morphology to greater values of a. Dashed lines are the unsmeared
data and solid lines have been smeared 10%. The inset presents the
same data (smeared only), zoomed out so as to better observe the
plateau. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 152

xvi

Figure 6.5

Simulated SANS profiles for the as cast and annealed morphologies,
showing fits to the experimental SANS data, as well as previous fits
performed with a polydisperse sphere model. Inset in the figure are
3-d views of the as cast and annealed morphologies. . . . . . . . . . 153

Figure 6.6

Density profiles of the as cast and annealed morphologies, with the
simulated morphologies compared to the experimental reflectometry
fits. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 154

Figure 6.7

Cross sections of the model morphologies, generated with experimentally measured density profiles, and those generated to an identical
characteristic feature size, a, but with a flat target density, which we
refer to as Ising-only. The left hand column has feature size a = 4,
while the right hand column has feature size a = 16. The top row are
morphologies generated with the corresponding reflectivity informed
density profile, while the bottom row is the Ising-only models, with
a flat target density profile. . . . . . . . . . . . . . . . . . . . . . . . 155

Figure 6.8

Results of exciton behavior comparing realistic morphologies to Ising
only morphologies. Plot is of exciton dissociation efficiency (χdiss )
vs. exciton decay length (Lx ). Maroon crosses represent the as cast
morphology, black circles the annealed morphology, green diamonds
the Ising only of feature size as cast (a=4) and yellow boxes the Ising
only of feature size annealed (a=16). . . . . . . . . . . . . . . . . . . 156

Figure 6.9

Local feature sizes for the as cast, annealed, and Ising only (a=4,16)
morphologies. Local feature size is calculated with the familiar form
a = 3V /S, but V is taken as the number of sites in a single slab
V = L × L and S is taken as the total number of BHJ interfaces seen
by the sites on a specific layer. . . . . . . . . . . . . . . . . . . . . . 157

Figure 6.10

Charge collection efficiency comparison between the as cast, annealed,
and Ising only (a=4,16) morphologies. Presented as the ratio of total
charges collected at their respective contacts versus charge recombination length. For the realistic morphologies, solid lines represent
electron collection and dashed lines represent holes. For the Ising
only morphologies, there was no difference between electron and hole
collection efficiencies, so they are presented as a single solid line. . . 158

Figure 6.11

Simulated SANS scattering for the four cases of x = 0, x = .1, x = .15,
and x = 2. fractional dispersed PCBM. The inset demonstrates the
zoomed out view of the same simulated scattering. . . . . . . . . . . 159

xvii

Figure 6.12

Exciton dissociation efficiency, χdiss , versus exciton decay length, Lx ,
for the four cases of figure 5.17. Inset is a magnified portion of the
data. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 160

Figure 6.13

Electron collection efficiency, χq , versus free charge recombination
length, Lq , for the four cases of figure 5.17. Inset is a zoomed out
version of the data, on a log-log graph. . . . . . . . . . . . . . . . . 161

Figure 6.14

Simulated SANS comparisons of the simple two-phase model (black)
and the crystallized P3HT model (dashed blue), as well as the experimental data (yellow boxes) and polydisperse sphere fits (black dots).
Shown inset are views of the c-P3HT crystal network on the left, and
a 3-d view of this morphology where PCBM has been colored brown,
c-P3HT is white, and amorphous P3HT is gray. . . . . . . . . . . . 162

Figure 6.15

Exciton dissociation efficiency for the as cast (maroon triangles), annealed (black circles), and c-P3HT (light blue diamonds) morphologies. Exciton dissociation efficiency, χx-diss , is defined as the ratio of
excitons successfully diffusing to and dissociating at a heterojunction
interface over total excitons generated. Exciton decay length, Lx , is
defined as the lattice constant times the ratio of the exciton hop rate
to the exciton decay rate. . . . . . . . . . . . . . . . . . . . . . . . . 163

Figure 6.16

Charge collection efficiency for the as cast (maroon up-pointing triangles), annealed (black circles), c-P3HT (light blue diamond) models.
Solid lines represent electron collection efficiency, while dashed lines
represent hole collection. The enhanced mobility hole transport is
also shown (dark blue right-pointing triangles). Charge collection efficiency, χq , is defined as the ratio of charges successfully collected at
their respective contacts over total charges generated. Data presented
as a function of charge recombination length, Lq . . . . . . . . . . . . 164

Figure A.1

Geometry of a single charge, e, adjacent to an infinite plane separating
materials of different relative dielectrics, 1 and 2 . . . . . . . . . . . 169

Figure A.2

Geometry of a pair of charges, with the potential on charge q given
by equation A.12 . . . . . . . . . . . . . . . . . . . . . . . . . . . . 173

Figure A.3

Single slit example . . . . . . . . . . . . . . . . . . . . . . . . . . . . 174

Figure A.4

Single slit example data . . . . . . . . . . . . . . . . . . . . . . . . . 175

Figure A.5

Double slit example . . . . . . . . . . . . . . . . . . . . . . . . . . . 179
xviii

Figure A.6

A simple random walker simulation with hops either confined to a
cubic lattice or allowed in all directions but with fixed hop distance.
Averaged over 107 walkers. On the LHS is the average modulus of
the radius traveled vs. time. On the RHS is the average variance
over time vs time, which according to equ.(A.65), will give 6D. . . . 193

Figure A.7

A demonstration of how desired hop rate is affected by attempt rate
in the FRM method. . . . . . . . . . . . . . . . . . . . . . . . . . . 194

xix

Chapter 1
Introduction

1.1

Motivation

Human beings have developed a rather addictive habit in the past century: the somewhat
thoughtless consumption of great amounts of energy for our comfort and enjoyment. Unfortunately, our capability to safely and efficiently generate energy has not grown nearly as
quick as our ability to consume it. As our population continues to grow (another habit we as
a species are rather fond of) it becomes clear that if we wish to maintain our modern lifestyle,
we must explore and develop sources of alternative energy. I have personally always thought
the sun seemed rather full of energy, and why don’t we try to harness some of that. Most
current forms of easily accessible energy can be traced back to the sun’s influence[10]. Coal
and oil, perhaps two of our worst habits, began their stories as ancient plants and animals
who millions of years ago found themselves sustained, directly or indirectly, by the energy of
the sun. Today we go to great lengths to discover the high energy chemical bonds of these
previous tenants, often found deep underground. After much industry and effort, we extract
these complex carbon chains hidden below, and proceed to combust them for our comfort
and entertainment or simply to pour them into our gas tanks. Wind and water both move
turbines for us, but the currents and motions utilized therein are consequences of the complex thermal and weather cycles in which our sun plays the rather central role of “heater”.
It seems justified, given it’s extensive past experience and outstanding track record, that
1

we directly consider the sun as principle energy provider for our species. I mean, it’s just
going to keep shining down on us for the foreseeable future anyway, so further development
of photovoltaic, photothermal, and solar to fuel technologies seems fitting.
Arguably, the greatest hurdle one faces in order to fully realize the concept of a solar
powered energy economy lies in the cost and efficiency levels of present day photovoltaic
(PV) technologies[11, 12]. One promising avenue of research in the last decade has been thin
film organic photovoltaic (OPV) devices. With devices typically composed of polymer and
nanoparticles arranged creatively on the length scales of 100-300 nm thick, OPV devices have
the potential to be both lighter and cheaper than conventional PV, not necessarily replacing
them, but becoming a supplementary energy source in designs where conventional PV could
not be used[13]. One of the remaining hurdles towards development of industrial quality
OPV devices is that their current maximum efficiencies are much lower than conventional
PV (about 10% compared to 40%)[14, 15]. It is a principle focus of this thesis to develop
tools that help characterize and explain the emergence of the performance limiting behaviors
in OPV devices.
Unlike their inorganic counterparts, absorption of a photon in OPV does not immediately
generate a free charge, but instead a tightly bound exciton state. In order for this exciton
to be dissociated into free charges an additional energetic asymmetry is required inside the
active layer of the device, such as the difference in energy levels between the acceptor and
donor material. A complication arises due to the relatively short lifetimes of excitons inside
these materials (≈ 10 ps), such that even though the exciton may be mobile throughout the
device it is only capable of successfully separating into free charges if it is generated within
close proximity of an acceptor/donor interface. The parameter known as the exciton decay
length (Lx-decay ), which corresponds to the approximate distance an exciton will diffuse
2

prior to it’s radiative decay, is useful in describing system morphologies. This requirement
has led to widespread use of the so called bulk-heterojunction (BHJ) design, in which acceptor and donor material are finely mixed on length scales of order Lx-decay such that most
photogenerated excitons easily discover an interface prior to decay. This ease of interface discovery increases the efficiency at which absorbed photons are converted into usable current,
known as the quantum efficiency (QE) of the device, by minimizing exciton decay losses.
An unfortunate downside to this design comes from the fact that once excitons dissociate
into free charges (electrons in the acceptor material and holes in the donor material) the
conductive transport path from the geminate location (site of exciton dissociation) to the
charge contacts will inevitably be more tortuous and complex in a BHJ structure compared
to a simple planar heterojunction. Thus, while the BHJ significantly enhances exciton dissociation efficiencies, it may lead to significant free charge loss. Great care must be taken to
design BHJ morphologies which optimize overall quantum efficiency, not simply maximize
exciton dissociation behavior while neglecting the effect on free charge.
This complex relationship between morphology and device performance[16, 17] is in no
way unique to the field of organic photovoltaics. Many next-generation energy applications
have material properties and device performance metrics which are strongly correlated to
their internal nanoscale features or morphologies. Batteries, LEDs, fuel cells, sensors, and
organic photovoltaics all rely on transport mechanisms intrinsically associated with the interpenetrating, percolative nanostructures employed by these devices[18, 19], and yet no models
exist which attempt to directly extract morphology from experimental data, with transport
calculations instead relying upon assumed or simplified models. We present a simple model
which generates three dimensional interpenetrating percolative nanostructural morphologies
that are simultaneously consistent with multiple types of experimental data. Though the
3

work emphasized in this thesis is on the P3HT/PCBM bulk heterojunction morphology,
modeled using small angle neutron scattering (SANS) and neutron reflectometry data, the
method is easily adapted to many other systems and forms of experimental data.
In order to characterize these morphological models, we employ transport calculations
based on modified dynamic Monte Carlo (DMC) device simulations. The DMC model is
excellent at exploring how subtle changes in the nanoscale features of a configurational
morphology will influence exciton and charge behavior in the active layer of the OPV. In
this thesis we focus on the effect of nanoscale morphology of a polymer/nanoparticle BHJ
OPV device. The morphology models used are extracted directly from and then fit to SANS
and neutron reflectometry data. We feel that by using morphological models which are more
directly correlated with available experimental data, any subsequent theories derived from
these models are inherently more accurate than those developed from assumed morphologies.
We will compare how the use of assumed morphologies[2] versus the presented more realistic
morphologies alters the outcome of transport calculations.
In this thesis, we focus on utilizing SANS and neutron reflectometry data as it pertains
to the OPV devices under investigation, but many other viable experimental datasets exist
which could readily be incorporated into future modeling efforts. A wide variety of techniques have been employed to better understand the complex systems involved in these BHJ
morphologies, including AFM studies probing local photocurrents[20] and space-charge limiting behavior[21], STM[22], SEM[23], and electron tomography[6] based studies which have
attempted to directly map the morphology, and attempts to directly observe carrier photogeneration with Kelvin probe force microscopy[24]. Each of these measurements is valuable
in it’s own right, however none so far have attempted to directly extract a morphological
model required for device simulation.
4

While the reflectometry data used in this thesis has previously been fitted to an exact
PCBM density profile[8], and is thus easily incorporated into our morphological model, SANS
data is not as directly invertible. While previous standard methods utilizing polydisperse
sphere fits have estimated aggregate PCBM cluster size[25], these fits have little physical
correlation with the actual system. In order to fit our theoretical models directly to the
SANS data, we have developed a novel algorithm we refer to as the distribution function
method (DFM), which calculates the SANS profile of a theoretical morphology. Although
not as conceptually straightforward as other methods of calculating a theoretical small angle
scattering profile[26], the DFM method allows for a more flexible and scalable calculation of
a simulated scattering profile, allowing us to easily tune the precision of the calculation while
exploring the influence of various features of a morphology on the overall SANS profile.
Utilizing our modified Ising model, novel DFM scattering algorithm, and DMC based
device simulation, we will present a step forward in characterization and understanding of
the effect of nanoscale morphology on overall device performance of P3HT/PCBM based
OPV devices.

1.2

Outline

The focus of this thesis is to describe the developed methods and results from an effort to
better characterize and understand the effect of nanoscale internal morphology on various
aspects of OPV performance. Though much work has gone into this question in the past[2,
27, 1, 28, 29], this thesis will present modifications to the dynamic Monte Carlo model
frequently employed to study OPV, which act to directly refine the morphological model
through the use of neutron reflectometry and small angle scattering data. Both the method

5

used to generate the model morphology itself and the calculations performed to evaluate
the theoretical performance will be explained in detail, and the importance of incorporating
previously underutilized experimental data will be demonstrated. The bulk of new work
presented is to be found in chapters 4 through 6.
Chapter 2 explains the basic principles of photovoltaic devices as well as the differences
between conventional PV and thin film OPV. This chapter will also introduce and motivate
the strong relation between nanoscale morphology and OPV device performance. Finally,
this chapter will present some of the attempts that have been made to enhance device
performance through morphology manipulation. This chapter is principally review.
Chapter 3 describes the dynamic Monte Carlo (DMC) transport model used to characterize the performance of a particular morphology. The details of the model will be presented, along with example calculations on a set of simple morphologies. The modifications
we introduce to the basic DMC model will then be motivated and explained in detail.
Chapter 4 presents a novel method to efficiently generate a simulated small angle scattering profile from any model morphology which we refer to as the distribution function
method (DFM). The benefits and limitations of the DFM method will be discussed, particularly in comparison to the more straightforward direct method employed by other groups[26].
Several examples of well known structures will be presented, demonstrating the validity of
the DFM method. This chapter reviews previous small angle scattering calculation methods,
as well as introducing our new procedure, and the results of this chapter will be published
in a techniques paper.
Chapter 5 describes the modified Ising model we use to generate theoretical P3HT/PCBM
BHJ structures that are consistent with neutron reflectometry and SANS data. Though the
reflectometry data is directly incorporated into the Ising model, SANS data is not as di6

rectly invertible, and thus we demonstrate how through the use of the DFM method we
select and refine our structures until they are consistent with the experimental SANS data.
Additionally, we will demonstrate how to expand the simple two-phase model for use in
more elaborate system descriptions. One such presented case is a model where one differentiates between aggregate PCBM and dispersed PCBM, as is suggested to occur in real
devices after annealing[30, 4]. Another presented model will incorporate P3HT crystallization. The modified Ising model we present, which allows for morphologies directly extracted
from experimental data, is entirely new and has recently been published[31].
Chapter 6 presents the results of DMC device simulations performed on the modified
Ising model morphologies introduced in chapter 5 and compares these results to those calculated from more primitive morphological models.
We will also discuss the ability to incorporate additional experimental data into this type
of morphologically driven device simulation. As the contents of this chapter will summarize
the bulk of the work, the results presented in this chapter have been recently published[31].
Chapter 7 presents conclusions of this work and directions for future work.

7

Chapter 2
The Basics of OPV

2.1

Solar Cell Devices

The essential mechanism for the generation of power in a photovoltaic (PV) device (solar
cell) is unsurprisingly, the photovoltaic effect[32]. Similar to the photoelectric effect, wherein
photon induced excitation ejects electrons from a material’s surface, in the photovoltaic effect
an absorbed photon’s energy is transferred to an electron within the absorbing material. By
absorbing this energy, the electron is freed from it’s previously bound state such that it
becomes mobile inside the material[33]. The sudden absence of an electron in the absorbing
material generates a mobile vacancy site which acts as a positively charged particle in the
system, called a hole[34]. In solar cell design, the goal is to collect these photoexcited free
charges (both the electron and hole) at the contacts of the device where they may perform
work on an external circuit. In this way a solar cell is not a voltage source like a battery,
but rather a current source[35]. Not all photoexcited charges will successfully be collected
at the contacts however, as there are a number of loss mechanisms inherent in photovoltaic
devices.
As charges are being generated in oppositely charged pairs (one hole and one electron)
a solar cell requires a mechanism to separate these opposite charges and drive them to the
opposing contacts. Note that photocurrent is often defined as flowing in the reverse direction
of a circuit (electrons moving to the cathode and holes to the anode) as the device is supplying
8

vacuum
-

Ef

+
valence band

n-type

p-type

anode

cathode

condunction
band

Figure 2.1 Schematic of a simple PN junction in a solar cell at short-circuit condition, as
evidenced by the energy levels of the anode and cathode being equalized. Absorbed photon
excites an electron from the valence band into the conduction band. The electron and hole
can then follow the chemical potential gradient to the cathode and anode respectively. The
Fermi energy, Ef , exists between the valence and conduction bands. For interpretation of
the references to color in this and all other figures, the reader is referred to the electronic
version of this dissertation.
power rather than consuming it, hence the convention that the photocurrent generated in
a solar cell is negative. The principle physics inherent to conventional photovoltaic devices
are well illustrated using a band structure diagram, such as the example presented in figure
2.1. In order to photoexcite charges, absorbed photons must have sufficient energy to excite
an electron in the valence band across the band gap (Egap ) into the conduction band.
Once excited, the now mobile electron and hole will follow the chemical energy gradient of
the device. This gradient, represented by the slope of the bands in figure 2.1, may arise
from a number of sources such as the depletion zone formed in a pn-junction, a charge
population concentration gradient, the difference between the work functions of the contacts
when connected to an external circuit, or an applied external potential. In band diagram
representations electrons will flow down a potential gradient and holes (because of their
opposite charge) will flow up. Both charges travel along the bands until they are either
collected at the contacts or otherwise lost. The total power delivered has the familiar form

9

P = I × V,

(2.1)

where P is the power, I is the current delivered, and V is the voltage at which this current is
collected. This potential difference across the device from anode to cathode arises due to the
quasi-Fermi level separation of the charges[36]. The maximum delivered photocurrent under
no applied load is called the short-circuit current, ISC . As photocurrent will be proportional
to the total area of the solar cell upon which light is incident, it is often convenient to describe
the current density rather than the absolute current. The short-circuit current density, Jsc ,
may be described by the function

JSC = q

dEγ bs (Eγ ) QE(Eγ ),

(2.2)

where the integral is taken over Eγ , which spans all photon energies incident upon the
solar cell, q is the charge of an electron, bs (Eγ ) is the incident photon flux as a function
of photon energy, and QE is the quantum efficiency, which describes the probability that a
photon with given energy (Eγ ) incident upon the device will successfully deliver a charge to
the external circuit.
At the contacts themselves, a slight energetic mismatch will likely exist between the active
material of the solar cell (the p-type or n-type semiconductor) and the electrode. If the
transition between the two materials results in an energetic barrier inhibiting charges from
leaving the active layer to the contact, it is referred to as a Schottky barrier. Charge must
tunnel through or be thermally excited over a Schottky barrier in order to be collected[37].
This delay in collection time can lead to the build-up of a space-charge layer near the contacts,
which can further impede charge transport. By contrast, if the resultant barrier enhances
10

charge transport across the barrier by nature of being at a lower energy state for the charge
relative to it’s energy in the semiconductor, the contact is referred to as Ohmic[38, 39].
In any device with a quantum efficiency of less than unity, not all photogenerated charge
is being successfully collected at the contacts. The three main loss mechanisms are radiative
recombination, Auger recombination, and trap assisted recombination[35, 40]. Radiative recombination occurs when an electron relaxes across the band gap, emitting an associated
photon of energy Egap . Radiative recombination is a common loss mechanism in a conventional solar cell, with it’s rate being governed by the width of the band gap. Auger
recombination also involves relaxation of a charge carrier across the bandgap, but the energy
released acts to push a different charge carrier to a higher energy level in the same band.
In terms of PV devices, this sacrifices one charge with no particular energetic gain for the
other, as the charge excited above the band edge will quickly thermalize back to it’s previous
energy level at the bandgap edge. Auger recombination is only a significant loss mechanism
in systems with a high free charge density, which rarely occurs in conventional devices under normal operating conditions. Trap assisted recombination occurs when charges become
trapped in energetic defects (traps)[41] which have an intermediate energy in the band gap,
such that the trapped charge gives up this difference in energy in the form of phonon radiation (heat). Trap assisted recombination is often the dominant form of loss near the surface
of a device, due to defects introduced by the bulk termination at the surface[38].
The potential difference applied by an external load will cause external charges present in
the contacts to be injected into the active layer of the solar cell. The injected charges flow in
the opposite direction to the flow of photogenerated charge current towards the contacts[35].
This injected opposing flow of charges is referred to as the dark current density, Jdark , as it
will be present whether the solar cell is under illumination or not. Though “dark” charges
11

(a)

(b)

(c)

(d)

Figure 2.2 Demonstration of applied voltage to a simple solar cell band diagram. As the
open circuit voltage, Voc is passed, the runaway diode behavior emerges.

12

are continuously thermally injected from the contacts into the solar cell[37], because of the
band structure design of conventional PV devices they will likely be unable to travel across
the device until sufficient external bias is applied, as is demonstrated in figure 2.2. This
diode-like behavior of the dark current is a common feature in solar cells, principally due to
the rectifying features of their internal energetics favoring charge transport in one direction.
One may describe the dark current density, Jdark , with an ideal diode equation

Jdark = Jo eqV /kB T − 1 ,

(2.3)

where Jo is a constant, V is the net voltage across the device, and kB T is the Boltzmann
constant times the temperature.
The equivalent circuit for a solar cell, seen in figure 2.3, demonstrates a photon induced
current source (JSC ), a dark current contribution modeled as a diode (Jdark ), additional
current losses modeled as a parallel and series resistance (Rp and Rs respectively), and the
contact points across which a solar cell would be connected to an external circuit with a
potential difference when a load is applied (V ).
Though a solar cell will deliver the maximum photocurrent at short-circuit, it is often
desirable to apply an additional external voltage which acts to oppose the flow of the photocharges. The reason for this applied voltage, Vapp , is to increase the total power delivered
by increasing the potential at which the photocharges are collected at the contacts. At a
certain threshold voltage, however, the flow of the photocurrent will be stopped entirely
rendering the solar cell useless. Because of the desire to maximize power delivered by the
device, solar cells are often evaluated with a current versus voltage (IV-curve) or current
density versus voltage (JV-curve), an example of which is seen in figure 2.4.

13

Figure 2.3 Equivalent circuit for a solar cell, showing photon absorption generating a shortcircuit current density, JSC , which has losses represented as a parallel and series resistances,
RP and RS respectively. Dark current, Jdark , is represented as a diode in the reverse
direction of the photogenerated current.
Though figure 2.4 is a function of applied voltage, Vapp , recall that the total voltage
across the device will be the combination of the applied voltage and the built-in bias, Vbi ,
such that the total voltage across the system is denoted Vnet = Vbi − Vapp . Three important
points of note on any JV-curve are the short-circuit current density (Jsc ), the open-circuit
voltage (Voc ), and the maximum power point[35]. The short-circuit current density describes
how the system behaves when no external voltage is applied (Vapp = 0) and is the maximum
possible collected current in the device. The open-circuit voltage is the point at which the
applied voltage has caused the net potential across the cell to be zero[42], at which point the
net current inside the device will be driven only by charge diffusion, which is typically much
smaller then the contribution of the drift current. The maximum power point is the location
on the JV-curve where total delivered power will be maximized, at voltage Vm and current
density Jm , which is subsequently the point of optimal device operation. In any regime of
the JV-curve outside this quadrant, either with applied voltages greater than Voc or under a
reverse bias with Vapp < 0, the solar cell will consume rather than generate power. A solar
14

Figure 2.4 Example of a current density (J) versus applied voltage (Vapplied ) graph,
demonstrating the diode like behavior of a solar cell. Curves are shown for both a device in the dark (Jdark ) and under illumination (Jlight ). Also shown on the graph are the
open circuit voltage, VOC , short-circuit current, JSC , thick dashed line to better understand
the fill factor, F F . Notable behavior points occur at short circuit, a , near maximum power,
b , near open circuit, c , and finally at runaway diode-like behavior, d .
cell’s ability to “fill-out” the area of it’s JV curve to the maximum possible power, Jsc × Voc ,
is described by the fill factor, F F , defined

FF =

Jm Vm
.
Jsc Voc

(2.4)

The total efficiency, η, of a solar cell describes it’s ability to convert incident light power
into usable work, defined as

η=

Jm V m
Jsc Voc FF
=
,
Ps
Ps

(2.5)

where Ps is the incident power density, delivered by the photons, to the solar cell.
Photogenerated free charges will transport throughout the solar cell by a combination
15

of electrostatic field induced drift and concentration gradient driven diffusion. The current
densities in a solar cell at equilibrium can thus be separated into a drift component and a
diffusion component, described (in one dimension) as[35]

Jn = qDn n + qµn Fn

(2.6)

Jp = −qDp p + qµp Fp ,

(2.7)

where Jn/p is the net current density for electrons/holes driven by the diffusion due to
the gradient of the charge population density,

n/p, the charge mobilities, µn/p , and the

diffusion constants, Dn/p , through the use of the Einstein relationship

µ=

qD
.
kB T

(2.8)

Note that this approximation of current is only valid in the regime where charge populations are at thermal equilibrium with the Fermi level, the material has a relatively low
density of defects or traps, the charge population is well described by Boltzmann statistics
(non-degenerate semi-conductors), and the gradient of the effective field (band gradient) is
relatively constant over the device[43]. This simple model does not explicitly incorporate
any morphological information, utilizing only an average mobility and charge population to
describe the photocurrent.

16

2.2

Organic Photovoltaics

Organic photovoltaic (OPV) devices have great promise for applications unsuitable for conventional solar cells, due to their light-weight and potentially less expensive fabrication
costs[15, 13, 2]. The difference in underlying mechanics and design of OPV systems compared to their inorganic counterparts stems from the fact that the energetics in an OPV
systems arise from disordered molecular systems, rather than continuous band structure
seen in conventional semiconductors. These molecular materials cannot strictly be thought
of as p-type and n-type semiconductors, as there are no inherent free charges to be exchanged and thus no depletion region will will form at pn-junctions. Instead, OPV systems
are described as a combination of acceptor and donor materials[44], requiring their energetic
design considered through molecular based HOMO (highest occupied molecular orbital) and
LUMO (lowest unoccupied molecular level) levels, rather than valence bands and conduction
bands. An example of a molecular based OPV HOMO/LUMO level diagram is shown in
figure 2.5. In these materials, absorbed light will not immediately generate free charges as
in conventional PV, but instead generate strongly bound Frenkel excitons. These photogenerated excitons must be dissociated prior to exciton decay if the charges are to be collected
at the contacts.
The exciton is a bound state of an electron and hole pair, Coulombically bound to one
another resulting in a bosonic particle (being the combination of two spin one-half particles).
In these materials, excitons exist for a very brief period of time compared to free charges
(on the order of picoseconds) before they decay[45], with the energy lost through radiative
(photon emission) or non-radiative (phonon emission) processes. On a band diagram, excitons exist as states slightly inside the bandgap, such that they are at a lower energy state

17

than the two free charges unbound from one another. In most conventional photovoltaics,
photogenerated excitonic states are readily separated via thermal excitation or internal field
gradients, such that most models of conventional solar cells do not even take the exciton into
account.
Unlike conventional PV, the internal energetics of an OPV device are often insufficient to
naturally facilitate exciton dissociation, requiring an additional energetic asymmetry to do
so. This asymmetry can be provided by the energetic contrast present at interfaces between
acceptor and donor material, called a heterojunction[35]. This naturally suggests fabrication
of devices in which exciton is likely to discover, via it’s drift motion, a heterojunction interface
prior to decay. The most widely studied morphology fitting this requirement is the bulkheterojunction[46] (BHJ), in which the active layer of the device is composed of a finely
mixed blend of acceptor and donor material such that the length scales of the absorbing
material domains are of the same order as the exciton decay length.
Once an exciton reaches a BHJ interface, the difference in energies of the HOMO/LUMO
levels between the acceptor and donor materials is sufficient to separate the tightly bound
exciton into it free charge carriers. In these OPV systems, the dissociated electron is considered bound inside the acceptor material, as transport back to the donor would involve
passing into a higher energy state and is therefore unfavorable. Similarly, the hole is bound
inside the donor material. The limitation of the BHJ design is that while it greatly enhances
exciton dissociation efficiencies, it decreases the overall charge collection efficiency[2, 47] due
to the tortuous path a free charge is forced to follow (electrons through the acceptor material
only, and holes through donor only) to reach the contacts. Optimization of the internal morphology such that it satisfies the transport requirements of both excitons and free charges is
thereby a principle concern in the development of higher efficiency OPV devices.
18

-

transport

anode

cathode

n

Voc

+

oto
ph

excitation

-

+

d

a

Figure 2.5 Example of a molecular PV device, showing photoexcitation occurring in the
donor material, then charge separation and transport to the contacts. Note the lack of any
band bending, but rather hopping between HOMO levels (hole transport) and LUMO levels
(electron transport).
Free charges in molecular devices are typically more localized than their inorganic counterparts, due to the energetic disorder inherent in the complex disordered energies and morphologies associated with OPV. This leads to charge mobilities in OPV devices which are
orders of magnitude smaller than those seen in their inorganic counterparts. Additionally,
charge transport is often modeled as a hopping process as compared to a continuous current,
with the hopping charges jumping between localized minima energy states on timescales
associated with random energetic fluctuations inherent in the system, as is conceptually
demonstrated in figure 2.6.
Free charge recombination is arguably a greater concern in organic systems, as evidenced
by their lower quantum efficiencies compared to inorganic devices. This is in part due to the
lower mobilities preventing free charge pairs from escaping the Coulombic attraction that
exists between them. Additionally, the lack of free charge within molecular materials leads to
a lack of any significant Debye screening effects. Consequently, the only energetic screening in

19

Figure 2.6 Qualitative illustration of how a charge hops on a length scale related to thermal
vibration of a local energy landscape.
these devices originates from the fluctuating random thermal contributions to the energetic
landscape, described by a Bjerrum length. Free charge recombination is often specified
as being either geminate or bimolecular in nature. Geminate pair recombination refers to
a pair of charges that were generated together (from the same photon absorption/exciton
dissociation event) recombining and thus being lost. Bimolecular recombination refers to
a charge which had successfully escaped it’s initial partner but then recombines with an
opposite charge some time later. Because of the low carrier densities present under normal
illumination conditions, geminate recombination seems likely the dominant form of charge
recombination, as opposed to bimolecular.

20

2.3

P3HT/PCBM OPV Device Design

The devices we model in this thesis are based on the work of Jonathan Kiel, of Michael
Mackay’s group at the University of Delaware. These solar cells were inspired by the popular
polymer/nanoparticle bulk heterojunction design[48], which uses the conductive polymer,
poly(3-hexylthiophene) (P3HT) as a donor material and the conductive nanoparticle, [6,6]phenyl-C61-butyric acid methyl ester (PCBM) as an acceptor material. A simple schematic
for the device is shown in figure 2.7.
The active layer of these devices is composed of a 1-to-1 by volume mixture of P3HT and
PCBM, which is then spin coated onto a prepared front contact. The front contact must be
transparent to incident photons, so that the photons may reach and be absorbed in the active
layer. Indium tin oxide (ITO) is a common choice for the front contact material, but other
common materials include fluorine doped tin oxide (FT) or doped zinc oxide. Prior to spin
coating, the front contact is often coated with an electron blocking layer of poly(ethylene
dioxythiophene):polystyrene sulfonate (PEDOT:PSS) so that only holes will be collected at
the anode (ITO) contact. Post spin cast annealement of the active layer was found to cause
an overall improvement in the final device performance[8, 25]. The back contact is typically
a metal, often aluminum, which is carefully evaporated onto the active layer.
In these P3HT/PCBM devices, photons are absorbed only into the polymer (P3HT) generating excitons. These excitons diffuse, with hop rates influenced by local site energies[49],
until they either discover a bulk heterojunction (BHJ) interface or decay. Though some
evidence suggests that decaying excitons emit photons which may later be reabsorbed in
the polymer and form a new exciton, for the purposes of our work we consider any decayed
exciton as lost. This assumption is justified because the decay-generated photons would

21

likely have wavelengths just below the HOMO-LUMO gap energy, and thus the likelihood
of re-absorption in these materials would be small.
Dissociated excitons form a pair of charges which are localized within separate materials
(electrons in the acceptor and holes in the donor) but are still strongly attracted to one
another via their Coulomb potential[50]. As the density of photogenerated free charges is
relatively low in these devices under normal operating conditions, no screening due to local
charge induced fields need be considered. We consider the attraction between these charges
to be a function of only their separation distance and the difference in dielectric constants of
the acceptor and donor materials. We may write the potential, φ, experienced by a charge
due to it’s geminate partner as

φ(r, z) =

2q
( 2 + 1)

1
r

−

( 1 − 2) q
( 2 + 1) 1

1
2z

,

(2.9)

where q is the charge sitting in a material of relative dielectric constant 1 at a distance of
z from the interfacial boundary, while it’s geminate partner (having opposite charge) resides
in the opposing material of relative dielectric 2 , at a total distance of r from the charge.
The full derivation of equation 2.9 is provided in appendix A.1. As the force on a charge q1
by charge q2 is given F = −q1 ∂φ , this will generate an overall attractive force due to the
∂r
opposing charge, given in the first term of equation 2.9. The second term of equation 2.9
will either be positive or negative, dependent on whether the image charge is generated in a
material of lower relative dielectric constant (resulting in a repulsive force) or higher relative
dielectric (resulting in an attractive force).
If the relative dielectric constants are equal, such that 1 = 2 = , then equation 2.9 will
simplify to the familiar form

22

φ(r, z) =

2q
( + )

=

2q
(2 )

=

1
r

−

( − )q
( + )

q
.
r

1
2z

1
r
(2.10)

Conformational variations and deep energy traps lead to a Gaussian like density of states,
σr , throughout the transporting medium with some width [43] which for these devices has
assumed values around σr ≈ 1 × 10−20 J ≈ 1/16 eV [51]. If the geminate pair separates to
a distance such that the Coulomb potential between the charges is less than the energetic
width of density of states, we may consider a pair now fully separated free charges. This
is not to say charges reaching this range are certain to be collected at the contacts, only
that any further recombination is to be considered bimolecular rather than geminate. The
average size of a geminate pair will thus be determined by the dielectric constants of the
materials and the hop rate of the charges.
Free charges are considered localized to 1 nm in these systems, with a single charge confined in a thermally oscillating energetic landscape. Charge transport is driven by a hopping
mechanism well described by Marcus theory[51, 43, 52], which will be further discussed in
section 3.3 of this thesis. Free charges within the confines of their respective material (electrons in acceptor and holes in donor) will hop throughout the morphology until they are
either collected at the contacts or lost due to recombination.

23

Figure 2.7 Schematic of a P3HT/PCBM based OPV device, with the physical layers shown
in the top figure and the corresponding energy levels on the bottom. In this diagram, light
passes in from the right hand side and is absorbed in the active layer. Electrons travel to
the back metal contact (typically aluminum) and holes, passing through the PEDOT:PSS,
are collected at the top contact (ITO).

24

Chapter 3
Monte Carlo Device Model for OPV
In this chapter, we begin by summarizing the different scales of modeling used in OPV
simulations and giving a background of the calculations that have been performed. We
motivate our use of the dynamic Monte Carlo (DMC) model for use in exploring the effects
of morphology on OPV device performance, summarizing how the DMC method has been
used previously. We will explain the DMC model in detail as well as the origin of the event
rates used for both excitons and charges. We will demonstrate an application of the DMC
model on simple bilayer systems with a number of different calculations, testing both the
effect of system energetics and morphology. Though we ultimately modified the presented
DMC method to better fit our focus on the direct effects of morphology on transport (see
chapter 6), we present the standard method in full detail here.

3.1

The Scales of Modeling and Literature Review

Modeling of organic photovoltaic devices is useful to both better understand how internal
mechanisms effect device performance and to help design higher efficiency devices[53]. The
principle factors involved in an OPV device’s performance, and thus the factors under investigation by simulation, may include photon absorption and conversion, exciton generation
and transport, charge generation and transport, carrier loss mechanisms, trapping, internal
field effects, HOMO/LUMO levels, charge extraction at contacts, and dark charge injec-

25

tion. With so many different mechanisms influencing OPV device performance, no single
model exists which simultaneously incorporates all these mechanisms on a device scale[54].
The models used in OPV device simulations can be separated into three scales, macroscopic
models, mesoscopic models, and atomistic models.
Atomistic models are those which occur on the smallest length scale (1 ˚ - 5 nm) and are
A
the most precise calculations of the processes occurring inside the device, often treating these
processes with the most physically precise theories. Examples of atomistic scale calculations
include quantum mechanical calculations of exciton transport dynamics across polymers[55],
photon absorption and localized exciton mechanics [56], calculations of charge transfer or
separation across a heterojunction[57, 58, 59, 60, 61], self-consisted field theory studies of
polymer phases[62], and Monte Carlo calculations of charge dynamics across HOMO/LUMO
gaps[43, 63]. Simulations on the atomistic scale, while more accurate than larger scale
models, are computationally infeasible to implement on length scales required for full device
simulations.
Macroscopic models are those which simulate the largest (100nm−1µm) and subsequently
coarsest length scales. These models are capable of simulating an entire device, but to do
so must approximate the internal features in some way. Additionally, these models rarely
simulate individual charges, instead modeling charge densities and net currents and often
producing full IV characteristics for a device. One of the most common types of macroscopic
models are the so called drift-diffusion calculations[64, 65, 28, 66, 67], which examine the
influence of effects such as trap states[68, 69], mobility dependence[70], material mixture
ratios[71], and doping effects[72]. Other examples of macroscopic models include simple
reaction rate calculations[73] or detailed balance studies[74].
While atomistic scale models are the most physically accurate, their small scale pre26

vents their use in exploring device-spanning phenomena. On the other end of the spectrum,
macroscopic models are quite capable of full device characterization, but only when utilizing
simplified morphological models which may wash out detailed nanoscale features. We feel a
third scale of simulation, the mesoscopic models (2nm − 400nm), are ideal for investigating
the influence of nanoscale morphology on a wide breadth of physical mechanisms of interest
in these OPV systems.
Mesoscopic models often simulate individual charge and exciton transport over a sample
or cross-section of the total device. The transport mechanisms in mesoscale models are often
abbreviated or simplified from those used in atomistic models. Many mesoscopic models
utilize the Dynamic Monte Carlo (DMC) method and include studies of the influence of
nanoscale structure on device performance[2, 51, 1, 75, 76, 77]. These models are readily
expanded to include many different effects, such as energetic disorder[27, 78, 79, 80], charge
separation and recombination[81, 82], surface wetting effects[83], localized trapping[84], electrode heterogeneity[85], material roughness[86], and light intensity[87]. Because of their
simplified transport mechanisms, mesoscopic scale models inherently lack the precise detail
provided by atomistic scale calculations. However, because of the much broader view accessible by models of this scale, the influence of internal nanoscale features may be explored
over lengths associated with charge lifetimes. We believe that mesoscale models, specifically
the dynamic Monte Carlo model, offer an ideal tool with which to bridge the information gap
between the accuracy of atomistic scale calculations and the macroscopic scale calculations
required for efficient full device characterization.

27

3.2

Basics of the Dynamic Monte Carlo OPV Model

In it’s simplest form, the dynamic Monte Carlo (DMC) device model simulates a solar cell’s
active layer, a nanoscale morphology composed of two materials: an electron acceptor (such
as PCBM) and donor (such as P3HT)[88]. The precise structure is defined on a cubic lattice,
with each site being defined as one material or the other. System dynamics are simulated via
modeling the transport of individual excitons, electrons, and holes; all of which are modeled
as single lattice-site occupying particles. The DMC model is viable in both two[89] and three
dimensions[2], though we will focus on to the three dimensional case for this thesis. The
assumption is made that free electrons are bound inside the acceptor material, and free holes
are bound inside the donor material. This assumption is justified because of the difference
in energy levels between the HOMO/LUMO levels of the acceptor and donor material, as
seen in figure 2.7.
Excitons are generated randomly throughout the bulk of the morphology’s donor sites,
and diffuse until they either discover a heterojunction interface (a donor site adjacent to an
acceptor site) where they may dissociate into free charges or decay. Free charges are generated
in electron-hole pairs, and proceed to drift/diffuse throughout the morphology until they
either discover an electrode and are collected, or are lost due to recombination. Real devices
typically have an electron blocking layer, such as PEDOT:PSS, between the active layer
and the hole collecting electrode in order to prevent photogenerated electron injection into
the wrong contact (as this would generate current flowing in the wrong direction, and thus
consuming power). In this simple model, we approximate this by only allowing charges to
be collected at their respective contacts (electrons at the cathode, holes at the anode).
The model morphology is mapped onto a cubic lattice with all sites being assigned as

28

either pure acceptor or pure donor, with a typical lattice spacing of al = 1 to 3 nm. Periodic boundary conditions are considered to occur in the x and y (horizontal) directions,
with reflective boundaries assumed in the z (vertical) direction. Additionally, electrodes are
assumed at the edges of the z direction, with the electron collecting electrode (a metal, such
as aluminum) at the air interface (z = 0) and the hole collecting electrode (transparent conductive contact, such as ITO) on the substrate (z = Lz ). Due to the assumed confinement
of electrons in acceptor material and holes in the donor material, there must exist a nearest
neighbor percolative path from any site to the corresponding contact if charges generated
upon such a site are to be collected.

3.2.1

First Reaction Method

The first reaction method (FRM) is a common algorithm used in Monte Carlo transport
simulations of OPV systems[2, 47], and is useful in efficiently handling a number of dynamic
particle trajectories simultaneously. The basic premise is that one keeps a chronologically
ordered queue of future events and their associated times. Every possible event in the
system (such as charge movement, recombination, collection at electrode, etc) is assigned
an associated rate. During each simulation step, the top event is removed from the queue
and, so long as that event is still possible, executed. After event execution, all candidate
next events for the involved particles are considered, with associated times generated from
a characteristic rate as

−1
τevent =
ln X,
ωevent

(3.1)

where τevent is the candidate event time, ωevent is the associated characteristic event

29

rate, and X is a uniform random number between 0 and 1. Equation 3.1 is taken directly
from the Gillespie algorithm[90], and will produce stoichiometrically correct results for the
rate given while still being dynamically unpredictable. For example, an event with a characteristic rate of ωa = 20 would be selected 4 times more often than an event with rate
ωb = 5. All candidate events have associated times generated in this manner, but only the
candidate event which produced the fastest time is selected and inserted into the queue at
it’s chronologically associated position.

3.2.2

Exciton Events

For excitons in our model, possible events and their associated rates include exciton generation on a donor (P3HT) site due to light absorption (νx-gen ), hopping movement to any
neighboring site (νx-hop ), exciton decay (νx-decay ), and dissociation into free charges when
adjacent to a heterojunction (HJ) site (νx-diss ). Commonly used values are shown in table
3.1.

3.2.2.1

DMC exciton generation

As excitons in this model are generated due to photon absorption, one could apply the BeerLambert law in order calculate the exciton generation rate. The Beer-Lambert law describes
photon absorption in a material with an exponential decay as a function of thickness and
density of the material, cumulatively describing a path-length. A simple derivation for
the Beer-Lambert law begins by dividing the absorbing sample into a series of thin slices,
perpendicular to the incident light.
The assumption is made that the beam intensity decreases as it passes through each slice
of the material due to molecular absorption and that all particles absorbing light in the slab
30

have a perpendicular cross-section to the beam, σ, such that any light not absorbed in a
thin slab passes through. The thickness of a layer, dz, is considered to be thin enough that
one absorbing particle does not obscure the interaction of the photons with another particle.
As each slab has a definable area, A, and we assume each absorbing particle has a specific
cross-section, we define each slab to have a particle density, n. The fraction of photons
absorbed by a specific slab will be proportional to it’s opaque surface area (σAN dz) divided
by the area of the slab such that one may write the intensity absorbed by a given slab as

dIz = −σnIz dz,

(3.2)

where we specified the intensity lost, Iz , due to a specific layer, dz. After rearranging
and integrating, we come to the equation

ln [Iz ] = −σnz + C,

(3.3)

where C is a constant due to integration. If we assume the difference in intensity for a
sample of thickness, l, is described such that it has intensity I0 at the surface (z = 0) and
final intensity Il at the bottom of the sample z = l, then the difference in intensity across
the sample may be written

ln [Il ] − ln [I0 ] = (−σln + C) − (0 + c) = −σln.

(3.4)

Defining the transmission through a sample, T , as the fraction of light intensity at the
bottom of the sample (Il ) over the intensity at the surface (I0 ), we see the expected exponential decay of intensity, as a function of thickness and other absorption parameters

31

I
T = l = e−σln = e−αl ,
I0

(3.5)

where we have combined the absorption parameters of number of particles and particle
cross-section into one parameter, α.
For many thin film OPV devices, however, the photon absorbing active layer is so thin (on
the order of 100 nm) that the Beer-Lambert law is not required to accurately simulate device
performance, as photon absorption may be treated as a constant throughout the thickness
of the active layer[2].
In the simple DMC OPV simulation, the exciton generation rate, νx-gen , is defined

νx-gen = Iincident αx Lx Ly ,

(3.6)

in which Iincident is the flux of photons incident upon the cell based upon the AM1.5
spectrum, αx is a material constant representing what percentage of photons incident upon
the device are absorbed (assumed 1 in the DMC simulations), and Lx/y are the simulated
morphologies dimensions in x and y, producing incident surface area. As noted previously,
our model assumes thin morphologies where absorption of photons does not have a depth
dependence.

3.2.2.2

DMC exciton hop rates

Once an exciton is generated, it proceeds to diffuse throughout the device with a rate νx-hop
until it decays or dissociates. The decay rate is based on experimentally observed lifetimes of
νx-decay = 1/τx-lifetime = 500 ps−1 . Note that due to the scales of morphological area being simulated in the DMC model (1000 nm2 to .16 µm2 ), and given that νx-gen
32

νx-decay ,

only rarely will there exist more than a single exciton in the model at any given time. Experimental measurements suggest exciton diffusion lengths on the order of ≈ 10nm. For
these P3HT/PCBM systems, we set the exciton hop rate to reproduce the desired mobility
through the use of the diffusion equation and the Einstein relationship. A detailed explanation on the selection of hop rates to achieve a desired mobility is found in appendix A.3 of
this thesis, but briefly the Einstein relation is defined

D=

q
kB T

µ,

(3.7)

which relates the diffusion constant, D, to the mobility, µ, via the charge, q, temperature,
T , and Boltzmann’s constant, kb . The relationship between the diffusion constant and the
average hop rate and average rate of a particle on a cubic lattice is defined

D = l2 ωi /6,

(3.8)

which gives the average hop distance, l, and average hop rate, ωi , required to generate a
specific diffusion constant. Note that the derivation of equation 3.8 is also found in appendix
A.3. Taking the average hop length used in our DMC model to be the lattice constant (l = a)
we may combine equations 3.7 and 3.8 to to final form

ωi =

3.2.2.3

6kB T µ
.
qa2

(3.9)

DMC exciton dissociation rate

Exciton dissociation may only occur when the exciton is directly adjacent to a heterojunction,
such that a geminate pair may form across the two materials. We take the dissociation rate
33

to be much faster than all other processes when the exciton is adjacent to a HJ (on the order
of 1000 × νx-hop ), but 0 otherwise. If an exciton is simultaneously adjacent to several HJ
interfaces (for instance, in a corner), the target interface to dissociate across is selected at
random.

3.2.3

Charge Events

Although charge diffusion can be correctly modeled in bulk material by selecting a single
time through the use of 3.1 and then assigning a random hop direction, we simulate the drift
behavior due to internal fields by dynamically modifying hop rates as particles move through
the morphology. These fields originate from a wide variety of sources including built-in bias
due to differences in work function, ∆φ, applied field, Vapplied , coulombic interaction with
nearby free charges, Ecoulomb , image charge effects near contacts, Eimage-c , and inherent
energetic disorder[91] throughout the morphology, Edisorder . In the DMC model, these
various field effects are maintained on a set of field maps, Ee/h/x (i, j, k), which specify the
current potential value of site (i, j, k). For convenience, we separate field effects into three
different maps for electrons, hole, and excitons (as Ee , Eh , and Ex respectively).
The random site energies, Edisorder simulate the effect of thermally induced energetic
disorder in the systems[92, 93], and are assigned at the start of any simulation with a
Gaussian distribution centered at 0 with a width of σrand . The net field across the device,
being the difference of the built in bias and applied bias (Vnet = Vbi − V app) is applied as a
static linear field gradient, Ee/h , across the system.
A free charge occupying site (i , j , k ) will coulombically contribute to nearby site (i, j, k)
with potential

34

V (i, j, k) =

q2

1

r R(i ,j ,k )

,

(3.10)

where V (i, j, k) is the contribution to the field, 1 is the elementary charge, r is the
relative permittivity of the material, and R(i , j , k ) is the distance from site (i, j, k) to the
charge generating the field. V (i, j, k) is added to the field map of the same-type charge
(being a repulsive potential) and subtracted from the field map of the opposite type charge
(being an attractive potential).
Image charge effects may occur near the contacts, and are treated by generating an
equivalent field as one due to a charge placed at the site (i, j, −k), that has the opposite
polarity as the original charge at site (i, j, k). In these systems, the image charge will act
to enhance charge collection at the contacts by attracting free charges, though this is a
function of material permittivity and is further discussed in appendix A.1 of this thesis. In
DMC models of these BHJ systems, the permittivity of the two materials involved are often
taken to be the same such that no image charge effects at the heterojunction itself need be
considered. This is an acceptable approximation, as the charge recombination rate is already
being fit to experimental results, which likely encompasses any enhanced recombination due
to image charges effects across the heterojunction.
While exciton dissociation generates a pair of free charges across a heterojunction, these
charges are still somewhat bound to one another due to the strong coulombic attraction,
and are thus referred to as a geminate pair. Though the precise mechanism of geminate
pair separation is somewhat debatable, from the standpoint of the DMC model a charge
must escape the energetic influence of it’s geminate partner in order to become a truly free
charge. The Debye length is commonly used in simulations of conventional PV devices for

35

the effective charge field cutoff length as it is brought about by the Coulombic fields of nearby
free charge, and previous OPV simulations have used the Debye length in a similar way[2].
However, as OPV systems have a relatively low density of thermally excited free charge, it
is not entirely correct to characterize geminate pair separation in this way. In discussions
of OPV devices it is more fitting to use the Bjerrum length as a field cutoff radius, being
the separation distance at which the electrostatic interaction between a pair of charges is of
equal magnitude to the local thermal excitations[94]. The Bjerrum length, λB , is defined as

λB =

q2
,
r kB T

(3.11)

where q is the elementary charge, r is the relative dielectric constant, T is the temperature. In our simulations, if a pair of charges is at distance λB or greater, we may consider
their electrostatic interaction effectively masked by natural thermal diffusive processes, and
as such, we need not calculate the effects of any single charge generated Coulombic fields
beyond the length λB . This allows for significant speed increases in the calculation.
Both the built in bias (∆φ) and the applied field (Eapplied ) are treated as static throughout the simulation. In most cases, this net static field is treated independently of BHJ morphology, because of the assumption that the relative permittivity, and thus relative dielectric
constant, are identical in the acceptor and donor. The means the internal static field is typically treated as a simple linear potential gradient across the morphology. In reality, the field
may be more complex, particularly near the contacts, but for most simple models a linear
field approximation is used. The common notation is that a positive ∆φ generates a field
which enhances charge transport to the charge’s respective contact (electrons to the anode,
holes to the cathode), while a positive Eapplied field acts to oppose charge transport. As

36

such, we write an effective vertically linear static field, Estatic

k
Ee-static (k) = Eapplied − ∆φ L
z

Eh-static (k) = Eapplied − ∆φ

Lz −k
Lz

:for electrons
(3.12)
:for holes,

where Lz is the largest lattice site value in the z-coordinate representing the substrate
interface, which is modeled as the hole collecting contact. At high applied field, diode like
behavior emerges. As has been previously discussed in section 2.2, this behavior is expected
as a high applied field will sweep the dark-charges being injected from the contacts across
the device to the opposing contact, thus providing a source for the runaway current in the
non-desirable direction.

3.2.3.1

DMC charge hop rates

Charge hop rates, νc-hop , are calculated dynamically between all nearest neighbor sites, m
and n using a simplified form of the Marcus rate equation

νc-hop

mn

= νc exp

−(Em − En + ER )2
,
4kB T ER

(3.13)

where νc is a constant, Em and En are the energies of the sites based on the sum of static
and dynamic site energies, ER is the reconfiguration energy, and kB T is the simulation
temperature. Note that the Marcus rate equation will be discussed in further detail in
section 3.3. Careful consideration must be taken in the selection of the prefactor hopping
constant, νc , such that the desired diffusive rate is produced. Additionally, the desired
effects of confinement on charge behavior must be considered, which are discussed in detail
in appendix A.3. For purposes of the simple model it is sufficient to calculate the hop rate
to all surrounding sites regardless of site assignment (acceptor or donor) and only later,
37

as the selected hop execution, reject movement to forbidden sites (hops out of confining
material). This mechanism will model the situation where charges that have become caught
up in tortuous morphologies have difficulty finding an acceptable hop site, and thus their
effective mobility will drop.

3.2.3.2

DMC charge collection and dark charge injection

The charge collection rate at electrodes, νc-collect is considered fast compared to all other
processes due to image charge effects, on the order of 100 × νc-hop .
Dark charge injection is modeled using rates derived from theoretical thermionic emission
from the contacts. The total electron current over a Schottky barrier can be described[39]

Jdark = A∗ T 2 exp −

qφB
,
kB T

(3.14)

where φB is the barrier height and the Richardson constant, A∗ , can be approximated
(neglecting tunneling and reflection)
∗ 2
∗ = 4πqm kB ,
A
h3

(3.15)

where m∗ is the effective mass. In the DMC simulations, dark charge injection is simulated by generating either an electron or a hole (chosen at random) at a random site of
the confining material (i.e. electrons in acceptor sites) which is adjacent to the injecting
contact. The simulation injection rate is based on the current density from equation 3.14
and the area of the contact being simulated. Under typical OPV operating conditions, dark
charges will almost always be recollected at their respective contacts. If the applied voltage
becomes great enough, however, the dark charges will be swept across the device and col38

lected at their non-corresponding contact. This mechanism is the principle reason why dark
charge injection must be modeled in order to calculate full JV-characteristics in a model
morphology, as they provide the runaway current seen at high applied voltage as in figure
2.4.

3.2.3.3

DMC charge recombination rate

The charge recombination rate, νc-recomb , is typically a fit parameter in these DMC simulations with a value selected such that the simulation reproduces actual device performance.
Charge recombination is only a candidate event when an electron and hole are directly adjacent to each other at a HJ site. Typical values used for νc-recomb are on the order of 100 s−1 ,
which is significantly slower than any other rate used in the simulation. The very slow rate of
νc-recomb can be understood because, despite referring to them as “free charges”, recently
generated geminate pairs are anything but free. The coulomb attraction, defined in equation
3.10, is typically much stronger than any other effect geminate charge pairs will encounter.
This tends to lead to the behavior of two “free” charges actually spending a majority of their
simulated existence bound to each other across an interface of the HJ. This pair of charges is
unlikely to drift apart, even when the interface runs parallel with the built in field. Because
of this, charge recombination rates must be low otherwise most geminate pairs would simply
recombine before ever having a chance to escape one another. This mechanism can be considered one of the weaknesses of this simple model. Note that if the selected next event for a
charge is a recombination action, the opposite charge must still exist and be present at the
time of recombination, and a new action will be selected via the standard FRM procedure.

39

3.3

Marcus Hop Rates

As charge hopping transport in this model is derived directly from Marcus theory, we will
briefly explain it’s origin in this section. Marcus theory was developed to explain the rates
associated with electron transfer reactions. In the case of the DMC model employed here,
the transfer reaction represents a charge hopping between localized energy sites. One may
conceptualize the energetic landscape inside this system as a series of local energy minimas
with the barriers between sites oscillating due to thermal excitations, as was illustrated in
figure 2.6. The motion of the localized charge carrier itself is a function of the charge’s
capability to hop between localized minima (lattice sites in the DMC model), and is limited
by the details of the current energy landscape. The oscillations are essentially phononmitigated interactions with nearby sites, such that when a potential crossing path does
come about (at some specific energy level) the charge may readily hop to an adjacent site as
charge transport is presumed to be much faster than the reorganization speed of the energetic
landscape. It is this subtle phonon interaction between the system’s energetic landscape that
creates the possibility for charge transfer, and not the reverse.
In order to calculate a hop rate, the Gibb’s free energy of activation, ∆Go , is used as
it represents the minimum energy required to excite a charge out of one localized energy
site and into another. The Marcus transition is conceptually visualized in figure 3.1 as a
function of free energy versus the reaction coordinate[52]. The reaction coordinate is an
abstract variable which describes the energy path involved in a reaction, with each point
along the axis corresponding to a different configuration of molecular distances and angles
related to the reactants and bonds. The standard thermally activated rate equation, known
as the Annherius equation[95], is given

40

Figure 3.1 Parabola geometry used to derive Marcus rate activation energy. The net energy change from the reaction is the Gibb’s free energy of reaction, ∆Go . λ is energy of
reorganization. Gibb’s free energy of activation is found by solving for the crossover point of
the parabolas. Hij is the interaction energy at the intersection of the two states. In the two
states shown, the state on the right has a lower overall configurational energy by ∆Go , but
the activation barrier must be overcome for charges to transport to this lower configuration.

kab = Aab exp −

∆Gab
kB T

,

(3.16)

where kab is the rate of a transition from state a to state b, Aab encompasses the probability of crossing from configuration a to b, and ∆Gab is the Gibb’s free energy difference
between the states. This exponential term can be thought of as expressing the probability of
forming a transitionary state between states a and b, which the charge may readily transition
through.
Solving for the activation energy given by the crossover point of the two parabolas in
figure 3.1, one finds that the free energy of activation is

(λ + ∆Go )2
∆G =
,
4λ

(3.17)

where λ is the reorganizational energy of the system and ∆Go is the total energy change

41

due to the reaction or hop. Typically, it is the reorganization energy which dominates the
activation energy, producing somewhat predictable results that hop rates to lower energy
sites will be faster than hop rates to a higher energy sites.
Note that an interesting consequence of Marcus theory is the prediction that hop rates to
configurations of much lower energy might be slower than the rate to a slightly higher energy
state (due to the squaring of the numerator in the exponential). This result is a phenomena
referred to as a Marcus inverted region, and has been experimentally observed[96].
For the DMC method, a simplified form of the Marcus theory is used to derive charge
hop rates. The prefactor term of equation 3.16 is set as a hopping rate constant, νo , the
exact choice of which is further discussed in appendix A.3. Incorporating the specific site
energies into the Gibb’s free energy term, the final form of the hopping rate is found

− Ej − Ei + ER
νij = νo · exp −
4ER kB T

2

,

(3.18)

where νij the hop rate from site i to j, νo is a prefactor, Ei/j are the site energies involved
in the transition, the difference of which generates the Gibb’s free energy of reaction, and
ER is the reconfigurational energy. This simplification of the prefactor is used to set the
hop rate to be consistent with the desired mobility for an isoenergetic system, set by the
Einstein relation.

3.4

Simple Bilayer Morphology Example

To demonstrate the effectiveness of the DMC model, we here present the case of a simple
bilayer[97]. Even though modern devices typically use a BHJ morphology[98], the bilayer
system allows us to explore the DMC model independent of complex morphological consid42

erations. We model the bilayer as two equally sized slabs of acceptor and donor material,
each 10 nm thick, that have periodic boundary conditions in the x and y dimensions with
contacts assumed on either side of the device in the z dimension, such that final system size is
35 nm × 35 nm × 20 nm. Excitons are generated uniformly throughout the morphology, with
a rate constant derived from the AM 1.5 light spectrum[82], scaled to the total volume being
simulated. Note that this example morphology is made intentionally thin, such that more
excitons will discover the interface and dissociate prior to decay and is not representative of
typical devices. The values of the various rates used in the simulation are given in table 3.1,
and are consistent with the values frequently used in DMC calculations of PCBM/P3HT
devices[47, 27, 2].

3.4.1

The effects of energetic disorder on a bilayer system

In order to demonstrate how the DMC model may be used to explore key features in OPV
devices, we first demonstrate the effect of increasing energetic disorder on the system. We test
this by widening the Gaussian distribution, σr , used to assign initial random site energies,
ER . Excitons are generated uniformly through the P3HT (donor) sites, and no dark charge
injection is modeled for this example. A difference in work functions of 1.2 eV is assumed
to exist across the morphology (Vbi ) which acts to drive charges toward their respective
contacts. We compare 4 values of increasing site disorder, σr = 0, 1/16, 1/8, and, 3/16 eV
at 21 different applied voltage values between 0 and 2 eV . For each combination of σr and
Vapplied , 20 simulations were averaged for each point simulated over 10 ms of illumination,
which was found to be sufficient to generate reliable statistics in the equilibrium region of
the device response. The results of these calculations are shown in figure 3.2, along with
error bars representing the standard deviation of each point.
43

(a)

(b)
Figure 3.2 Example of how energetic roughness would effect a bilayer. Data presented as
collection efficiency versus applied voltage (Vapplied ), demonstrating the well known diodelike behavior. Note that no dark current is being simulated, so no runaway current in the
opposite direction is possible. Figure 3.2a is the charge collection efficiency, χc-collect , while
figure 3.2b is the exciton dissociation efficiency, χx-diss .
44

Table 3.1 Parameters used in the bilayer example, except when otherwise specified.

Parameter

Value

exciton generation rate
exciton hop rate
exciton decay rate
exciton diss rate
charge hop rate
charge recombination rate
charge collection at electrode rate
φ0
EB (contact injection barrier)
a (lattice constant)
T (temperature)
ER (reconfigurational energy)
σr (width of random site energies)
Lx/y/z (lattice dimensions)

5 × 1027 m−3 s−1
0.02 ps−1
0.002 ps−1
100 ps−1
.676 ps−1
5 × 10−7 ps−1
.00646 ps−1
1.2 eV
.4 eV
1 nm
298K
.25 eV
.062415 eV
35 / 35 / 20 nm

Note the emergence of a dramatic turn-off voltage occurring between 1 and 1.5 V (depending on σr ), at which point charges are no longer capable of efficiently escaping the
interface and are lost due to recombination. As the energetic roughness is increased, the
voltage required to effectively pin the charges to the interface drops, such that they are more
likely to recombine. This result is not surprising, as the charge hop rates, given by equation
3.18, have a clear dependence on the difference of site-to-site energies. As the disorder in
random site energies was increased, so was the average difference between adjacent site energy values and thus the hop times between them. Longer average hop times meant charge
recombination became a more common event candidate and thus, occurred more often. Additionally, energetic traps formed by random low energy wells become more prevalent, which
act to inhibit efficient charge transport.
Exciton behavior for this bilayer system is shown in figure 3.2, demonstrating that while

45

applied voltage has no effect on the exciton behavior (as they are unaffected by the electric
field in this model), the dissociation efficiency dropped as the site disorder was increased.
Similar to the charge behavior, this is a consequence of decreased average hop rates, such
that exciton decay become more frequent.

3.4.2

The effects of sample depth on a bilayer system

To demonstrate that the DMC model can capture effects introduced by the morphology
itself, we look at a bilayer system in which the total length of the morphology, in the
z dimension, is gradually increased. We will look at bilayer systems with total vertical
dimension length Lz = 20, 40, 60, 80, and 200, with the heterojunction occurring in the xyplane at depth Lz /2 in each case. All other parameters of the simulation are shown in table
3.1. Once again, calculations are performed at 21 different values of Vapplied between 0 and
2 V , simulating 10 ms of illumination, with the data presented being the average of 20 full
simulations. Excitons are generated through the bulk of the donor sites and no dark charges
are simulated in this example. The results of these simulations are shown in figure 3.3.
As one can see, lengthening the z-dimension of the morphology (normal to the contacts)
has a dramatic effect on exciton collection. This illustrates the need in OPV for excitons to
be generated in close proximity to a HJ interface if they are to be successfully dissociated
into free charges. As the size of the morphology was increased, so was the average path
length from an acceptor site to the HJ, and thus excitons were capable of diffusing to the
interface prior to decay.
The charge collection efficiency, by contrast, had a much more subtle change across
different device lengths, as the data shows only a slight decrease in charge collection efficiency,
most notable near the cutoff voltage of approximately 1.5 V . This illustrates that so called
46

“geminate” charge recombination is the dominant loss mechanism in the DMC model, as
most charge pairs that escape the interface (and subsequently each other) are successfully
transported to their corresponding electrode, no matter what the path distance might be.
The minimum charge escape distance is dictated here by the Bjerrum length, as was discussed
in section 3.2.3. The effects of the Bjerrum length are will illustrated in this bilayer example,
as there are essentially no confinement effects influencing charge transport to the contacts
(other than the reflective wall of the interface).

3.5

Conclusions

In this chapter, we have introduced a characterization of the different scales of modeling
used in OPV simulations; these scales being atomistic (1 ˚- 5 nm), mesoscopic (2 - 400 nm),
A
and macroscopic (100 nm - 1 µm) models. We discussed previous calculations which were
performed with each of the model scales, and motivated our selection of a mesoscopic model
in order to study the effects of nanoscale morphology on OPV device performance.
We then introduced the dynamic Monte Carlo OPV device model, and explained it’s
ability to simulate exciton and charge transport behavior in detail. Derivations of a number
of the key mechanism rates were presented. Example calculations were presented on a simple bilayer morphology, demonstrating the effects of energetic disorder and layer thickness.
Clearly, the DMC model presents a tool to explore how both energetics of a system and
morphology would influence the internal transport properties of a device.
Deficiencies of the DMC model, as it was presented in this chapter, are computational
efficiency and the reliance on fitting parameters. Additionally, the demonstrated effects
of both model parameters (such as energetic roughness) and specific morphology (such as

47

(a)

0.4
L = 20
L = 40
L = 60
L = 80
L = 200

0.35

χx-diss

0.3
0.25
0.2
0.15
0.1
0.05
0
0

0.5

1

Vapp

1.5

2

(b)
Figure 3.3 DMC calculations showing the effect of varying total bilayer thickness on charge
collection efficiency, χc-collect (top), and exciton dissociation efficiency, χx-diss (bottom).
Presented are the results for systems of thickness L = 20, 40, 60, 80, and 200, with total
morphology size 35×35×L, and with periodic boundary conditions in the x and y dimensions.
No dark charge is simulated here. All other parameters used in the simulation are shown in
table 3.1.
48

the dimensions of model used) suggest that it’s crucial to use the most realistic model
morphologies possible in order to improve the reliability and use of the DMC simulation
results.

49

Chapter 4
A Novel Method to Simulate Small
Angle Scattering
Many materials under study for modern energy and industrial applications are composed
of percolative, interpenetrating, complex nanostructures, such as Nafion networks in fuel
cells[99, 100, 101], structural materials[102, 103], organic photovoltaics[16, 18], batteries[104,
105], OLED’s[106, 107], and sensor materials[108]. Though the exact mechanisms influencing device performance in many of these systems are not yet well understood, the transport
properties are known to be strongly correlated with the details of the complex nanoscale
material networks, and as such much recent work has involved structural characterization of
these morphologies. One such measurement is small angle scattering, which can inform on
nanometer scale internal structure. Of particular interest to the study of organic systems is
small angle neutron scattering (SANS), which often can provide contrast between material
components that are indistinguishable with small angle x-ray scattering (SAXS). The challenge of SANS and SAXS measurements lies in interpreting the data, as standard analysis
requires fitting experimental data to well defined analytical form factors and structure factors. Although informative, these standard fitting procedures do not directly extract a 3-D
model, which is necessary for device simulation. Reverse Monte Carlo methods are capable
of handling small scale structures (on the order of 5-10 nm) but larger structures required for

50

full scale device modeling are computationally difficult to refine using conventional methods.
To extract large scale nanostructural models from SAS data we developed a novel algorithm to efficiently calculate the theoretical scattering structure function for morphological
models of arbitrary scale and resolution that we refer to as the distribution function method
(DFM). We will present the DFM algorithm here, discuss the limits of it’s applicability,
and demonstrate it’s usefulness through a range of simple to complex morphological model
examples.

4.1

Principles of the Fourier Transform

Before delving into the principles of scattering, it is worthwhile to discuss the principles of
the Fourier Transform as they relate to scattering. The Fourier transform (F.T. for short)
is a decomposition of a function into the sum of it’s complex periodic frequencies. One
commonly defined form of a FT is
∞
1
√
F (k) =
f (x) e−ikx dx.
2π −∞

(4.1)

To give an example, consider a simple cosine with frequency ν, defined in real-space as
f (x) = cos(νx). Recalling that cosines may be written as cos(x) = (eix + e−ix )/2, the
Fourier transform into wave number space, k, is
∞
1
F (k) = √
cos(νx) e−ikx dx
2π −∞
∞
1
= √
e−ix(k−ν) + e−ix(k+ν) dx
2 2π −∞
π
δ
+ δ(k−ν)
=
2 (k+ν)

51

(4.2)

where we have used the definition of the Dirac Delta function

∞

2πδ(k−ko ) =

−∞

e−ix(k−ko ) dx =



 1

 0

if k = ko

(4.3)

otherwise.

The result expresses the function, f (x) = cos(νx), in wave number space, F (k) =
π
2

δ(k+ν) + δ(k−ν) , which demonstrates that this function is composed only of periodic

frequencies with wavenumbers k = ν and k = −ν. It’s important to note that f (x) and F (k)
carry identical information, only in different forms, and that performing a Fourier transform
from one form to the other will not cause any inherent loss of information. Also, note the
choice of Fourier function used in this example is not unique, as many different forms are
acceptable so long as they conform to principles of a Fourier Transform. Here, we used a
form that generated functions in terms of wave number, k, from a real space coordinate, x,
but could have just as easily used a form which converted functions in time, t, to functions
in frequency space, f . Wave number is a natural corollary to frequency, as it expresses the
number of waves which may fit in a distance of 2 pi (units of cycles per unit length), being
defined as k = 2π/λ, where λ is wavelength.
A more complex FT example is given in figure 4.1, which demonstrates how a complex
function (in this case, a series of broad delta-functions) would generate a series of peaks
which are symmetric about zero. As is explained in detail in appendix A.2, if the real
space function represented an aperture (series of slits), the Fourier transform would be the
transmitted intensity of an incident wave through the aperture.
Some important principles of Fourier transforms are that first, they work through decomposition of periodic functions into their symmetric (cosine) and asymmetric (sine) components. Though this might at first seem a narrow band of functions (periodic only), realize
52

(a)

(b)

(c)

Figure 4.1 A demonstration of the Fourier relationship. Figure 4.1a shows an arbitrary
function composed of superimposed Gaussian distributions of varying heights. Figure 4.1b
demonstrates the real component of the Fourier transform of the function in figure 4.1a.
Figure 4.1c demonstrates the auto-correlation of the function in figure 4.1a.
as the length of the periodic cell grows to infinity, no repetitions become necessary and thus
any function may in principle be Fourier transformed. Second, Fourier transforms of real
and symmetric functions produce a function which is real and even, while Fourier transforms
of functions which are real and antisymmetric produce a function which is imaginary and
odd. Finally, the value of the Fourier transform at the origin, F (k = 0) is proportional to
the total area under the curve of the original function, f (x).
A useful mathematical concept in Fourier transforms and their applications is the convolution theorem, which states that two functions, g(x) and h(x), may be convoluted by one
another (defined with the symbol ⊗) as

g(x) ⊗ h(x) =

∞
−∞

g(t)h(x − t) dt.

(4.4)

The effect of a convolution is to “blur out” the function g(x) and h(x) by each other, as the
two functions are essentially dragged across one another, causing the resultant convolution
to be significantly different than a simple product of the two functions. An example is given
in figure 4.2, showing a series of delta function peaks, of various amplitudes, convoluted with
a Gaussian distribution.
53

(a)

(b)

(c)

Figure 4.2 Convolutions of Gaussian functions of varying widths with a series of delta
function peaks. The function describing the delta functions is f (x) = δ(x + 10) + 2δ(x + 5) +
3δ(x) + 2δ(x − 5) + δ(x − 10). The Gaussian is described as g(x) = exp −x2 /2σ 2 . Widths
of Gaussian distributions used in these three examples are σ 2 = .1 in figure 4.2a, σ 2 = .5 in
figure 4.2b, and σ 2 = 3 in figure 4.2c.
One convenient aspect of convolutions relates to their effect on the Fourier transforms of
resultant functions, giving the relation

f (x) = g(x) ⊗ h(x) ⇔ F (k) =

√

2π G(k) × H(k).

(4.5)

Notice that functions g(x) and h(x) have been convoluted (⊗) to produce f (x), but the
Fourier transform of this new function, F (k) is the result of the product of the Fourier
transforms of the two initial functions, G(k) and H(k). This result is useful, as many times
taking the convolution of two systems can be complex, where as the Fourier transform and
subsequent product is simple.

4.2

Basics of scattering

Scattering measurements are characterizations of how an incident particle, say an x-ray or a
neutron, is effected by a sample. We describe the incident particle through it’s energy and
momentum, with the basic equations

54

λ=

E = hν =

h
p

(4.6)

hc
= eV = ω = c|k|,
λ

where λ is particle wavelength, h is Planck’s constant,

(4.7)

is the reduced Planck’s constant

( = h/2π), p is a particle’s momentum, E is it’s energy, ν is it’s frequency, and c is it’s
velocity (the speed of light). It is common to express the energy of a wave in units of
electron volts (eV ), a measurement of the energy required to take a single electron of charge
e ≈ −1.6 × 10−19 C across a potential of 1 volt.
The incident particle’s momentum can be expressed in terms of wavenumber, k, as

k=

2π
2π
=
p.
λ
h

(4.8)

Wave number describes the number of times a wave oscillates in a given length, which is
directly proportional to momentum through a constant.
In the case of neutrons, their kinematic energy may also be described via the wavenumber
as

E=

2k2
p2
h2
=
=
,
2m
2m
2mλ2

(4.9)

where m is the mass of the neutron. Scattering particles are described through their
waveform, ψ. The equation of a propagating plane wave is written

i k·r−ωt

ψ = Ae

55

(4.10)

A = |A|eiφo ,

(4.11)

where A is the complex amplitude of the wave which has absolute amplitude, |A|, and
phase offset, φo , with frequency ω at time t.
When a particle scatters from a sample, the interaction may be characterized by the
change in momentum and energy. The momentum change, P is

P = ki − kf = q,
where

(4.12)

= h/2π is the reduced Planck’s constant, ki/f is the initial and final momentum

of the particle, and q is the transferred momentum. The energy transfer from the incident
particle to the sample, ∆E, is

∆E =

ωi − ωf .

(4.13)

When no energy is transferred (∆E = 0), the scattering is referred to as elastic. All
scattering events discussed in this thesis will be in the elastic regime. The implications of
elastic scattering and equation 4.7 is that the wavelength, λ, cannot change before and after
a scattering event, and by extension neither can the magnitude of the wavevector, |k|. This
condition leads to the result

|kf | = |ki | =

2π
.
λ

(4.14)

Thus we have introduced momentum transfer, q, which is a common variable in many
scattering measurements.

56

4.2.1

Single Scattering Event

An example of a simple elastic scattering event is shown in figure 4.3, where an incident
wave with momentum ki is scattering off a sample at an angle 2θ, such that it has final
momentum kf . Described using momentum transfer, q = kf − ki , we see that the amplitude
of q is

q=

4π sin θ
.
λ

(4.15)

Figure 4.3 Single scatter event example
We may describe the spherical scattered wave, ψf , as

ψf = ψo f (λ, θ)

eik·r
,
r

(4.16)

where ψo describes the incident plane wave, f (λ, θ) is a function describing the fraction
of an incident wave of wavelength λ that is deflected in the direction, θ. Note that there is
no azimuthal (φ) dependence on the function f (λ, θ) due to the spherical symmetry. The

57

scattering strength of a particle, f (λ, θ), explicitly describes how strongly the incident wave
interacts with the scatter. Note that the term “scatterer” is used here, as although the
scattering is principally caused by interactions with atoms, it is sometimes more convenient
to consider average scattering from a molecule or small volume of material. The scattering
interaction itself is caused by a number of different specific mechanisms, dependent on incident wave and sample. For instance, electrons may scatter due to their Coulomb interaction
with an atom’s electric field, hence larger atoms (with greater Z number) scatter x-rays
stronger than smaller atoms. Neglecting magnetic interactions, neutrons interact only with
the nucleus of the atom and thus their scattering strength has no simple correlation with
atomic size. As these scattering interactions can be dependent on incident particle wavelength and angle, we have up until now included their explicit dependence in f (λ, θ). As
we shall discuss later, in the q-range used in many small angle scattering measurements the
strengths of these interactions can often be approximated as a constant.

4.2.2

Multiple Atom Scattering

Scattering from an ensemble of atoms can be calculated by applying the principles of superposition to the methods introduced for scattering due to a single atom. Assuming the total
scattered wave, ψf , will be composed of the sum of the scattering due to each individual
atom, δψf j , where the index, j, identifies the scattering atom, one may rewrite equation
4.16 as
ikf ·(r−rj )

e
ik ·r
δψf j = ψo e i j fj (λ, θ)
,
|r − rj |

(4.17)

where the j-th atom is at coordinate, rj . Taking the sum over all N atoms in this

58

ensemble, with each atom located at rj we find

N

ψf =

δψf j
j=1
N

=

ik ·(r−rj )
e f
fj (λ, θ)
|r − rj |

iki ·rj

ψo e
j=1

= ψo e

i kf −ki ·rj

N

ikf ·r

fj (λ, θ)
j=1

e

|r − rj |

,

where as we have rearranged the terms, it’s straightforward to rewrite the exponential
in the sum in terms of momentum transfer, q. Additionally, we may assume that the interatomic distances, |rj − rk |, in the system is much smaller than the distance to the detectors,
|r|, and thus the far-field approximation can be made, |r − rj | ≈ |r| = r. Taking this limit
into account, we may write

N

ikf ·r

ψf = ψo e

iq·rj

fj (λ, θ)
j=1

e

|r|

.

(4.18)

Finally, as we are interested in the net intensity (I ∝ |ψf |2 ) measured at the detectors,
we take the modulus squared of ψf to find

2
ψo
|ψf |2 = 2
r

N

2
iq·r
fj (λ, θ) e j

.

(4.19)

j=1

As we shall see later, this result is of particular use when we derive the Debye formalism
in section 4.3.2.

59

4.2.3

Scattering strength

The principle mechanism by which matter scatters an incident wave is dependent on the
incident scattering particle. In the case of x-rays and electrons, scattering is primarily
dominated by electrostatic interactions with the electrons surrounding an atom, hence the
stronger interaction with atoms of greater Z number. Electrons primarily scatter due to
the Coulomb interaction with the electron cloud, whereas x-rays primarily scatter due to
interaction with the electron cloud itself[109]. For this reason, x-rays may penetrate deeper
towards the nucleus than an electron. For neutrons, scattering may occur either due to the
nuclear interaction at the nucleus of the atom or via a dipole-dipole interaction referred to
as magnetic scattering[110].
In the classical image of scattering, we may imagine the incident wave as a single particle
and the scattering atom as a single spherical body of finite size. In this picture, the relative scattering strength of an atom may be represented by a larger overall scattering area.
Although this is a simple classical interpretation of scattering, it inspires the common unit
of relative scattering strength called the scattering cross-section, σ. Because of the strength
of the scattering interactions, electrons will encounter much larger overall scattering crosssections than x-rays, and neutrons will encounter the smallest cross-sections overall. For this
reason, neutrons are capable of probing much deeper into many materials than electrons or
x-rays, particularly those with high Z number such as metals.
The value of σ may be measured by looking at the rate of incident particle scattering
over all angles, R(s−1 ), over the total amount of incident flux on the scattering sample,
Φ(s−1 m−2 ). The scattering cross-section is related as σ = R/Φ, with units of m2 .
In principle, the scattering cross-section will be dependent on both the incident angle (θ)

60

and wavelength (λ) of the incident wave, as the driving scattering interaction mechanisms
are so dependent. However, in the regime that small angle scattering measurements are
taken the relative q-values of incident scatterers will be such that little variation occurs in
the σ(θ, λ)[111], and as such we may refer to a constant scattering cross-section, σ. The
scattering length, b, is another quantity which may describe material scattering strength.
Scattering length is related to the scattering cross-section, similar to the relationship between
the radius of a sphere and it’s surface area

σ = 4π|b|2 .

(4.20)

As the measurement taken in any small angle scattering experiment is the number of
incident scatterers deflected to a specific solid angle defined by the measurement detector, a
useful quantity to discuss is the differential cross-section

dσ
dΩ

=

R(2θ, φ)
,
N Φ∆Ω

(4.21)

where the number of scatterers deflected per solid angle, R/∆Ω, is normalized by the
total number of scattering particles incident on a sample, Φ, and the number of scattering
units (such as atoms) in a sample, N . All scattering due to internal structural features inside
a sample will be present in the differential cross-section[111, 110]. Although scattering from
a sample will always be composed of an elastic (∆E = 0) and inelastic (∆E = 0) component, the inelastic scattering will generally be small compared to the elastic component, and
constant across the q-ranges of interest here. As such, we treat total differential cross-section
and elastic differential cross-section as one-in-the-same for the purposes of this thesis.
The elastic scattering to a specific solid-angle due to a sample may be written
61

R(2θ, φ) = ψf

2

δA,

(4.22)

where δA is the area subtended by the solid angle defined by (2θ, φ).
Along with the definition of scattered wave due to an ensemble of scatterers provided by
equation 4.19, we may describe the differential cross-section as

dσ
dΩ

1
=
N

2

N

bi eiq·ri

,

(4.23)

i=1

where the sum is taken over all N scatterers of the system, each of which has scattering
length density, b, and exists at vector ri in the system. Note that here, the differential crosssection has been normalized by total number of scatterers, N , but the exact normalization
may vary between disciplines to include such values as mass, volume, average scatterer
strength, and so forth[110]. The overall shape of the differential cross-section in all cases is
always a function of the particular arrangement of scatterers in a given sample, as well as
the scattering strength of each of these components. This is why all structural information
of a sample is found in the differential cross-section.

4.2.4

The Full SAS Picture

Up to now, we have discussed scattering due to individual scatterers, each described by
a scattering length, b. However, we may instead treat scattering from a large number of
atoms in a finite region with a scattering length density, β(r). The scattering length density
(SLD) may be defined by assuming that scattering from a point source of constant scattering
strength, b, is described as

62

b = β(r)dV,

(4.24)

where the vector, r, defines the coordinate of the total scattering volume dV . Using β(r)
and equation 4.23, we may generalize the differential cross section as
2

dσ
∝
dΩ el

β(r)eiq·r d3 r .

(4.25)

V

Notice that by defining the differential cross-section in this way, we see that the structure
of a sample (explicitly defined by the function β(r)) is related to the differential cross-section
by a Fourier transform.
The absolute measured intensity due to elastic scattering in a SAS measurement may be
written

I(q) = Io (λ)η(λ)T (λ)∆Ω

dσ
⊗ Rinst (q) + B(q),
dΩ

(4.26)

where the incident flux, Io , efficiency of the detector, η, and transmission through the
sample, T , are all dependent on incident wavelength. ∆Ω represents the angle subtended by
the detector at momentum transfer q. All these factors are convoluted with the instrument
resolution, Rinst , except for the background signal, B.

4.2.5

Form Factor vs. Structure Factor

Often times, small angle scattering data is discussed in terms of a form factor, P (q), and a
structure factor, S(q). In principle, a form factor (or atomic form factor) is the scattering
profile due to the scattering from an individual scatterer, such as an atom or a spherical

63

nanoparticle. If a sample is composed of an ensemble of scatterers, the arrangement itself
contributes to the scattering profile as a structure factor. The overall scattering profile will
be the product of the form factor and structure factor

I(q) = P (q) × S(q).

(4.27)

Because of the Fourier relationship between real space structural information and qspace scattering data, scattering from small physical features (small r information) will be
represented inversely by large q-scale features and vice verse. Thus, if a system is composed
of well dispersed scatterers in medium transparent to the incident wave (the dilute limit),
then any contributions due to long range ordering of individual scatterers will occur in the
low-q region of the data, well separated from the form factor contributions of the individual
scatterers.

4.3

Scattering Simulation

In this section, we introduce two general methods of calculating the scattering profiles of
model morphologies, which we refer to as the direct method and the Debye method. Both
methods rely on the inherent Fourier relationship between real space data and scattering
data. We will describe both these methods in detail, and discuss the strengths and weaknesses
of each. While other forms of structural scattering simulations exist[112, 113, 114, 115, 116,
117], we limit our comparison here to the direct and Debye formalisms because of their
underlying mathematical symmetry.
Both the direct method and the Debye formalism require discretization of a scattering

64

Figure 4.4 Comparison of methods used to simulate small angle scattering, highlighting the
overall mathematical similarities of the two methods.
structure into essentially point scatters and rely on the Fourier relationship between a morphologies’ real space function and it’s scattering profile. The differences between the two
methods are qualitatively illustrated in figure 4.4. Note that both of these methods assume
the scattering data involved is angle averaged[111].

4.3.1

Direct Method

The direct method of simulating small angle scattering begins by defining a morphological
model on a discrete lattice and then taking a full 3-d Fourier transform of the structure,
before finally projecting the resultant 3-d complex lattice down into one dimension[26].
In detail, the values of the typically cubic lattice sites on the real space model correspond
to the scattering length densities of the associated material. This results in a discretized lattice of SLD values, ρ(x, y, z). The full 3-D Fourier transform is then calculated on ρ(x, y, z),
which results in a a 3-D lattice of scattering intensities, I(qx , qy , qz ). In order to compare

65

directly to experiment, it is necessary to project this scattering information down into one
dimension, resulting in I(q). This projection involves first summing cylindrical shells of
I(qx , qy , qz ), centered about the z-axis (considered parallel to the direction of the detector
from the sample), which results in a 2-D rings of scattering intensity. These rings are integrated, much as they are in experiment, to generate the final angle averaged 1-D scattering
intensity plot, I(q). Finally, in order to account for the natural resolution limits of experimental measurements, a smearing function is often applied to I(q). This smearing may also
be used to mask the effects of the lattice spacing.
In order to alleviate the effects of discretizing a morphology into point scatterers, one can
convolute a sinc function in 3 dimensions with the 3-D scattering data prior to projection.
The reasoning for this is that the Fourier transform of a sinc function is a square wave,
and considering the relationship described in equation 4.5, this is equivalent to calculating
the scattering from a structure generated by point scatters multiplied by cubic scatterers,
thus filling the void space between the point scatterers so long as the width of the sinc
function is selected accordingly[26]. The downside to this approach is that, while one does
effectively eliminate void space, one is left with a structure of sharply transitioning nonphysical densities which thus must be later smoothed out with a smearing function regardless.
The only information added to the structure is unphysical, and thus we feel this step is
superfluous.
The advantage of the direct method is that it is a straightforward calculation which scales
as N ln N , which is relatively efficient for small systems. However, the need to perform a
3D Fourier transform and convert these results into a 1D form can become computationally
intensive for large systems.

66

4.3.2

Debye Method

The Debye method involves construction of a function containing structural correlation information in real space, followed by Fourier transforming this function to generate scattering
information in reciprocal space. The derivation of the Debye formalism begins by defining
the scattered amplitude due to an ensemble of atoms

A(q) =

bk e −i q·rk ,

Ak = C
k

(4.28)

k

where C is the scattering amplitude of a single atom. The measured scattering intensity
from the ensemble of atoms, I(q) = A∗A, may thus be written

I(q) = A∗A
= C2

bj bk e
j

e i q·rk

k

= C2

bj bk e
j

−i q·rj

−i q· rk −rj
.

k

(4.29)

Every segment in this arrangement of atoms, (ri − rj ), will have a corresponding reverse
segment of equal magnitude but opposite direction, (rj − ri ), we may combine both these
contributions as follows

67

−ix· a−b

e−ix·(a−b) =

1
2

e

=

1
2

e

−ix· a−b

+e
+e

−ix· b−a

ix· a−b

= cos (x · rab ) = cos (xrab cos φ) ,

where rab = a − b, and φ is the angle between a and b, φ =

(4.30)

a,b . Using the results of

equation 4.30 in equation 4.29, we may write the intensity as

I(q) = C 2

bj bk cos q · rjk
j

k

= C2

bj bk cos qrjk cos φjk ,
j

(4.31)

k

where φjk here defines the angle between the incident scattering wave, q, and the pairposition vectors, rjk , where φ ∈ [0, π]. Under the assumption that scattering from all angles
are equally likely, in the words of Guinier, “the mathematical problem” is in finding the
average value: cos q · rjk .
To calculate this average, we begin by stating that if all angles, φjk , are equally likely,
1
then the probability that an angle exists between φ and φ + dφ is given by: 2 sin (φ) dφ. The

average of the phase function, cos q · rjk , can thus be written

68

cos q · rjk

π

=
0

cos q rjk cos φ

sin φ
dφ
2

π/2

=

cos q rjk cos φ sin φ dφ

0

π/2
1
cos q rjk cos φ d q rjk cos φ
q rjk 0
0
q
=−
cos (u) du
q rjk hr

=−

=

1
sin q rjk
q rjk

,

(4.32)

where u = h rjk cos φ . Simplifying variables such that rjk = r, this produces the result

cos (q · r) =

sin (qr)
,
qr

(4.33)

such that when utilized in the angle averaged form of equation 4.29, we generate the well
known Debye formula for scattering intensity

I(q) =

bj bk
j

k

sin(qrjk )
.
qrjk

(4.34)

While a straightforward calculation of equation 4.34 is acceptable for small scale systems,
it becomes computationally difficult as the number of scatterers increases. This can be
somewhat circumvented by first calculating a weighted pair density, ζ(r), defined as

bj bk δ(r − rjk ),

ζ(r) =
j

(4.35)

k

such that a sine-transform (often referred to in the scattering community simply as
the Fourier transform) of ζ(r) needs only to be calculated once (for each value of r seen),
69

generating the scattering intensity

∞

I(q) =

dr ζ(r)
0

sin(qr)
.
qr

(4.36)

Note that the Fourier relationship between scattering data and the real space structural
data, as defined by ζ(r) here, may be written in a variety of forms. Any function which
contains real-space structural pair data, such as the atomic pair distribution function, radial distribution function, etc., can be Fourier transformed into scattering data, with the
differences being only in normalization and scaling. We compare a number of different relationships commonly used in the literature in appendix A.4 of this thesis.
As compared to the direct method, the Debye method requires the additional step of
calculating pair-structural information (ζ(r)), however this step circumvents the need to
later project the scattering data down into one-dimension, as ζ(r) is one-dimensional. It
is the calculation of ζ(r), requiring a sum over N (N − 1)/2 pair distances, which is the
dominant computational step. Thus, the Debye method scale as N 2 , making them less
efficient than the direct method that scaled as N ln N . However, we found neither method
to be acceptably efficient for use in reverse Monte Carlo-like structural refinement of large
scale morphologies.

4.4

DFM Algorithm

We have found that for many systems, the full calculation of the structural correlations
is not necessary to estimate the Debye form factor to high precision, as one may instead
randomly sample the structure to generate a weighted distribution function, w(r), which
acts as the real space structural data (defined by ζ(r) previously) for the morphology. We
70

thus refer to our algorithm as the distribution function method (DFM). The number of pairs
sampled per calculation may be scaled according to the demands of the application and
morphology itself, with higher sampling rates corresponding to greater accuracy and lower
rates with faster calculation time. Once generated, a sine-transform of w(r) generates an
approximated scattering profile, I a (q). This allows the calculation to scale as f × N 2 , where
f is the ratio of random pairs sampled, Np , over the total number of pairs possible in the
system such that 0 < f < 1. Thus, the DFM approach can be adjusted such that it is more
efficient than the direct method when f < ln N/N by sacrificing some level of accuracy.
For our systems, we find that a w(r) which has been calculated with as little as 0.001% of
the total pairs to be found in the system will generate an accurate scattering profile. Though
further sampling refines high-q details and may increase confidence in the simulated scattering structure, for purposes of empirical reverse Monte Carlo-like refinement, the abbreviated
sampling works remarkably well.
We will begin by detailing the algorithm and sampling method used in the DFM algorithm, demonstrating it’s validity through a series of simple examples. We then examine
the limitations of the model as well as methods to improve results as demonstrated with
examples of more complex morphologies, such as systems of monodisperse and polydisperse
spheres and an interpenetrating two-phase structure typical of the nanoscale morphologies
found in many energy applications.

4.4.1

Methodology

The DFM algorithm assumes a morphology defined by a discrete set of points, with each
site assigned a scattering strength based upon the scattering length density (SLD) of the
represented material. A sites which is considered a mix of different materials should be
71

assigned an appropriately weighted SLD value. The approximated scattered intensity, I a (q),
is thus

Rmax

I a (q) = C

w(r)
r=Rmin

sin qrij
,
qrij

(4.37)

where the constant, C, is typically adjusted to compare the experimental system of
interest. For direct comparison to the full Debye calculation (over all pairs), one can set
C = (N 2 − N )/Ns , where N is the total number of scatterers and Ns is the number of pairs
sampled. The sum is taken over all binned distances in the system, Rmin to Rmax , and we
have incorporated the weighted distribution function, w(r). Here we define w(r) as

Ns

bi bj δ(r − rij ),

w (r) =

(4.38)

pairs
where the sum is taken over Ns sampled pairs of sites, i and j, which are distance rij
apart with corresponding scattering weights bi and bj . Typical material SLD values (for
−2

neutrons) are −1 × 10−6 ˚
A

to 5 × 10−6 ˚
A

−2

. To simplify the sum of w(r), if the SLD values

used in the model have integer ratios to each other it is convenient to reduce the weights
to integer values for computational efficiency. The difference between using the real values
and representative integer values can be later adjusted in the normalization. The weighted
distribution function, w(r), is defined similarly to the atomic distribution function, ρ(r),
with differences in the normalization (see appendix A.4).
We find that C × I a (q) calculated from w(r) quickly converges to I(q) (the results of the
full calculation) for many systems, particularly those with strong contrast between different
scattering species such as the 1-to-0 scattering values assigned in the examples here. As
in our simulation w(r) is not a continuous distribution but rather a histogram, care must
72

be taken to select an appropriate number of bins, Nbins = D/dr where D is the longest
length measured and dr is the width of a single bin. Note that in a cubic system of width
√
√
L, the longest possible distance is 3L in a non-periodic system, and 3L/2 in a periodic
system. Sufficient binning of the sampled data is necessary to avoid violating the Nyquist
limit of dr ≤ π/qmax and related aliasing effects. By contrast, a much finer grid than that
dictated by the Nyquist limit is unnecessary, as information on distances smaller than dr is
not representative of any actual morphological features. Though acceptable if smoother data
is desired, over binning will slow down the Fourier transform and thus the overall algorithm
efficiency. Additionally, when over-binning care must be taken not to misinterpret features
smaller than the Nyquist limit as these are unphysical. We typically use 10,000 bins in
producing w(r) in equation 4.38, and 30,000 when transforming to I a (q) in equation 4.37,
as we find this number is sufficient for the size L = 400 systems we typically model.
While it’s possible to use different sampling methods in order to weight different regions of the morphology preferentially, here we sample all pairs with equal probability. A
demonstration of the accuracy of the method is shown in figures 4.5 and 4.6, which compare
the analytical form factor and DFM calculated structure for a sphere and set of cylinders
respectively, where the structures have been modeled with coordinates on a cubic lattice.
The analytical form factors shown in these two examples were generated using Irena Igor
macros[118]. As can be seen, the form factor of these systems is recovered to a high precision by sampling a total of 107 total pairs, as compared to the total number of pairs in the
models, being approximately 1.7 × 1010 for the sphere in figure 4.5, 1.9 × 1010 for the smaller
cylinder and 3.1 × 1011 for the larger cylinder in figure 4.6.
In each of these examples the constant C is selected such that the simulated scattering
intensity, I(q), is normalized to unity as q → 0. Note also that only in the case of the large
73

Figure 4.5 Comparison of scattering from a single R = 100 ˚ sphere calculated using the
A
analytical form factor and the DFM algorithm. Analytical form factor provided by the
IRENA SAS Macros in Igor, where sphere to solvent contrast is taken 1-to-0. DFM calculation performed over a randomly sampled set of 107 pairs, or .05% of the total model pairs.
The lattice spacing here is taken as al = 1 ˚.
A

74

˚
Figure 4.6 Simulated scattering from cylinders of length, L = 100 A, and radius, R =
˚, calculated using the analytical form factor for a cylinder and compared to the
25, 50 A
scattering calculated with the DFM algorithm. Analytical form factors provided by IRENA
SAS macros in Igor, where cylinder to solvent contrast is taken as 1-to-0. DFM calculation
performed over a randomly sampled set of 107 pairs, or .005% of the total system pairs.

75

cylinder, R = 50, L = 100, constructed from Np = 785524 discrete points, do we see any
deviation from the analytical from factor. The deviation we do see corresponds to small
scale features which are not as well represented due to the number of pairs sampled. We
find that when the number of pairs sampled is set equivalent to the total number of pairs in
the system, there is no discernible difference in the resultant simulated scattering implying
that the sampling method itself does not introduce any additional structure-like features in
the scattering profile.
As small angle scattering data is useful for a number of polymeric system studies such as
observing phase transitions[119, 120], we next look at an example of scattering from an ideal
polymer chain. For this calculation, we simulate a polymer with one million monomers,
generated through an unbiased (non-self avoiding) random walk in real-space, with bond
length of 1˚. For each calculation, 1000 polymers are generated and then each sampled to
A
generate w(r), and the results are shown in figure 4.7. The measured radius of gyration for
the simulated polymers was found to be rg = 398 ˚, and we compare the calculated scattering
¯
A
to the Debye form factor[121] for this value of rg . To explore the effects of sampling rate,
scattering profiles were generated using both 106 pairs and 100 pairs per configuration.
We see that the results from the DFM algorithm match the Debye form factor to a high
degree in both sampling rates, though the errors introduced by the lower sampling rate are
evident at high-q. Note that this result is of particular interest, for while the higher sampling
rate (.0002% total pairs sampled) captures the analytical solution exactly, the low sampling
rate manages to captures all critical information up to the high-q limits of the model where
it begins to deviate. This result demonstrates the usefulness of the DFM algorithm as even
though the low sampling rate is looking at 10 orders of magnitude less total pairs than the
full Debye formalism, a nearly identical scattering profile is generated.
76

Figure 4.7 Simulated scattering from an ideal chain polymer model with 106 monomers with
resultant radius of gyration, rg ≈ 398 ˚. The DFM calculation is performed by sampling the
A
given number of random pairs on 1000 random polymer configurations. Note that for the
size of polymer simulated here, the total calculation would require a sum of nearly 5 × 1011
total pairs per polymer configuration. Even with only 100 pairs sampled per polymer, the
analytical solution is exactly recovered up to high-q values.

77

4.4.2

Limits of method

Though the examples shown in figures 4.5 - 4.7 are presented as evidence of the validity of
the DFM algorithm, it’s primary use is in the efficient calculation of scattering profiles for
complex structures which are not well described via analytical formalism. However, it is
important to understand the natural limits of this method if the profiles generated are to
be considered reliable. The two primary limitations on any DFM algorithm are finite size
effects and sampling rate errors.

4.4.2.1

Finite size effects

Finite size effects are seen as comparatively high-frequency oscillations in q−space, where the
frequencies correspond to the length scales of the sample box. In cases where the relevant
structural information is contained at higher q-values, which are well separated from the
low-q finite size effects, interpretation of the data may seem straightforward. However, the
oscillations induced by the finite size effects may still interfere with correct data analysis due
to the hard-wall like nature inherent in the abrupt termination of the structural information
at the box edge. To alleviate this effect, we suggest a Gaussian “fuzzying” routine be
performed on the morphological models prior to scattering profile calculation, which acts to
roughen the box edges. Every site in a morphology is assigned a random number, X ∈ [0, 1],
which is then compared to the value of a Gaussian distribution, G(σ), with origin at the
center of the morphology with a width of σ. If for a given site, if X > G(σ), the site is pruned
from the list of sites to be used in generating w(r). Thus, the surfaces of a morphological
model are roughened such that the contributions to the scattering profile will be damped as
compared to a hard-wall sample box.
An example of the Gaussian “fuzzying” routine is given in figure 4.8, which demonstrates
78

Figure 4.8 DFM calculation of scattering from a system of 697 spheres of radius r =
3
20 ˚ placed randomly inside a box of size 8003 ˚ which were selected to have a minimum
A
A
distance between sphere centers of L = 100 ˚. Demonstrated are the effect of calculating the
A
scattering on this morphology while points were randomly removed from the structure using
a “fuzzying” routine via a Gaussian probability with width σ. This procedure reduces the
inherent hard-edged box effects seen as high frequency oscillations at the expense of total
signal contrast and sample size.

79

the routine for a system of monodisperse spheres of radius r = 20 ˚ which have been ranA
domly placed at a minimum center-to-center distance of D = 100 ˚ from one another in a
A
cubic sample box of size L = 800 ˚. Presented in figure 4.8 are the full data (no “fuzzying”
A
routine performed) for comparison with data from σ = 30, 60, 120, and 240 ˚. The Guinier
A
regime of the plot, seen as a plateau at q → 0, will have an intensity dependent on the
total number of scatterers. We thus expect higher intensities for systems which have larger
values of σ, as less points will have been removed from the structural model. For purposes of
comparison, we normalize all the presented data such that I(0) = 1, leading to the apparent
rise in intensity of the scattering peaks as σ is decreased. Note also that as the width of the
Gaussian is decreased, the shrinking overall model size is evident in the shift of the left most
bend in the scattering profile, which is increasing in q-value as the system size is shrinking
in real space. In large box simulations, it is important that the box size contributions to the
scattering are kept well separated from the sample features of interest.
By applying an appropriate σ to the fuzzying routine, one is able to dampen the finite size
effects at the expense of contrast between structures found at long range in the sample and
overall model size. For example, in figure 4.8, one might select 120 < σ < 240 ˚ (between
A
15% to 30% the total box width) in order to better resolve both the average center-to-center
−1

structure factor peak seen at q ≈ (2π/100) ˚
A
−1

at q > .2 ˚
A

as well as the spherical form factor seen

. Smaller σ values than this, correlating with more sites pruned, do little to

further clarify the data and may actually begin to obscure the information of interest. This
trade-off between total information sampled and clarity of the signal must be evaluated for
each morphology.
If using a lattice model, the lattice constant employed (a) will contribute a strong signal at q = 2π/a, which typically occurs at much higher q-values than the features one is
80

investigating. Undesirable contributions of the lattice constant will rarely be a problem as
long as the lattice resolution employed by a model is sufficient to capture morphological
characteristics of interest. Peaks due to the bin size used in generation of w(r) or I(q) will
also occur, but as long as the bin width is less than the lattice constant, these will not occur
until much higher-q.

4.4.3

Sampling rate errors

When using the DFM algorithm, choice of sampling rate is critical in the generation of a
w(r) that will lead to an accurate I(q). The ideal rate to use, however, is highly dependent
on individual morphology. If the structure is isolated and a simple fit is all that’s required,
such as was the case with the polymer in figure 4.7, surprisingly few pairs are required to
accurately generate a scattering profile up to high values of q. When insufficient sampling
occurs, the effects will be most evident in the high-q regions, as this corresponds to the
smallest length-scale features in a morphology and are thus less precisely sampled than
those features which exist on longer length scales. It is recommended that when employing
the DFM algorithm on an unfamiliar system, a careful examination of how sampling rate
effects resultant scattering profiles be performed.
One of the most useful aspects of the DFM algorithm is that once a system’s minimum
sampling rate is established, dramatic speed increases are possible for many classes of morphologies. If one is interested in larger scale morphological features only (i.e. q < 1 a−1 ),
we have found that applying a smearing function to I a (q) can alleviate high frequency oscillations introduced by sampling rates too low to accurately capture the high-q information.
As a demonstration, we generate large (4003 ) lattice box models of 3-d interpenetrating
two-phase morphologies where the goal is to characterize the growth of average cluster size
81

in these models under Ising model simulated annealing[46, 2].
We initially generate structures by randomly assigning sites as one of two materials
assignments, s = 1 or −1. A simple nearest-neighbor site energy Hamiltonian is employed,
Hi = −

j si sj , where a site i with material assignment si will have an energy based on the

sum of it’s nearest neighbor site assignments, sj . Sites are selected at random and evaluated
with Metropolis Monte Carlo, driving the morphology towards increasing phase segregation.
We quantify the degree of phase segregation with the characteristic feature size, a = 3V /S,
where V is total morphology volume and S is total surface area between the two phases.
Further details of the Ising model employed are given in section 5.1 of this thesis. We fix the
planar density along one direction of the morphology so that the system cannot completely
phase segregate, thus keeping a somewhat homogeneous mix throughout the model. DFM
calculations were performed as the system was annealed from a = 4 to a = 16, sampling
with a rate of 106 pairs (approximately 10 orders of magnitude less pairs than required by
the full Debye calculation), using the Gaussian “fuzzying” routine with σ = 120. Once
calculated, the smearing function suggested by Kline for use in neutron scattering data[118]
was employed. The results are shown in figure 4.9. The smearing function used is defined
by

∞

Is (q0 ) =

R(q, q0 ) =

0

R(q, q0 )I(q)dq

1
2
2πσq

1/2

exp

− (q − q0 )2
.
2
2σq

(4.39)

(4.40)

Here, Is , is the smeared intensity, σq is the standard deviation of the resolution in q. In
practice, this smearing function encompass all effects that cause resolution losses in a system.
Values of σq between five and ten percent were often found to best match experiment.
82

Figure 4.9 A demonstration of how simulated scattering can be used to study the effects of
annealing a complex, interpenetrating structure such as the Ising model shown here. Also
demonstrated is the use of applying a smearing function to somewhat under-sampled data.
The raw data is presented as dashed lines and the σq = 10% smeared data is presented as
solid lines. Inset figures show cross-sections of these two-phase morphologies as they are
annealed from a = 4 to a = 16.

83

As demonstrated in figure 4.9, the scattering profile before smearing (dashed lines) and
after a smearing function (solid lines) show the same large trends in plateau position and
overall profile, but the unsmeared data contains significantly more noise and would thus be
difficult to fit to experimental data. If DFM were being employed to observe the drift of the
characteristic plateau, as it shifts to higher intensity and lower-q as feature size is increased,
it is clear that the sampling rate used here would be sufficient and the presented high-q
inaccuracies would be acceptable. The smearing effect itself is less necessary as a increases,
as there are few small scale features left in the morphology.

4.4.4

The effects of Interaction

A novel aspect of the DFM approach is that structure factors and form factors need not be
treated independently. In theory, one considers the form factor, P (q), to be the contribution
to the scattering due to the shape of the individual scatterer and the structure factor, S(q),
to be the contribution due to the arrangement these individual scatterers. The overall
scattering is a product of these two terms

I(q) = S(q) × P (q).

(4.41)

In it’s simplest interpretation, the DFM algorithm treats all scatterers as points, and
thus it is the structure factor that is being calculated. However, if the model resolution
is sufficiently high we may construct “scatterers” from a series of points, and then build a
larger superstructure out of these scatterers to generate a unique structure factor. Here,
we demonstrate this by looking more closely at the case of monodisperse spheres randomly
placed inside a large box with different minimum center-to-center distances.

84

Similar to the morphology modeled in figure 4.8, we randomly place spheres of radius
r = 10 ˚ inside a cubic sample box of side length L = 400 ˚, varying the minimum allowed
A
A
sphere center-to-center distance as D = 20, 40, and 100 ˚. Given the resolution we employed
A
(1 ˚) and that we constructed the morphology on a cubic lattice, each sphere was composed
A
of 4169 individual scatterers. The number of spheres used in each case was varied such that
each configuration was filled to high density. DFM calculations were performed on these
systems, sampling 109 pairs and the resultant scattering profiles were smeared by 5%. For
direct comparison, the scattering profiles were normalized such that I(0) = 1. The results of
these calculations are shown in figure 4.10. For comparison, the form factor for a r = 10 ˚
A
sphere is also shown.
In all cases, the average center-to-center distance, D, is clearly observable in figure 4.10
as a small peak at q = (2π/D), with subsequent damped oscillations. Note that while both
the spherical form factor (clearly visible at q = 2π/r) and structure factor due to the average
center-to-center distance, ≈ D, are clearly distinguishable for the D = 100 and D = 40 ˚
A
cases, the form and structure factor contributions become somewhat indistinguishable when
the average center-to-center distance becomes comparable to the diameter of the spheres,
D ∝ 2r = 20 ˚. While no overlap of spheres is occurring, the effects of the structure factor
A
peak in the D = 20 ˚ case are clearly interfering with the normally well defined r = 10 ˚
A
A
spherical form factor, highlighting the inherent difficulties present in interpreting high density non-crystallized systems, even when they are composed of very regular nanostructures.
This simultaneous modeling of both a form factor and structure factor is possible because
of the fine-grained mesh we are employing here, and such a high-resolution model is computationally simple to probe with the DFM algorithm.
As a final example we model a fixed number of polydisperse spheres inside a large box.
85

˚
Figure 4.10 Spheres of constant radii, r = 10 A, randomly placed inside an L = 400 ˚ box
A
˚. For comparison, the
with a minimum center-to-center distance of D = 20, 40, and 100 A
form factor for a single sphere of radius r = 10 ˚ is also shown. The DFM calculation is
A
performed by sampling 109 pairs, and smearing resultants by 5%.

86

We enforce a minimum center-to-center distance between spheres of D = 100 ˚ using a
A
cubic sample box of side length L = 800 ˚, and fix the number of randomly placed spheres
A
in each case to be N = 370. For each sphere, a radius is selected based on a polydisperse
Schulz distribution with an average of r = 20 ˚ and width σ = 5, 10, and 15. The Schulz
¯
A
distribution being defined[122]

z + 1 z+1 z
P (r) =
r exp −
r
¯

z+1
r
r
¯

1
,
Γ (z + 1)

(4.42)

where r is the average radius, z = r/σ, and Γ(z + 1) is the gamma function defined at
¯
¯
z + 1. The polydispersity index (P DI) of each resultant morphology was measured.
The results of the DFM calculated scattering profiles, along with an equivalent system
with monodisperse spheres (P DI = 1.0) are shown in figure 4.11. Of note is the emergence
of clear polydispersity effects as soon as the sphere radii is permitted to vary, seen as a
natural blurring of the sphere peaks visible at q > .2 in the r = 20 ˚ form factor case. As
A
the polydispersity is increased, the bend over point continues to decrease, as is consistent
with polydisperse sphere form factors. Similar to the results of monodisperse spheres in
figure 4.10, the structure factor contribution of the average center-to-center distance, D =
−1

100 ˚, is plainly visible at q = (2π/D) ≈ .06˚
A
A

for all cases except the most polydisperse,

P DI = 8.15, at which point the polydispersity induced bending occurs at small enough q
that it obscures the structure factor peak.

87

˚
Figure 4.11 Scattering from a series of 370 spheres randomly placed in a L = 800 A box
˚. The sphere radii are selected from
at a minimum center-to-center distance of D = 100 A
a Schulz distribution shown in the inset to the figure. For each distribution, the average
radius is r = 20 ˚, but the deviation used in each case is σ = 5, 10, and 15, which generated
A
resultant polydispersity of P DI = 1.36, 3.27, and 8.15 respectively. For comparison, also
shown is the scattering from spheres placed under similar conditions with constant radius,
r = 20 ˚, corresponding to P DI = 1.0. DFM calculation performed by sampling 109 pairs,
A
and smearing results 5%. The inset shows the Schulz polydisperse sphere distributions used
to select radii.

88

Chapter 5
Nanoscale Morphologies
Incorporating Reflectometry Data
and Assessing Effects on Percolative
Processes
Both characterization of the BHJ structure and a greater understanding of how morphological
effects influence overall OPV device performance are active areas of research[123, 116, 124,
125, 126, 127, 128]. Arguably, the most successful models that relate morphological features
to device performance have relied on Ising models[86, 84, 2, 129], published examples of which
are shown in figure 5.1. While the Ising model is in no way the only acceptable method for
addressing the effects of morphology on device performance, it’s simplicity and versatility
make it ideal for studying the effects of morphology on transport mechanisms present in an
OPV device. In this chapter, we will begin with a detailed description of the simple lattice
Ising model as it has been used for generating BHJ morphologies.
In order to generate more realistic model morphologies, we introduce modifications to
the simple Ising model which incorporate experimental data by refining the structure using
Monte Carlo methods. Specifically, we will demonstrate models which incorporate neutron

89

reflectometry data taken on P3HT/PCBM BHJ films. We will present derivations of the
basic principles of reflectometry and demonstrate how a density profile may be derived from
the reflectometry data[130]. The derivations presented are taken from texts[110, 121] unless
otherwise noted.
We will present 3-d interpenetrating, percolative morphologies consistent with the experimental reflectivity data, and compare DMC performance calculations performed on these
realistic morphologies to results from previous simple “Ising only” models, which contain
no reflectivity data. We will then demonstrate modifications to the model to examine the
effects of P3HT crystallization and degree of PCBM sequestration.
Finally, we will show that reflectivity data on it’s own is insufficient to select a morphology
for use in DMC calculations, as a large number of qualitatively varying morphologies may
be generated to fit the reflectivity data to high precision.

5.1

Simple Three Dimensional Ising Model

Originally formulated to model the statistical dynamics of ferromagnetism, the Ising model
presents a straightforward and efficient method to generate two-phase systems which are
formulated with a simply defined configuration based Hamiltonian. The principles of a spinflip Ising model are that each site of a lattice is assigned a designation (say, 1 or -1) which
define the site as being one of two possible states/spins/etc. We use the word “spin” interchangeably with “material” in this discussion, as they are only meant to designate site
assignment. Model morphologies are generated through the use of Metropolis Monte Carlo
(MC) procedures, wherein potential site flips are evaluated based on their overall effect on
a defined Hamiltonian, H(a), where the evaluated energy will be based upon the specifics

90

(a)

(b)
Figure 5.1 Figure 5.1a demonstrates the work of Steve Forrest’s group using a simple Ising
model to roughen the interface between a set of bilayers[1]. Figure 5.1b demonstrates the
work of Alison Walker’s group using an Ising model to study the effect of BHJ feature size
on device performance[2].

91

of the configuration. The system is modified via careful acceptance or rejection of spin
flips based on randomly chosen sites. Candidate spin flips generate altered configurations,
b, whose ideality is evaluated in comparison to the initial configuration, a, via the defined
Hamiltonian as a change in system energy ∆Eab = H(b) − H(a). Candidate flips which
generate a more energetically favorable configuration such that ∆Eab < 0 are always accepted. Candidate spin-flips which generate less favorable configurations (∆Eab > 0) are
only accepted with a certain probability based on the magnitude of ∆Eab and the simulation temperature. Through the use of the well known Metropolis algorithm[131]. The overall
acceptance probability can be summarized as



Paccept (∆E) =

1


 exp − ∆E
kb T

if ∆E < 0
(5.1)
otherwise,

where Paccept (∆E) is the probability of accepting a candidate spin-flip, ∆E is the energy
change of a system due to this candidate spin-flip, and kb T is the simulation temperature. In
the high-temperature limit, almost all flips will be accepted where as in the low-temperature
limit, no energetically unfavorable flips from the ground state will occur.
The nearest neighbor spin 1/2 Ising Hamiltonian function, Ho , is defined as
n.n. sites
Ho (Si ) = −
JSi Sj ,

(5.2)

i=j

such that the energy contribution of a site, i, which has spin assignment Si , is related
to the interaction strength, J, with nearest neighbor sites j, that have spin assignment, Sj .
Sites may have spin designation of S = +1 or S = −1. With such a Hamiltonian, values of
J > 0 will coarsen the morphology, minimizing surface area between the two spin-phases. If

92

J < 0, the system will tend towards a mixed configuration, maximizing surface area between
the two spin-phases. A common metric used to describe the current state of a two-phase
morphology is the characteristic feature size, a, defined as

a = 3V /S,

(5.3)

where V is the total volume and S is the total surface area between the two spin-phases.
Smaller values of a correspond to a greater level of mixing between the two materials, whereas
larger values of a designate a coarser grained-structure.
A time-step in an Ising model MC process is defined to be attempting to make a number
of random site spin flips equal to the number of sites in the lattice, N = Lx × Ly × Lz for a
3-d cubic lattice. Candidate sites are selected at random for each potential flip, rather than
a constant sweep through the morphology as this would bias the structure. This randomselection method means that in a single MC step, some sites may not be given an opportunity
to flip while others could potentially flip many times. One can correct for this by initially
generating a list of all sites in a lattice, randomly ordering this list every time step, and going
through this list to select candidate sites; but the results of the two methods are essentially
the same.
The total number of time-steps required to anneal an Ising morphology to a desired
characteristic feature size, atarget , will be dependent on the specific Hamiltonian used and
the simulation temperature. For values of J > 0 the morphology feature size evolution
will grow as T 1/2 where T is the number of Ising steps, as is demonstrated in figure 5.9.
Higher values of T will permit the morphology to more quickly explore configurational phasespace, but one must be careful as too high a temperature will prevent the morphology from

93

converging on the ideal energetic minima configuration.
Two versions of the Ising model MC algorithm are the Kawasaki (spin-swap) dynamics method and the Glauber (non-conserved) Monte Carlo (GMC) method. In the GMC
method[132], a single site’s contribution to the Hamiltonian is considered for each spin flip,
and if that flip is found acceptable (passes the Metropolis test) then the site’s spin-assignment
is flipped. In the Kawasaki method[133], two sites of different spin values are considered simultaneously, having their individual candidate spin-swap energy contributions summed to
generate ∆E. If the candidate swap is accepted, then the two site spins are exchanged,
effectively exchanging their spin assignments. In the Kawasaki method total spin is thus
conserved, which is not guaranteed in the GMC method. By nature of single site evaluation
compared to pairs of sites, the GMC method will anneal a structure to coarser grained structures (greater a) more rapidly, but the lowest energy configuration will be one in which all
sites are assigned a single spin value. However, by adding a term to the Hamiltonian defined
by equation 5.2, which is a function of total population density, we may overcome this and
find structures with a target population. An example of such a modified Hamiltonian is
n.n. sites
p
Hconfig (Si ) = −
JSi Sj + Jpop νconf ig − νtarget ,

(5.4)

i=j

where the prefactor Jpop and the parameter p constrains the current population of a spin
species, νconf ig , to the target population, νtarget .
As Ising models have been used extensively to describe nanomorphology present in OPV
devices, our goal is to improve the techniques in order to incorporate experimentally provided information about the morphology. Reflectometry is one experimental method which
lends itself easily for incorporation into an Ising model, and we introduce the principles of

94

reflectometry measurements here.

5.2

Principles of Reflectometry

z
i

i

i

i

x

t
t

y

Figure 5.2 Schematic of reflectometry measurement, with incident wavelength, λi reflecting
at an angle of θi off of the surface, with a transmitted wave having associated wavelength
λt at angle θt . Note that this figure demonstrates the reflection and transmission angles
associated with only the first layer of a multilayer system. Further reflections/transmissions
are implied but not shown.
Reflectometry is a measurement which probes the scattering properties of a very flat
sample, typically described by a 1-d vertical function normal to the surface of the sample.
In a reflectometry experiment, a beam of particles (neutrons or x-rays) is incident nearly
parallel to the surface of a sample (see figure 5.2). Based on the properties of the sample,
a fraction of the incident beam intensity will be specularly scattered off the sample, being
reflected symmetrically forward at the incident angle, θi , which is then measured with a
detector placed far from the sample. That fraction of the beam which is not reflected is
transmitted into the sample, at angle θt , and based on the material will have an associated
transmitted wavelength, λt . The measurement thus characterizes reflectance from a sample
with values ranging 0 ≤ R ≤ 1, where 0 represents total absorption in the sample and 1
95

being total reflection. The incident beam is characterized by it’s momentum transfer, q,
which is defined

q=

4π sin θi
.
λi

(5.5)

The measurement will thus inform on a sample’s reflectivity as a function of momentum
transfer, R(q). As the beam is incident over a large surface area of the sample, all structural
information due to lateral (x-y dimension) features will be averaged together and only those
details in the vertical (z-dimension) are discernible.
The reflectivity measured off a sample will be related to the scattering strength of the
materials in the sample, described most directly by their scattering length density (SLD).
Standard reflectivity data analysis involves fitting an experimentally measured reflectivity
profile to a modeled SLD profile, β(r), which due to the lateral averaging of the sample may
simply be considered as β(z).

5.2.1

Kinematic Derivation

To begin with, let us examine the Fourier transform relationship between the elastic differential cross-section and the SLD function. Recall the relationship derived in equation 4.25,
we may write this more straightforwardly as
2

dσ
∝
dΩ

dr β(r) eiQ·r ,

(5.6)

V

noting again that the proportionality is used instead of an equals sign as the RHS has
not been normalized by a constant such as the number of scatterers, mass, or volume, as
this normalization varies across disciplines but this core concept relation remains the same.
96

If we define a sample to have dimensions (2Lx , 2Ly , 2Lz ) (with lengths selected to simplify
the calculation), and assume β(r) is equal to β(z) inside the sample, and 0 outside, we may
rewrite equation 5.6 as

dσ
sin2 (Lx qx ) sin2
∝ 16
2
dΩ
qx

∞

Ly qy
2
qy

−∞

dz β(z)eizqz

2

,

(5.7)

where the two sinc-squared prefactors come from the scattering of a rectangular slab of
width 2Lx × 2Ly . Thus, the scattering will have a maximum at qx = 0 and qy = 0 of the
value 16L2 L2 and be focused in the small region qx/y < L π .
x y
x/y
The relationship between elastic differential cross-section and the specular reflectivity
measured in a reflectivity measurement is best understood by comparing the definitions

dσ
Number of particles deflected an angle (2θ, φ) per unit solid angle
=
dΩ
Number of incident particles per unit area of beam
and

R=

Rate of specular reflective scattering
.
Rate of incidence

Clearly, both values are based on the ratio of scattered particles to the illuminated portion
of the sample, which in the case of the reflectometry is 4Lx Ly sin θ, whereas the differential
cross section is integrated over a solid angle, ∆Ω. Thus, in the range of specular condition,
we may write the reflectivity in terms of the differential cross-section

R(Q) =

1
4Lx Ly sin θ

dΩ
∆Ω

Using the substitution
97

dσ
.
dΩ

(5.8)

dΩ
∆Ω

dσ
∆Ω
dσ
≈
×
dΩ
4
dΩ qx =0,qy =0,qz =−q

(5.9)

and approximating ∆Ω as

16π 2 sin θ
∆Ω ≈
.
Lx Ly q 2

(5.10)

Assuming that the area under illumination is large compared to the wavelength of the
incident particle, this yields the relation

R(q) =

16π 2
q2

∞
−∞

dz β(z)e−iqz

2

.

(5.11)

Notice also that unlike the differential cross section definitions, which were left as proportionalities, the reflectivity is a direct equality as the former was a consequence of a variety
of conventions used in normalization, but this is not required for the SLD as one may define

β=n b ,

(5.12)

where the SLD, β, is defined as the product of the atomic density (atoms per unit
volume), n, and the coherent scattering length, b . Thus, utilizing integration by parts we
can redefine the reflectivity data as a function of the derivative of the SLD

16π 2
R(q) ≈ 4
q

∞

dβ −iqz 2
.
dz
e
dz
−∞

(5.13)

Here we have written the result as only an approximation, not direct equivalence. As we
will show, this is because any result derived using the kinematical assumption is not valid

98

for reflectometry, particularly at low values of q.

5.2.2

Dynamical Derivation

While an excellent first order approximation, the Fourier relationship derived in equation
5.13 was formulated under the Born or kinematical approximation, which assumes that
the scattering process is weak and that the scattered wave has a negligible effect on the
incident beam. In such a kinematical picture, multiple scattering is ignored, hence the
simple superposition of the scattering due to individual scatterers. In neutron reflectometry,
however, measurements are taken in a regime where all incident neutrons are specularly
reflected from the surface of a sample, which is not weak scattering and thus, not in the
Born regime.
The perturbation of the incident wave by the scattering process is considered negligible
in the Born approximation. However, a “thick” sample is used in neutron reflectometry
due to the path the neutrons take being L = T / sin θ where T is the thickness and θ is the
incident angle. This leads to the failure of the approximation of equation 5.13, particularly
as θ becomes close to 0 or the thickness, T gets large. Deviations at higher λ values come
from the reciprocal relationship between momentum of an incident particle and wavelength,
meaning that the more energetic a particle is, the less it will be perturbed.
In order to generate accurate models of reflectivity, one must use a dynamical theory of
scattering, rather than a kinematical (Born based) one used previously. We present that
derivation here.
Using the standard definition of a refractive index

99

i

t

z

n

air
r
(a)

-

t

-

-

N

2

+

+

+

N

1

i

1

2

z
2

N

S

zN

z2

1

z1

r

air

0

(b)
Figure 5.3 Two qualitative 1-d examples of an incident wave, ψi , which originates in air
(βair = 0) interacting with a sample. Figure 5.3a demonstrates a incident wave ψi incident
on a sample with index of refraction n, with subsequent transmitted wave ψt and reflected
wave ψr . Figure 5.3b demonstrates the incident wave on a multilayered sample, with N
different layers atop a substrate layer S. Each layer has a propagating component, ψ − to
the left and a reflected component ψ + to the right.

100

n(λi ) =

cos(θi )
λ
= i,
cos(θt )
λt

(5.14)

where the index of refraction is considered a function of incident wavelength, λi , it is
defined as the ratio of the cosine of the incident angle, θi , over the cosine of the transmitted
wave angle, θt . This is equivalent to the ratio of the incident wavelength over the transmitted
wavelength, noting that the frequency of the wave will be unchanged through the boundary.
To derive the dynamical scattering expression we begin with a simple one dimensional
single boundary example, as defined in figure 5.3a, with the interface at z=0. The incident
wave is defined

ψi (z) = ψo e−izki sin θi ,

(5.15)

where ψi is the incident wave as a function of coordinate z, ψo is the incident amplitude,
ki is the incident wavenumber defined as ki = 2π . At the boundary, a fraction of the wave
λi
will be transmitted and the rest reflected. We may define these transmitted and reflected
waves similarly as

ψt = te−izkt sin θi

(5.16)

ψr = reizki sin θi ,

(5.17)

where the coefficients, t and r, define the amplitude of transmitted and reflected wave
respectively. The sign convention denotes propagation direction, as is seen in figure 5.3a.
Using equation 5.14, the transmitted wave vector can be written kt = 2πn . Similarly, the
λi
101

relation of the incident angle to reflected angle can be written as 1 − sin2 θt =

1−sin2 θi
. The
n2

boundary conditions which must be satisfied are thus

ψi + ψr = ψt

(5.18)

d
d
(ψi + ψr ) = ψt ,
dz
dz

(5.19)

the solution of which can be found to be

t=

2ψo
and r =
1+α

1−α
1+α

ψo

(5.20)

where
kt sin θt
n
(1 − sin2 θi )
α=
=
1−
ki sin θi
sin θ
n2

1/2

.

(5.21)

Clearly, we can see that total reflection (r = 1) will occur when α is equal to 0. Additionally, if the term in brackets of equation 5.21 is less than 0, then α will be imaginary, which
cannot occur if any wave is to be transmitted. We may thus define the fraction of reflected
wave as the ratio of the amplitude reflected to the incident amplitude, squared, as

R=

r 2 1 − 2 (α) + |α|2
.
=
ψo
1 + 2 (α) + |α|2

(5.22)

Total external reflection will occur when R = 1. Thus, in order to have total external
reflection the necessary condition exists

sin2 θi ≤ 1 − n2 .

102

(5.23)

Complete reflection thus occurs for materials with n < 1 (for an interface with the
incident wave passing through air) below a critical angle of

θc = sin−1

1 − n2 .

(5.24)

Recalling the definition of momentum transfer, q = 4π sin θi , and the condition given by
λ
equation 5.23, we may rewrite equation 5.21 as

α=

n
(1 − sin2 θi )
1−
sin θi
n2

1/2

1/2
n
n2 − (1 − sin2 θi )
n sin θi
1/2
1
−(1 − n2 ) + sin2 θi
=
sin θi

=

1/2
sin θi
n2 − 1
=
1−
sin θi
sin2 θi
1 n2 − 1
=1−
2 sin2 θ
8π 2 (1 − n2 )
≈1−
.
λ2 q 2
i

(5.25)

Note that the last step of this derivation takes advantage of the binomial expansion,
√

1 − A = 1 − A/2 − . . ., and this approximation will be more accurate as q increases.

Noting that α ≈ 1, we may rewrite equation 5.22 in the form

R=

2
1 − α 2 16π 2 π(1 − n2 )
,
≈ 4
1+α
q
λ2

(5.26)

which when compared to equation 5.13, implies that the term in brackets must be equated
with the SLD, such that
103

n2 ≈ 1 −

βλ2
i.
π

(5.27)

Note that this is the same result Fermi derived for the transmission of a signal through
a planer slab of uniform material in 1950. Referring back to the critical angle as defined in
equation 5.24, we see that we may define a critical value for q where total external reflection
will abruptly cease

qc = 4 βπ.

(5.28)

Having derived the dynamics of reflection and transmission for a single boundary, we
may generalize to treat multiple boundaries. Having a series of boundaries, defined by their
SLD profiles


 βs




β(z) =
 βj




 0

for z < zN
for zj < z < zj−1

(5.29)

for z > 0,

where j = 1, 2, 3, . . . , N are the indices defining the N layers of the structure, with an
air interface at z = 0, as seen in figure 5.3b. The wave function in a single slab, j, may be
written

izkj
−izkj
+ Bj e
,

ψj (z) = Aj e

(5.30)

where Aj and Bj are the amplitudes for waves in the reflected (positive in z) and transmitted (negative in z) direction, respectively. The wavenumber for each component may be
written
104

kj =

2π
λj

sin θj =

2πn)j
λi

1−

1 − sin2 θi
.
n2
j

(5.31)

In a similar method to the single layer, we may rewrite the wavevector for the j-th layer

kj =

1
2

q 2 − 16πβj .

(5.32)

The boundary conditions for this set of wavefunctions is

ψj+1 (zj ) = ψj (zj ) and

d
d
ψj+1 (zj ) =
ψ (z ).
dz
dz j j

(5.33)

The solution, similar to equation 5.20, has the form

Aj =

ikj ψj (z) + ψj (z)
izkj

2ikj e

and Bj =

ikj ψj (z) − ψj (z)
−izkj

.

(5.34)

2ikj e

The reflectivity is measured as the ratio of reflected amplitude to transmitted amplitude
at the surface interface (j = 0), will be given by |A0 /B0 |2

ik0 ψ(0) + ψ (0)
R(q) =
ik0 ψ(0) − ψ (0)

2

.

(5.35)

As transmission occurs when AN +1 = 0 (implying no reverse wave inside the substrate),
we may write that at the substrate the wave is defined

ψ (zN ) = −iks ψ(zN ),
which has the solution

105

(5.36)

ψ(zN ) = 1 and ψ (zN ) = −iks .

(5.37)

Starting with this solution, one may recursively solve the wave function inside each
previous layer, ψ(zj − 1), as a function of the following layer with the relationship

ψ(zj − 1) = cos(kj Tj )ψ(zj ) +

sin(kj Tj )
ψ (zj )
kj

(5.38)

and
ψ (zj − 1) = −kj sin(kj Tj )ψ(zj ) + cos(kj Tj )ψ (zj ),

(5.39)

up to j = N where the surface of the multiple layers is reached at z = 0. This iterative
process for calculating reflectivity was proposed by Abel´s and Parratt, and this is the
e
method used by many modern reflectivity calculations, such as those demonstrated in figures
5.4 to 5.7.
Examples of the relationship between SLD density profiles and associated reflectometry
data are shown in figures 5.4 to 5.7, which have all been generated with the commonly used
motofit macros[134] for IGOR. In all examples of figures 5.4 to 5.7 and discussions here, SLD
−2
values have units of 10−6 ˚ . In the first example, figure 5.4, a single step of thickness 25
A

˚ (βstep = 5) between an air interface (βair = 0) and a substrate (βsub = 5) is shown with
A
different step “roughness” of 0, 5, and 10 ˚, where roughness is a somewhat counterintuitive
A
term used in the motofit macros to denote a blending transition between two adjacent layers
of an SLD density profile. The first important features of note in figure 5.4 is the sudden
√
drop from total external reflection (R = 1) at the critical angle, qc = 4 πβ, at which point
the reflectometry takes on the R(q) ∝ q −1/4 dependence as seen in equation 5.13. At angles

106

Figure 5.4 A demonstration of the simulated reflectometry data from a simple step function,
with various amounts of “roughing” of the edges to smooth the transition between layers.

Figure 5.5 An example of various step functions with combinations of sharp and smooth
transitions between steps, demonstrating how identical reflectometry data can be generated
from structures with different SLD density profiles.

107

Figure 5.6 An example of simulated reflectometry data from layers of alternating SLD
values, in order to demonstrate how SLD contrast effects reflectometry data.

Figure 5.7 A demonstration of the effect of varying the background layer of an otherwise
fixed SLD density profile.

108

lower than the qc , neutrons will be entirely reflected from the interface, thus no information
on the internal structure is probed. This initial drop at qc often provides a clear picture
on the SLD value of the first layer the neutrons pass though, so long as the SLD value
of this material is positive. The second feature to note from figure 5.4 is that increasing
roughness causes a further Gaussian like decay from the R(q) ∝ q −1/4 relationship, which
can be described as

R(q) ≈ Ro (q) × exp(−σ 2 q 2 ),

(5.40)

where Ro (q) is the reflectometry due to a sharp interface (no roughing) and σ is the
thickness of the associated Gaussian which has been convoluted with the sharp interface
structure to generate a smoothed transition between two adjacent layers.
The second example, seen in figure 5.5, exhibits different transitions between an air interface, 25 ˚ thick layer of βstep = 2.5, and substrate (βsub = 5). There are essentially
A
two transitions which must be made, from air to the step layer (∆β = 2.5) and from the
step layer to the substrate (∆β = 2.5), with the effects of different combinations of these
transitions (sharp interfaces or smooth) being shown. Note that while the two sharp interfaces (black line) decay the least at higher q, while the two smooth interfaces (green line)
decay the most, the mixed transition states of sharp-smooth or smooth-sharp (blue and
orange, respectively) are indistinguishable from one another. This highlights an important
shortcoming of the reflectometry method, that different SLD-density profiles may generate
identical reflectometry graphs, as the relationship between them is based on the integral of
all SLD derivatives (dβ/dz) throughout the sample, as seen in equation 5.13, and thus R(q)
contains no inherent information about the ordering of the layers, nor the SLD values inside

109

the system, only the transitions between layers.
The third example in figure 5.6 shows the reflectometry for a system of 6 total layers between air (βair = 0) and a substrate (βsub = 5) which have three different sets of alternating
SLD values for the layers, of 3 and 3 (black line), 1 and 3 (blue line), and 4 and 3 (green
line). This example demonstrates the importance of contrast difference between components
in a complex system, as the greater the difference in SLD values of the layers is, the more
significant the deviation of the R(q) profile from the uniform case. Note that the sign of the
difference doesn’t matter, only the magnitude, a result of the loss of phase information in
the scattering measurement. Contrast is often a principle concern in scattering experiments,
leading to the frequent use of deuterated samples in neutron based studies as deuterium has
a significantly different SLD value than hydrogen (≈ .67 to ≈ −.37 respectively).
The final example in figure 5.7 is a specific demonstration of altered contrast in a system,
wherein an arbitrary complex layer between the air and substrate is held constant, but the
substrate SLD is altered from β1 = 5 to β2 = 2. This subtle change to the system can be seen
to cause a dramatic shift in R(q) at all values of q. This highlights a technique of backinglayer contrast change, in which an unknown morphology is probed using a variable backing
layer, such that two measurements for the same system are taken with a well known difference
between them. The resultant two reflectometry profiles, R1 (q) and R2 (q), may then be solved
with the constraint that the density profile must be consistent with both measurements, as
only the backing SLD value will have altered. This method produces density profiles of
significant confidence, as the morphologies are fit not just to a single measurement, which
we know to be underconstrained (as shown in fig. ??) but two measurements simultaneously.
This is the method utilized by Jon Kiel to generate reflectometry fits from the annealed and
as cast P3HT/PCBM solar cells[8], the results of which are shown in figure 5.10.
110

5.3

Generation of BHJ Morphologies Consistent with
Reflectometry Data

First and foremost, the models we generate should be consistent with the general picture
of these OPV BHJ morphologies. We first motivate our structures by presenting recently
published images and models of various polymer/fullerene BHJ morphologies, shown in figure
5.8. Though subtle differences between sample preparation or model picture will vary from
group to group, a clear consensus is that the nanomorphology present in these BHJ structures
is one of 3-d interpenetrating, percolative structure.
The bulk heterojunction morphology has been modeled extensively in previous work
through the use of a simple spin-flip Kawasaki lattice Ising model[2, 47], with sites being
assigned as either pure PCBM (acceptor) or pure P3HT (donor). For our initial models, we
employ a methodology similar to that used previously, with the addition of a constrained
Ising model to incorporate an experimentally provided target density profile, ρtarget . At this
point we impose no other constraints on the morphologies, such as P3HT crystallization or
PCBM sequestration. The Hamiltonian we use is

n.n. sites
p
Hconfig (Si ) = −
JSi Sj + Jpop νconf ig − νtarget 1
i=j
L

+Jdens

0

dz |ρconf ig (z) − ρtarget (z)|p2 ,

(5.41)

where the Hamiltonian evaluating a specific site i with material assignment, Si , is defined
by three terms. The first term defines the phase-segregation energy, J, looking at nearest
neighbor material assignments, Sj . The third term, defined by the energy Jdens , evaluates
111

(a)

(d)

(b)

(e)

(c)

(f)

Figure 5.8 Various images and models of polymer/fullerene BHJ morphologies. Figure 5.8a shows bright field TEM images of an polymer/fullerene BHJ OPV film composed of PTB7:PC61BM/DCB + DIO, with artificial black lines added to highlight
differences between the polymer-rich and nanoparticle-rich domains[3]. Figure 5.8b is
a cartoon illustrating the nanomorphology proposed in BHJ OPV devices under the
“Rivers and Streams“ model[4]. Figure 5.8c is an image of the x-y plane in a polymer(PF10TBT)/nanoparticle(PCBM) film generated from electron tomography data, with
the polymer regions highlighted in light and PCBM regions in dark[5]. Figure 5.8d shows
a sample cross-section of reconstructed data taken on a P3HT/ZnO OPV device obtained
through electron tomography with ZnO being yellow and P3HT being transparent[6]. Figures 5.8e and 5.8f show sulfur maps and carbon maps, respectively, taken through energy
filtered electron tomography (EFTEM) for 1:1 weighted P3HT/PCBM films which were annealed for 30 minutes at zero defocus and with dimensions approximately 400 nm across.
Light regions in the sulfur map correspond to higher concentrations of P3HT-domains[7].

112

Figure 5.9 Demonstration of Ising annealement over time on an example morphology annealed with parameters J = 2.0, kB T = 2.0, Jdens = 15.0, and Jpop = 4.1, fitting the
annealed reflectivity profile shown in figure 5.10. As the morphology was being annealed,
exciton transport calculations were performed for each integer value of feature size reached
(a = 1, 2, 3, . . .) using a simplified form of the DMC model as is described in chapter 6.
The results of these calculations are presented in the inset as exciton dissociation efficiency
versus characteristic feature size (a). The calculation was performed using decay of lengths,
Lx = [.01, .1, 1, 10, 100], which are shown. Four specific morphologies are highlighted, with
their corresponding results also highlighted in the inset.
how closely a candidate configuration’s in-plane density, ρconf ig (z) adheres to the target
density profile, ρtarget (z). The final term, defined by the energy Jpop , evaluates how closely
the candidate configuration’s target population density, νconf ig , adheres to the target as
provided by integration of the target density profile, νtarget . The terms p1 and p2 are fitting
parameters which for the models presented here we set to unity, but which can be varied to
tune the strength of the constraints.
With the addition of the second term (defined by Jpop ) we may change from the previously utilized pair-swapping Kawasaki Ising method to the non-conserved single-site flipping
113

Glauber Metropolis Monte Carlo method. The introduction of this third term adds a constraint such that the population of the entire morphology will not deviate far from the target
population value, ρtarget .
So long as Jdens is sufficiently large, all generated morphologies will fit the reflectivity
profile to high precision. For our models, we utilize PCBM density profiles derived from
neutron reflectometry measurements taken by Jon Kiel on spin-cast P3HT/PCBM BHJ
OPV mimics before and after annealing, the fits for which are shown in figure 5.10. At this
stage, the model contains no information which informs on the target characteristic feature
size (a), so we generate many morphologies with a wide range in a. 2-D cross-sections of
some morphologies generated using this method are shown for the as-cast morphology in
figure 5.11 and the annealed morphology in figure 5.12. Clearly, this method generates
morphologies which are qualitatively similar to previous models, but also consistent with
neutron reflectometry data.

Figure 5.10 PCBM density profile fits derived from neutron reflectivity data and an associated cartoon theorizing the morphology[8].

114

The cross-sections illustrate that as characteristic feature size, a, is increased, the model
morphologies become increasingly phase segregated. The effects of the imposed density
profile are immediately visible, in that while phase-segregation is clearly perceivable with
increasing values of a, the material is tethered to span the entirety of the morphology because
of the target density profile. The details of these reflectivity informed density fluctuations
are seen most dramatically in the sudden drop in PCBM (acceptor) sites on the top electron
collecting (metal) side of the as-cast morphologies, which is clearly seen comparing the topdown views of the as-cast and annealed morphologies in figures 5.14a and 5.14b respectively.
More subtle effects of the imposed density profiles are perceptible, as higher density PCBM
sites (brown) correlate with higher PCBM densities in the target profile.
For a straightforward comparison between the models generated here, and the experimental picture of BHJ morphology, we take a 400 nm × 400 nm section of a published image of
electron tomography of a BHJ morphology[5] which is similar to those used in the Kiel reflectometry study, artificially color it similarly to how we have presented our models, and present
the results in figure 5.13. Note that the polymer used in the tomography study was poly[2methoxy-5-(30 ,70 -dimethyloctyloxy)-1,4-phenylene vinylene] (MDMO-PPV), and not the
P3HT used in the neutron reflectometry studies here, so this example is shown only for
qualitative comparison of a polymer/fullerene BHJ nanostructure.
While these models are similar to previous Ising model morphologies, generating 3-D
interpenetrating percolative morphologies, the method presented provides a simple and easily expandable framework on which to incorporate additional details into the model. To
demonstrate this, we first explore the theoretical effect of PCBM sequestration, as proposed
in the “Rivers and Streams” model[4].

115

5.3.1

The Three-Phase BHJ Picture

Though much work has been done on the P3HT/PCBM bulk heterojunction solar cell with
a so-called two-phase model, describing the morphology as either a site of acceptor (PCBM)
or donor (P3HT) material, recent work has suggested that the BHJ internal morphology
is better described by a three phase model. Yin and Dadmun[4] suggest a picture of the
morphology they refer to as as the “Rivers and Streams” model, with three distinct components; a phase of pure aggregate PCBM particles, a phase of pure crystallized P3HT, and a
mixed phase of non-crystallized P3HT with dispersed PCBM. A principle assumption in this
model is that charge transport is dominated by transport through the pure material phases
(non-mixed), with hole mobility greatly enhanced in the crystallized polymer and electron
mobility higher in the aggregate PCBM regions. The colloquial name, “Rivers and Streams”
comes from Yin and Dadmun’s theory that charges flow through small “streams” of pure
material percolating through the bulk of the mixed phase, connecting with large “rivers” of
pure material which connect more directly with the contacts. Yin and Dadmun propose that
the fraction of dispersed PCBM is around 20% the total PCBM in a morphology.
Here, we explore the effects of the three-phase model picture using our constrained Ising
morphologies. To begin with, we reduce the complexity of the three-phase model by treating
the crystallized P3HT and mixed phases as one and the same. Yin and Dadmun themselves
suggest that crystallized P3HT phases likely do not occur within films as thin as those used
in OPV devices[4], requiring films an order of magnitude thicker in order to be observed with
small angle x-ray scattering. Additionally, the difference between crystallized and amorphous
P3HT are imperceptible in neutron scattering studies, as the SLD for the two are very similar.
For these reasons, we model two-phases for the “Rivers and Streams” morphologies; a phase

116

of pure aggregate PCBM and a phase of dispersed PCBM in P3HT.
The density profiles, provided by the Kiel neutron reflectometry data, are assumed to
represent the features of aggregate PCBM throughout the film and we assume that PCBM
present in the mixed phase is evenly dispersed. We define a variable, x, which defines
the fraction of total PCBM in a film that has become dispersed in the amorphous P3HT,
and reduce the annealed PCBM density profile evenly at all sites by this dispersed PCBM
fraction, x, as shown in figure 5.15. Using these modified aggregate PCBM density profiles,
we generate morphologies using the modified Ising model as described in section 5.3 and
select those morphologies which have a characteristic feature size of a ≈ 16 (selected for
reasons explained in section 6.1.1) in order to compare to the previous annealed morphology.
We generate morphologies with dispersed PCBM fractions of x = 0, x = .1, x = .15 and
x = .2. Cross-sections of the generated morphologies are shown in figure ??, which upon
initial inspection reveal only subtle differences. However, closer examination reveals that as
the faction of dispersed PCBM (x) is increased, a discontinuity develops between the network
of agglomerate PCBM sites and the electron collecting electrode. To visualize this effect,
we superimpose cross-sections 15 layers thick, labeling those sites which have a connected
aggregate PCBM path to the electrode of interest in green, those sites which contain no such
path and are thus disconnected from the contact are shown in white, and removing all sites
which are the mixed PCBM/P3HT phase, the results of which are shown in figure 5.17.
Qualitatively, the cross-sections in figure 5.17 reveal several initial features of note, foremost being that the majority of the aggregate PCBM in the x = .2 morphology is disconnected from the electron accepting contact. Comparing the x = .15 and x = .2 cases, we
see that a sudden and dramatic transition occurs, implying that a percolative threshold for
this sort of morphology was crossed. Note also that the x = 0 case does have aggregate
117

PCBM sites which are disconnected from the contact (marked as white in figure 5.17) but
that the size of these non-viable PCBM aggregates are much smaller, typically on the order
of a single site, as compared to all other cases where x > 0 in which the non-viable PCBM
has formed larger aggregates.
The qualitative differences between the cases are a direct consequence of having less
PCBM available as the value x is increased, because the models compared here were all
generated to have a feature size of a ≈ 16. In each case, the volume is fixed by the dimensions
of the morphology (V = 4003 ) yet there are a decreasing number of aggregate PCBM sites
overall, meaning on average, each site must be more likely to be an interfacial site and
subsequently less of the sites can be interior bulk of aggregate PCBM (directly adjacent to
only other aggregate PCBM sites). This can be seen somewhat in the “beading” of aggregate
PCBM as x increases. For a similar reason the distance between aggregate PCBM clusters
must increase.

5.3.2

Crystallized P3HT Fiber Networks

Recent studies have shown that transport through crystalline polymer regions are an important factor in overall OPV BHJ device performance[135]. Extensive study of crystalline
P3HT (c-P3HT) has been performed[136, 30, 4, 137, 138, 139], which found crystalline
structures forming with cross-sectional areas of 20 nm × 4 nm and lengths of up to 400 nm.
The molecular structure and packing of individual fibers is well understood. The packing
of fibers is highly dependent on processing conditions, with annealing typically leading to
higher crystallinity. Here, we explore the morphological effects of generating models which
directly incorporate crystallized fiber networks through a templating procedure.
We begin by seeding the morphology with semi-aligned non-overlapping c-P3HT fiber
118

regions. The fibers are added layer by layer with each layer being 4 lattice units thick (in
the z-dimension). The orientation of each layer alternates between spanning the x and y
dimension. For each layer, an initial propagation angle is randomly selected, θ ∈ − π to π
4
4
radians, though each individual fiber may deviate by a random angle of ∆θ ≤ .3 radians. If
a fiber would overlap with a previously placed fiber, it is rejected. All fibers are placed down
segment by segment, with each segment having the dimensions of 20 × 4 lattice units2 . As
segments are being placed, there is a random chance of breaking the crystal fiber, taken here
to be .25%. A break causes a gap to form in a specific fiber, with this gap having a random
width based on a Gaussian distribution with average size 60 and width 20. The c-P3HT fiber
segments resume placement once the break gap is completed. This break behavior leads to
fibers with an average length of 250 nm. Figure 5.18b demonstrates the c-P3HT fiber lattice
structure, from a top-down perspective. c-P3HT fibers are placed down in this way up to
a predetermined maximum density for a given layer set by the P3HT profile derived from
the reflectivity data in figure 5.10 up to a maximum of 30%, generating a density of c-P3HT
such as that seen in figure 5.18a.
Once the initial c-P3HT template is constructed, the remaining P3HT sites required
to match the density profile of the model with the target density are randomly placed on
each layer. A model morphology seeded in this way is shown in figure 5.19. At this point,
Ising annealement is performed on the morphology in order to generate structures with
greater characteristic feature sizes. This annealement is similar to the procedure described
in section 5.3, however those sites which are marked as c-P3HT (shown in blue in figures
5.19) are frozen. The interaction energy, J, between c-P3HT and other sites is taken to
be the same as the interaction between amorphous P3HT sites with some adjustment to
the parameters used in the Hamiltonian (given in equation 5.41), principally using higher
119

temperatures in the simulation, we are able to again generate model morphologies with the
characteristic feature size of a = 16 (this value was selected for reasons discussed in section
6.1.1 of this thesis). For calculations of surface area required to define a, we calculate the
surface area between all PCBM and P3HT sites. This annealed c-P3HT morphology is shown
in figure 5.20, again with PCBM sites shown in brown, amorphous P3HT shown in white,
and crystallized P3HT shown in blue. The c-P3HT templating procedure has introduced a
clearly visible anisotropy into the system, whereby the static framework of the crystallized
polymer network forces additional constraints into how the system may anneal.

5.4

Conclusions

In this chapter, we have introduced a simple method to generate 3-d interpenetrating, percolate nanostructures which are consistent with neutron reflectivity data. Though we choose
to use a constrained Ising model, this method could easily be expanded to other morphological models such as Cahn-Hilliard, phase-field models, interpenetrating sphere models, and
so on. We have also demonstrated how these sorts of morphological models may be used
to qualitatively explore the effect of various system features, which we demonstrated by
comparing morphologies generated with a simple two-phase model, a more complex threephase sequestered PCBM model (like Yin and Dadmun’s “Rivers and Streams” model) and
a crystallized P3HT fiber network model.
As noted earlier, many model morphologies are consistent with reflectometry data as
this system is underconstrained. These model morphologies are being developed for use
in transport dynamic simulations, but we can assume that any model in which morphology
plays a significant role could produce very different results for two different morphologies that

120

are fit to the same target density profile, but significantly differ in other internal features.
A clear example of this is shown in the inset to figure 5.9 which shows a sample DMC
calculation, looking at exciton dissociation ability, as a function of not only parameters used
in the model, but also morphology used. In order to further constrain the model, in chapter
6 we take advantage of small angle neutron scattering (SANS) measurements, which were
also performed on the annealed and as cast BHJ OPV devices by Jon Kiel[25].

121

(a)

(b)

(c)

(d)

Figure 5.11 Model morphologies (L = 400 × 400 × 400) fit to the unannealed (as-cast)
density profile in figure 5.10. Cross sections of the x-y plane are shown for 4 different feature
sizes, a = 4 (5.11a), a = 8 (5.11b), a = 12 (5.11c), and a = 16 (5.11d). Brown represents
PCBM sites and white represents P3HT sites.

122

(a)

(b)

(c)

(d)

Figure 5.12 Model morphologies fit to the annealed density profile seen in figure 5.10. Cross
sections of the x-y plane are shown for 4 different feature sizes(figure), a = 4 (5.12a), a = 8
(5.12b), a = 12 (5.12c), and a = 16 (5.12d). Brown represents PCBM sites and white
represents P3HT sites.

123

(a)

(b)

(c)

Figure 5.13 A comparison of electron tomography[5] performed on PPV:PCBM BHJ solar
cells to simple morphologies generated using the constrained Ising model method. First,
figure 5.13a demonstrates a zoomed in 400 nm × 400 nm section of electron tomography
taken of a BHJ morphology[5]. Next, in figure 5.13b, is an artificially colored two-phase
view of this image. Finally, figure 5.13c is a sample side-cut of an example morphology
generated with our model with comparable scale, containing feature size a = 8.1. The
interpenetrating morphology of figures 5.8b and 5.8c are evident.

124

1

Fraction PCBM

0.8

x=0
x = .1
x = .15
x = .2

0.6
0.4
0.2
0
0

100

200

3

300

400

Depth in 400 System
Figure 5.15 Assumed aggregate PCBM density profile from the three-phase model.

126

(a)

(b)

(c)

(d)

Figure 5.16 2-d cross-sections demonstrating the sequestered PCBM morphologies, with
figures 5.16a, 5.16b, 5.16c, and 5.16d representing x = 0, x = .1, x = .15, and x = .2 reduced
volume fractions respectively. Brown sites represent aggregate PCBM, white represents the
mixed PCBM/P3HT phase.

127

(a)

(b)

(c)

(d)

Figure 5.17 Superpositions of 15 cross-sections of the reduced aggregate PCBM morphologies, where only aggregate PCBM is displayed as green for sites connecting to the electron
collecting electrode and white for sites which are disconnected. The four cross sections shown
are for the x = 0, x = .1, x = .15, and x = 2. dispersed PCBM cases.

128

(a)

(b)

Figure 5.18 Figure 5.18a shows the P3HT density profiles used in model generation here,
demonstrating both the crystallized fraction (c-P3HT) and the total fraction. Figure 5.18b
shows a top down view of c-P3HT fiber network. Superposition of 20 total layers shown,
with each fiber being 4 layers deep.

(a)

(b)

Figure 5.19 Views of the unannealed c-P3HT fiber-network based morphologies. Figure
5.19a shows a 2-d cross section of the initial configuration of a model morphology, which
has had a c-P3HT fiber network inlaid but has not yet undergone annealing, such that the
PCBM (brown) and amorphous P3HT (white) are still well mixed. Figure 5.19b is a 3-d
view of the same morphology.

129

(a)

(b)

Figure 5.20 Views of the annealed (a=16) c-P3HT fiber-network based morphologies. Figure 5.20a shows a 2-d cross section of the final configuration of a model morphology, which
has undergone annealing of the amorphous P3HT and PCBM (but left the c-P3HT network
frozen) to characteristic feature size, a=16, such that the phases are now well segregated.
Figure 5.20b is a 3-d view of the same morphology.

130

Chapter 6
Transport in nanostructures
consistent with SANS data
As was highlighted in the end of chapter 5, reflectometry data on it’s own provides no information about the nature of the internal structure in-plane, only informing on an averaged
density profile in the xy-plane, at a given height through the morphology (z-dimension).
Unfortunately, this information may yield models with dramatically different performance
characteristics, as was demonstrated through the wide variety of models fitting the neutron
reflectometry data in chapter 5. This wide variation in modeled morphologies comes from
the simple fact that the system is underconstrained. Clearly, if we are to generate more
realistic morphological models for use in transport calculations we must incorporate experimental data which is complementary to the neutron reflectivity provided density profiles.
Fortunately, Jon Kiel also took small angle neutron scattering (SANS) measurements of
these P3HT/PCBM BHJ systems[8], which informs on internal structural information we
may use to further constrain our models.
SANS provides angle averaged information about the internal structure of a film in all
dimensions, rather than composition as a function of depth as was the case with reflectivity
measurements. Previous analysis of SANS measurements on P3HT/PCBM BHJ morphologies employed standard fitting procedures, which modeled the BHJ morphology as a system

131

of polydisperse spheres. Though these polydisperse sphere fits do provide information on the
nature of the nanoparticle agglomeration and highlight average changes in agglomerate size
due to system annealing, this type of analysis does not provide exact morphological models
required for use in morphologically driven device simulations. Indeed, while the use of a
polydisperse sphere model does fit the SANS data, it’s use is in no way meant to imply that
the morphology itself is composed of spherical agglomerates, as we know these BHJ contain
3-d percolative interpenetrating nanostructures.
We find that by incorporating small angle scattering data into our models, we are able
to select an appropriate characteristic feature size, a, to use for the as cast and annealed
cases. By simultaneously fitting to both the reflectivity data and the SANS data, new
morphological insights are possible in these complex systems. In this chapter, we look at
what effects the use of models consistent with both SANS and neutron reflectometry data
has on DMC device performance calculations which rely on model morphology.
Unlike the reflectometry data, we do not directly incorporate the experimental SANS
data into the Hamiltonian used to generate morphologies. Instead, we refine the model
morphologies to fit the SANS data though the use of the simulated small angle scattering
algorithm introduced in chapter 4. In this chapter, we will describe the reasons for this
procedure, as well as describing the methods used to fit to the SANS data in detail. We will
first demonstrate this procedure on the commonly used two-phase model, but then expand
to the sequestered (dispersed) PCBM model and the crystallized P3HT fiber model, all
introduced in chapter 5. In each case, after refining the models to the experimental SANS
data, relevant DMC transport calculations will be performed to understand what effects
different morphological models would have on overall device performance. The DMC model
we employ in this chapter is slightly modified from the traditional form which was presented
132

Table 6.1 Parameters extracted from SANS based on polydisperse sphere model
Morphology
PCBM agg. radius (nm)
Total PCBM volume in film
Dist. between PCBM agg. (nm)

As Cast
1.0 ± 0.75
52%
9.2

Annealed
5.9 ± 2.4
51%
16.3

in chapter 3. We will explain the reasons for the modifications, and go through our modified
DMC model in detail.

6.1

Nanostructures consistent with SANS and neutron
reflectometry data

Jon Kiel performed small angle neutron scattering (SANS) measurements on P3HT/PCBM
BHJ OPV systems on the NG-1 reactor at NIST[8], fitting the data to a set of polydisperse
spheres with a Schulz distribution to study the effects of annealing on the BHJ morphology.
The results of these fits are summarized in table 6.1, which suggest that annealing the film
caused the PCBM agglomerates to increase in radius by a factor of 6, while at the same time
the agglomerates were separating from one another.
Unlike neutron reflectometry data, SANS data is not directly invertible as there are
potentially countless morphologies that would generate seemingly identical SANS profiles.
Indeed, the original polydisperse sphere fits should not be interpreted as implying that the
nanostructure itself is composed of agglomerate spheres, only to highlight how the dominant
scattering structure (agglomerates of PCBM) increased in relative size. We generated a wide
variety of candidate morphologies by varying the Hamiltonian, simulating SANS profiles for
each. We tuned the parameters used in the Ising model until a morphology is generated
whose simulated scattering profile fits both the neutron reflectometry data and experimental
133

SANS data to high precision.
For the SANS simulation, the site weights used (the values used for b in equation 4.38)
are based on the corresponding scattering length densities (SLD) as the pure real materials,
−2

SLDP 3HT = .74 × 10−6 ˚
A

and SLDP CBM = 3.7 × 10−6 ˚
A

−2

. A low-q upturn will be

ever present in the calculated SANS profile due to the sample box size itself, which may
not directly fit the experimental data. However, as we are looking to model the internal
morphological characteristics of the BHJ thin films, and not model the entire film to scale,
this is acceptable so long as the features of interest are distinguishable from the effects of
sample box size. Previously described methods in section 4.4 of convoluting the scattered
data or reducing the box size contributions are employed, with typical smearing values of
10%.

Figure 6.1 Demonstration of how a lattice constant, blattice , can be extracted from the
Bragg-like nearest neighbor peak generated by the DFM method. Once final scaling is
performed such that the calculated scattering profile is best it to experimental data, these
lattice peaks (as highlighted here) will correspond to 2π/blattice . Shown are the fits for the
annealed and as cast two-phase morphology seen in figure 6.5.

Once model morphologies with SANS profiles comparable to those seen in experiment
134

are generated, a final fit to the experimental data is performed in which the scales of the
simulated SANS plots are adjusted until the features of interest best match the data. As we
are using a cubic lattice model, there will be a very strong Bragg-like peak due to the lattice
spacing itself seen at the q = 2π/blattice in the calculated scattering, where blattice is the
lattice constant, as is demonstrated in figure 6.1. This final scaling fit allows us to precisely
extract the lattice constant for a specific model morphology, which is useful in transport
calculations utilizing this morphology.

6.1.1

Two phase model

We first discuss the simple two-phase model, in which the nanoscale morphology is described
as a lattice of sites, each assigned either pure P3HT (donor site) or pure PCBM (acceptor
site).
The two key features of interest in the original SANS data are the amplitude and position
of the plateau-terminating bend, shown in figure 6.2. For the as cast BHJ film, this plateau
˚−1
has an amplitude of approximately 8 cm−1 , with the bend cresting at q ≈ .05A . For
the annealed BHJ film, the plateau has an amplitude of approximately 600 cm−1 and the
˚−1
bend cresting at q ≈ .015A . This plateau shift to a higher intensity and bend at a lower
q-value is indicative of an increase in characteristic feature size. As the generation of w(r)
from equation 4.38 involves the product of the SLD values of pairs, the contrast in the
intensity contribution between those features with the strongest SLD values (bS ) and the
weakest (bw ) will go as the difference of squares (b2 − b2 ). In the case of the P3HT/PCBM
w
S
system under investigation in this thesis the contrast in intensity between PCBM features
to P3HT features will be 3.72 − .742 ≈ 13, a significant signal difference. For this reason,
characteristics seen in SANS profile will be predominantly caused by PCBM features of the
135

morphology.
Simple Guinier analysis of the intensity change of the plateau informs on the aggregation
number[110], defined as the number of basic scattering units per scattering particle. In
this case, the scattering unit is taken as an individual PCBM molecule and the scattering
particle is the agglomeration of many PCBM molecules. If we assume the density of these
particles remains constant before and after annealing, then the number of scattering particles
is proportional to the volume of the scatterer (N ∝ r3 in the case of spheres). Noting that the
relative intensity of the as cast and annealed SANS plateaus are 8 and 600 cm−1 respectively,
this allows for a quick estimation of average agglomerate growth due to annealing

rannealed
=
ras cast

Nannealed 1/3
≈
Nas cast

600 1/3
≈ 4.2.
8

(6.1)

This back-of-the-envelope calculation is consistent with the results of the polydisperse
sphere fits, which showed an aggregate radius increase of 5.9 ± 5.0 (taking fractional error
propagation into effect). Though these error bars seem large, recall that a Schulz distribution
has a particularly long tail which and can be described as
R2
r·R
r σ2
exp − 2 R
σ
P (r) =
2
r · Γ( R2 )
σ

R2
σ2

R2
σ2

,

(6.2)

where P (r) is the probability of a sphere with radius r, the average sphere size is defined
by R, the width of the distribution is defined σ, and Γ is the gamma function. An example of
a polydisperse Schulz distribution is given in figure 6.3, which demonstrates the distributions
found by Jon Kiel in his fitting of the annealed and as cast morphologies. A polydisperse
Schulz sphere function is often used to fit SANS data of polymer systems[121] not because

136

the system is thought to be composed of spheres, but because this wide distribution of length
scales is thought to be similar to the length scales present in the actual systems.
For our modeling here, we use the two phase morphologies as described in section 5.3 of
this thesis. Figures 5.11 and 5.12 show examples of the as cast and annealed morphologies
we tune using our simulated SANS algorithm. Sites are assigned scattering weights based
on the pure material SLD values (found by Jon Kiel in his previous fitting methods[8]) of
bdonor = .74 and bacceptor = 3.7. We roughen the edges of the sample box by randomly
removing sites via the algorithm described in section 4.4.2.1 of this thesis, with a Gaussian
width of 120 lattice sites from the center of the box (thus covering 60% the box edge width),
as this length was found to remove a majority of high frequency sample-size oscillations
without effecting the larger scale internal features of interest.
We find best fit models to the experimental SANS data to have features sizes of ≈ 4b
for the as cast cast and ≈ 16b for the annealed case, where b is the model lattice constant.
Final scaling to the experimental data extracts lattice constants of b = .667 nm for the
as cast model and b = 1.25 nm for the annealed model. Taking the lattice constants into
effect, these results imply the annealed morphology has increased in average aggregate radius
by a factor of approximately 7.5 over the as cast morphology, which is consistent with the
standard polydisperse Schulz spheres fit, which had an average feature size growth of 5.9±5.0.
Comparison of the experimental and modeled SANS data is presented in figure 6.5.
Having selected the characteristic feature sizes to be used in these percolative, interpenetrating BHJ structure models (aas cast = 4 and aannealed = 16), we now generate a
complimentary set of morphologies with identical feature sizes but incorporating a flat target
density profile rather than one provided by the neutron reflectivity data. In this way, we
may explore the effects of the density profile on model predictions independent of feature
137

size. These so called “Ising Only” morphologies are shown in figure 6.7c and 6.7d.

6.1.2

Modified DMC Method

As we are interested in characterizing how nanoscale morphology effects device performance,
independent of the particular parameter set used in the simulation, we alter our dynamic
Monte Carlo calculations from the previously employed method presented in section 3.4 of
this thesis. Rather than specifying a specific set of model parameters for use in each calculation, such as applied voltage, electron mobility, hole mobility, work function, etc, we instead
attempt a more general characterization by sweeping over a broad range of recombination
lengths for all carriers, which we define as:

L=b

νhop 1/2
νloss

(6.3)

where L is the relevant recombination length, b is the lattice constant, νhop is the associated particles hop rate, and νloss is the associated particles recombination or decay rate.

6.1.2.1

Exciton behavior in the 2-phase model

We will first look at how morphology effects exciton dissociation efficiency when comparing
the as cast and annealed morphologies. The results of the DMC calculation are summarized
in figure 6.8. Exciton dissociation, χdiss , is defined as the ratio of excitons successfully
dissociating at a heterojunction over the total number of excitons generated. The exciton
decay length is defined Lx = b

νx-hop
νx-decay

1/2

. Excitons are generated with equal proba-

bility at all P3HT sites, and then diffuse with hop/decay rate combinations which provide
the specific decay lengths given on the axis of figure 6.8. The data presented is taken at

138

100 different decay length values, with each point being averaged over one million exciton
lifetimes. Exciton dissociation at a heterojunction is treated as exceptionally fast compared
to all potential DMC actions, such that excitons will almost always dissociate upon reaching
a site which is a heterojunction (adjacent to an acceptor site).
As the results in figure 6.8 show, the as cast morphology with it’s much finer characteristic
feature size (a ≈ 4), is much more efficient at dissociating excitons at any given exciton decay
length compared to the annealed morphology (a ≈ 16). This is a predictable outcome, as the
excitons are generated in the donor material (P3HT) and must diffuse to a heterojunction
interface in order to dissociate. Any morphology which offers shorter average paths to a
heterojunction will naturally outperform systems with longer paths, and thus smaller feature
sizes are always preferable from the standpoint of exciton dissociation efficiency. In the case
of the as cast morphology, it is observable that once the decay length is of order equal to the
feature size, essentially all excitons generated will successfully dissociate. While this effect
is not as clear in the annealed model data, this is only due to the bounds of Lx used in the
simulation and the trend is still seen.
Now, we look at the differences between the so called “realistic” morphologies, being those
generated with the reflectivity informed density profile, and the “Ising only” morphologies.
This comparison data is also shown in figure 6.8. We note that in both cases, the Ising
only morphologies exhibit better exciton dissociation than their realistic counterparts, other
than a small deviation from this trend for the a ≈ 4 morphologies for Lx-diss < 1. To
understand this low Lx-diss discrepancy, one can look at the average characteristic feature
size of each individual layer in the morphologies, shown in figure 6.9. While the Ising
only models are generally consistent at maintaining constant feature size at every depth
in the morphology, the realistic morphologies contain significant feature size fluctuations,
139

particularly in regions where the target density profile was far from 50% such as near the
edges. These regions of anisotropic feature size must exist in order to fit any non-uniform
density profile, as some regions will have greater average P3HT agglomerates (and thus,
sites to form a heterojunction) to the rest of the model, and others less. As excitons are
generated in P3HT, those regions containing lower PCBM volume fractions present, on
average, a longer path excitons must travel in order to discover a heterojunction interface.
Additionally, as excitons are generated uniformly through the bulk of the donor material
(P3HT), those regions with more total donor material will have a greater number of exciton
generation events, and thus the behavior of excitons in these regions will be slightly more
significant in the DMC averaged results. This subtle difference, entirely due to PCBM density
profile fluctuations, is what leads to a slight diminishing of exciton dissociation efficiency for
the “realistic” morphologies. For a fixed characteristic feature size, density fluctuations are
detrimental to exciton dissociation efficiency in these morphological device calculations.
At very low recombination lengths (Lx-diss

1) essentially no excitons will be able to

hop from the site they originate at prior to decay. Those excitons that do dissociate must
be generated on sites of P3HT which are already adjacent to PCBM (at the heterojunction), where due to the selected high exciton dissociation rate, successful dissociation may
occur. The crossover behavior we see in the a ≈ 4 cases at Lx ≈ .3 nm has been traced to
the small difference in BHJ surface area between the two models, as careful inspection of
each morphologies reveals that the as cast has an exact feature size of aas-cast = 4.007997,
where as the Ising only has an exact feature size of aIsing = 4.023306. As feature size is
inversely proportional to surface area, this shows that the as cast morphology has more total
P3HT sites sitting directly next to PCBM sites, which coupled with the understanding of
how density fluctuations are detrimental to overall exciton dissociation efficiency, explains
140

the crossover behavior seen in the data. For comparison, the morphology used in the annealed morphology, compared to the Ising only a=16 morphology, had exact feature sizes of
aannealed = 16.014643 and aIsing = 16.000523 respectively, which fits with our understanding as there is no such crossover for the a ≈ 16 morphologies.

6.1.2.2

Charge behavior in the 2-phase model

We now compare the morphologies from figure 6.7 in terms of their ability to collect free
charges as they diffuse throughout the morphology. The results of the simulations are presented in figure 6.10, where charge collection efficiency, χq , is defined as the ratio of free
charges collected at their respective contacts over total charges generated due to exciton
dissociation. The charge recombination length, Lq = b

νc-hop
νc-recomb

1/2

, is defined as the

lattice constant (b) times the ratio of the free charge hop rate to the free charge recombination
rate.
The results in figure 6.10 show that for all morphologies under consideration, free charge
is collected with roughly the same efficiency, with a notable exception of electron collection
efficiency in the as cast morphology, which was significantly reduced compared to all other
morphologies. This poor charge collection efficiency is due to the abrupt drop in PCBM
density which occurs at the electron collecting electrode visible on the left-hand side of the
as-cast density profile in figure 6.6. This leads to a situation where a bulk of the PCBM
sites in the morphology have no percolative paths which would allow free electrons to reach
the electron collecting electrode, regardless of the recombination length used. This is clearly
seen by comparing the PCBM sites (brown) which reach the contact layer in the top-down
views of the as cast (figure 5.14a) and annealed (figure 5.14b) morphologies.
We see a subtle (≈ 15%) drop in the maximum possible charge collection efficiency of
141

electrons in the annealed morphology compared to the Ising only morphologies, and again
find this effect is due to electrons which are generated on regions of acceptor (PCBM)
material which had no percolating paths to the electron accepting electrode.
If one compares electron collection to the hole collection efficiency (in figure 6.10) for the
as cast and annealed cases, hole collection efficiency is nearly the same in both cases, and
very similar to the electron collection efficiency for the annealed case. While this implies
that differences in morphology between the as cast and annealed morphologies did effect
hole collection efficiency, the majority of the P3HT sites here are connected to the hole
collecting contact and thus no significant disconnect in the hole current occurs. For the Ising
only data we found the electron and hole collection efficiencies to be essentially identical,
such that we present it here as a single solid line. This is understandable, as the Ising only
morphologies are isotropic at all layers, due to their flat target density profile, and thus there
is no morphological differences between the behavior of electrons and holes. Looking at the
Ising only data from Lq = 10 to 400, we see that charge transport is more efficient in the larger
features size, a = 16, compared to the finer grained a = 4 morphology. This observation is
consistent with the conventional understanding that larger aggregates of material allow for
easier diffusion of the charges, or to put it another way, morphologies with smaller feature
size offer more tortuous paths which diffusing charges must transverse. At recombination
lengths Lq > 400 the charge collection efficiency of the two different feature size models
merge together, as this is outside the regime in which the tortuous path for the diffusing
particles is the limiting factor. The only upper limits on charge collection efficiency when
Lq gets sufficiently large are therefore the limits based on charges being generated on sites
which do not connect to the associated charge’s electrode, or so called “islands”.
Overall, we see that charge collection efficiency is enhanced by greater characteristic
142

feature sizes, but that subtle effects introduced by using an experimentally measured density
profile are significant enough to alter the results of a calculation, regardless of the set of device
model parameters used.

6.1.3

Sequestered PCBM model

We now look at the effects of the so called three-phase or sequestered PCBM model on
DMC device calculations. Recall that, similar to the “Rivers and Streams” model[4], the
sequestered PCBM picture models the system as sites assigned to be either a pure phase
of aggregate PCBM or a mixed phase of P3HT and dispersed PCBM. In order to directly
compare to the previously employed annealed two-phase model, we generate models with
feature size a ≈ 16. Using the model morphologies generated in section 5.3.1, simulated
SANS profiles were generated using the distribution function method (DFM) with associated
scattering strength weights proportional to the SLD value of pure PCBM for the aggregate
PCBM phase (bP CBM = 3.7) and a weighted mixture of pure PCBM and pure P3HT
(bP 3HT = .74) for the mixed phase, the results of which are shown in table 6.2.
Table 6.2 SLD values for the four different fractions of dispersed PCBM morphologies seen
in figure 5.17, used in simulated scattering in figure 5.17.
x

bd-PCBM/P3HT

bagg-PCBM

0
.1
.15
.2

0.74
1.04
1.18
1.33

3.7
3.7
3.7
3.7

Simulated SANS calculations were performed on the morphologies shown in figures 5.16
and 5.17, with identical Gaussian fuzzying and smearing routines as those used in the previous two-phase model. The results of the SANS simulation is shown in figure 6.11.
As each of the three-phase models shown here have been generated to feature size a = 16,
143

it is not surprising that the simulated SANS profiles are quite similar. The primary difference
in the simulated scattering profiles is that as the fraction of dispersed PCBM (x) increases,
the scattering intensity decreases slightly. This can be attributed both to the decrease in
PCBM present and the decreasing contrast between the two materials (PCMB and mixed
phase) as x increases. However, all volume fractions presented generate SANS profiles which
are quite comparable to one another based on fits to the experimental SANS data, and thus
no strong conclusion can yet be drawn on which model is more accurate.
We now compare these four sequestered PCBM model morphologies through the use of
our modified DMC transport calculation, as we did previously with the two-phase model.

6.1.3.1

Exciton behavior in the sequestered PCBM model

The results of the DMC exciton transport calculations are shown in figure 6.12, with results
presented as a function of a exciton decay length, Lx = b

νx−hop
νx−decay

1/2

. These graphs were

generated using 66 different values of Lx , with each point being the average of one million
exciton lifetimes.
We note that the most efficient exciton dissociation occurs in the original x = 0 case
(no PCBM sequestered), and that progressively less excitons successfully dissociate as the
fraction of dispersed PCBM is increased. As excitons are being generated uniformly though
the bulk of the mixed phase (as photons are absorbed only in P3HT) an exciton will have,
on average, a further distance it must travel to discover a heterojunction as x is increased
because the average distance to an aggregate of PCBM increases. Note that while we already
demonstrated that increasing feature size decreases exciton dissociation efficiency, this is not
the same effect here as feature size between all four morphologies studied here is held at
a ≈ 16. Instead, this is an effect of the decreasing average surface area each donor site
144

(mixed phase) will have, as the total number of donor sites is increasing with x, but feature
size (and thus surface area) is held constant.

6.1.3.2

Electron collection in the sequestered PCBM model

Next we look at the results of DMC based electron collection efficiency calculations in figure
6.13. The most striking feature of the data is the x = .2 result, which shows electron
collection efficiency of essentially zero. This result comes from the apparent transition of the
system to a non-percolative network (in regards to reaching the electron collecting contact)
which occurs in the x = .2 morphology. It is clear that in these morphological models,
some critical threshold exists between x = .15 and x = .2, wherein the PCBM network
becomes too disconnected in order to function successfully as an OPV device. This result is
consistent with our qualitative analysis of the morphology, which showed an overwhelming
volume of disconnected PCBM sites as is demonstrated in figure 5.17 (in white). While
all other morphologies with x < .2 show relatively similar electron collection efficiency, it
is worth noting that the x = .1 case exhibits slightly improved electron collection. We
have investigated this, and found it to be an effect of well defined “rivers” of aggregate
PCBM which extend from the electron collecting contact into the bulk of the morphology.
Though one would intuitively consider more aggregate PCBM to produce a more favorable
electron transport profile, we found that due to the morphologies specific exciton dissociation
dynamics, excitons were more likely to dissociate onto sites which had a further average path
length to the electron collecting contact in the x = 0 case as compared to those excitons which
dissociated successfully in the x = .1 case. This result is directly related to the anisotropic
density profile introduced by the neutron reflectometry data, and would not have come
about without it’s inclusion. Note that hole collection efficiency is not discussed here, as the
145

principle effects of sequestering PCBM in this model will be on electron transport.

6.1.4

Crystallized P3HT fiber model

Finally, we look at morphologies generated with a crystallized P3HT (c-P3HT) template, as
were shown in figure 5.20. To directly compare to previous models, the morphologies studied
here have been Ising annealed to characteristic feature size, a = 16. Simulated SANS profiles
were calculated for the c-P3HT model using site scattering weights of bP 3HT = .74 and
bP CBM = 3.7, just as with the two-phase model. Sites assigned as c-P3HT or amorphous
P3HT are weighted the same in the simulated SANS, as the neutron SLD of the two materials
is identical. Once again, we utilize a fuzzying procedure with the same parameters as the
two-phase model, and apply a smearing function of 10% to the resultant profiles, shown in
figure 6.14.
Comparing the simulated SANS profile of the c-P3HT model to the simple two-phase
model, we see the two are very similar. The fiber template method introduced a structural
anisotropy to the system, yet as the model morphology was annealed to characteristic feature
size a = 16, it’s angle averaged scattering profile (which is sampled by our method) is nearly
indistinguishable from those structures generated without the fiber template.
DMC device performance calculations were performed on the c-P3HT model, looking
at both exciton, electron, and hole dynamics, similar to the method used to characterize
the simple two-phase morphologies. Transport through crystallized polymers is known to
exhibit significantly higher in-plane mobility, while their through-plane mobility is thought
to be significantly lower[140]. In order to examine the effect of enhanced hole mobility due to
c-P3HT regions, we also perform so-called “fast” hole transport calculations, where mobility
inside c-P3HT is enhanced in the planer (xy) direction.
146

6.1.4.1

Exciton behavior in the c-P3HT fiber model

DMC device calculations were performed on the c-P3HT model morphology and compared
to those results found for the two-phase model, the results of which are presented in figure
6.15, 66 values of Lx were modeled, with each data point being the averaged behavior of
one million exciton lifetimes. Recall that the feature size for the as cast model was a = 4,
compared to the features size of a = 16 for the annealed model, which is also the feature
size we use in the c-P3HT model. As was explained in discussions on exciton transport
dynamics in the two-phase system, this feature size difference explains why the as cast case
exhibits significantly better exciton dissociation than either of the a = 16 morphologies.
Comparing both a = 16 morphologies, however, we see the c-P3HT model had consistently
better exciton dissociation efficiencies than the two-phase annealed morphology. This can be
attributed to the introduction of a fiber network in the c-P3HT morphology, which lead to
a significant anisotropy not present in the two-phase model. In the c-P3HT morphology, a
greater fraction of the total surface area (acceptor/donor bordering sites) is stretched across
the xy-plane, rather than evenly distributed in all three dimensions.
Excitons in the c-P3HT model can thus, on average, diffusively discover heterojunctions
more quickly because of the smaller average feature size in the z-dimension. This interesting
result demonstrates that characteristic feature size alone is not sufficient to predict transport
dynamics, as any inherent anisotropic factors must also be considered in any morphologically
dependent calculation.

6.1.4.2

Charge behavior in the c-P3HT fiber model

We now perform charge transport calculations on the c-P3HT model. In addition to the
standard calculations of charge collection efficiency (χq ) versus charge recombination length
147

(Lq ), we also perform a so called “fast” transport calculation, where we examine the effect
of enhanced hole transport through crystallized P3HT regions. For these fast-transport
calculations, we treat the hole hopping rate on sites assigned c-P3HT as 100× faster than
the given Lq . As this enhanced mobility is considered only in-plane[140], hop rates through
the plane (in the z-direction) are not enhanced. The results from these calculations as well
as the two-phase model for comparison are shown in figure 6.16.
The c-P3HT model displays very similar results for electron and hole transport, with both
charges’ collection efficiency being comparable to the annealed two-phase morphology. When
examined closely, we see the measured difference between the electron and hole collection
efficiencies is in fact even smaller than the two-phase model, which is likely due to the cP3HT templating, which induced the generation of a more connected network for not only
the P3HT (donor) sites, but also the PCBM (acceptor).
The hole collection efficiency in the “fast” c-P3HT model was significantly enhanced,
which is somewhat surprising considering that hole transport dynamics through the plane (zdirection) were in no way explicitly enhanced, yet it is the transport through the plane which
must occur efficiently in order for holes to quickly reach and be collected at the contact. This
result can be understood because holes in the “fast” system may now very quickly explore
the xy-plane such that “dead-end” regions of the morphology, which previously would have
presented a significant hazard to transport, are now easily explored. This allows the holes
to more efficiently explore the morphology and quickly discover those percolative paths to
the contacts which, based on our studies of the sequestered PCBM morphologies, we know
can enhance charge collection efficiencies.
These final results are of particular interest, as they demonstrates that both exciton and
charge transport would be enhanced by a crystal fiber network from a purely morphological
148

standpoint, and even more so when considering the increased mobility granted by crystallized
materials.

6.2

Conclusions

In this chapter, we have demonstrated how the inclusion of experimentally measured data can
alter the results of device performance calculations in morphologically dependent systems,
such as the BHJ OPV systems. In our work, we utilized small angle neutron scattering data
(SANS) from PCBM/P3HT systems to refine and more precisely select a realistic model
morphology which would otherwise be underconstrained. Although we utilized two forms of
complementary neutron scattering data, the techniques presented could readily be adapted
to a variety of other experimental data such as AFM, SEM, STM, x-ray scattering, etc. We
have shown how the inclusion of experimentally measured data can have a significant effect on
device calculations when compared to predictions made with simple assumed morphologies.
The DMC method we applied here makes no assumptions about specific device parameters, instead sweeping over a wide range of potential material properties to give an overall
picture of the morphological effects on the system. Indeed, our results show that logical modifications to the simple Ising model, such as the inclusion of sequestered PCBM or a network
of crystallized P3HT fibers, can have significant effects on device performance calculations
entirely independent of the parameter set used.
While we are unable to make claims on the accuracy of one particular model over another, as every morphology used herein has been simultaneously consistent with both neutron
reflectometry and SANS data, we suggest that the methods used here (inclusion of experimental data into a modified Ising model) are an ideal framework to explore the effects of

149

many other forms of data and theoretical model systems.

150

(a)

(b)
Figure 6.2 Small angle neutron scattering (SANS) data of the P3HT/PCBM BHJ thin film
mimics, from the thesis of Jon Kiel[9]. Figure 6.2a shows the original SANS data for the as
cast and annealed systems, as well as the polydisperse sphere fits. Figure 6.2b is a cartoon
by Jon Kiel, suggesting how the agglomerate PCBM size increases might be interpreted.
151

Figure 6.3 Demonstration of Schulz Distributions of the polydisperse fits for the as cast
experimental SANS data(R = 1.0, σ = .75) and the annealed experimental SANS data
(R = 5.9, σ = 2.4).

Figure 6.4 Demonstration of how the plateau shifts to greater intensity and the bend shifts
to lower q as the constrained Ising model anneals a BHJ morphology to greater values of
a. Dashed lines are the unsmeared data and solid lines have been smeared 10%. The inset
presents the same data (smeared only), zoomed out so as to better observe the plateau.

152

Figure 6.5 Simulated SANS profiles for the as cast and annealed morphologies, showing
fits to the experimental SANS data, as well as previous fits performed with a polydisperse
sphere model. Inset in the figure are 3-d views of the as cast and annealed morphologies.

153

PCBM Volume Fraction

0.8

0.6

0.4
Experimental : As Cast
Simulated : As Cast
Experimental : Annealed
Simulated : Annealed

0.2

0
0

100

200

300

400

Distance From Air Surface (lattice units)
Figure 6.6 Density profiles of the as cast and annealed morphologies, with the simulated
morphologies compared to the experimental reflectometry fits.

154

(a)

(b)

(c)

(d)

Figure 6.7 Cross sections of the model morphologies, generated with experimentally measured density profiles, and those generated to an identical characteristic feature size, a, but
with a flat target density, which we refer to as Ising-only. The left hand column has feature
size a = 4, while the right hand column has feature size a = 16. The top row are morphologies generated with the corresponding reflectivity informed density profile, while the bottom
row is the Ising-only models, with a flat target density profile.

155

1

χdiss

0.8

Ising only (a=16)
As Cast
Annealed
Ising Only (a=4)

0.6
0.4
0.2
0

0.1

1

10

Lx(nm)
Figure 6.8 Results of exciton behavior comparing realistic morphologies to Ising only morphologies. Plot is of exciton dissociation efficiency (χdiss ) vs. exciton decay length (Lx ).
Maroon crosses represent the as cast morphology, black circles the annealed morphology,
green diamonds the Ising only of feature size as cast (a=4) and yellow boxes the Ising only
of feature size annealed (a=16).

156

30

Local feature size

As Cast (a = 4)
Annealed (a = 16)
Ising Only (a = 4)
Ising Only (a = 16)

20

10

0

100

200

300

400

Depth through morphology (lattice units)
Figure 6.9 Local feature sizes for the as cast, annealed, and Ising only (a=4,16) morphologies. Local feature size is calculated with the familiar form a = 3V /S, but V is taken as
the number of sites in a single slab V = L × L and S is taken as the total number of BHJ
interfaces seen by the sites on a specific layer.

157

1

0.1

χq

As Cast
Annealed
ising only : a=16
ising only : a=4

0.01

0.001
10

100

1000

10000

Lq (nm)
Figure 6.10 Charge collection efficiency comparison between the as cast, annealed, and
Ising only (a=4,16) morphologies. Presented as the ratio of total charges collected at their
respective contacts versus charge recombination length. For the realistic morphologies, solid
lines represent electron collection and dashed lines represent holes. For the Ising only morphologies, there was no difference between electron and hole collection efficiencies, so they
are presented as a single solid line.

158

Figure 6.11 Simulated SANS scattering for the four cases of x = 0, x = .1, x = .15, and
x = 2. fractional dispersed PCBM. The inset demonstrates the zoomed out view of the same
simulated scattering.

159

Figure 6.12 Exciton dissociation efficiency, χdiss , versus exciton decay length, Lx , for the
four cases of figure 5.17. Inset is a magnified portion of the data.

160

Figure 6.13 Electron collection efficiency, χq , versus free charge recombination length, Lq ,
for the four cases of figure 5.17. Inset is a zoomed out version of the data, on a log-log graph.

161

Figure 6.14 Simulated SANS comparisons of the simple two-phase model (black) and the
crystallized P3HT model (dashed blue), as well as the experimental data (yellow boxes) and
polydisperse sphere fits (black dots). Shown inset are views of the c-P3HT crystal network
on the left, and a 3-d view of this morphology where PCBM has been colored brown, c-P3HT
is white, and amorphous P3HT is gray.

162

Figure 6.15 Exciton dissociation efficiency for the as cast (maroon triangles), annealed
(black circles), and c-P3HT (light blue diamonds) morphologies. Exciton dissociation efficiency, χx-diss , is defined as the ratio of excitons successfully diffusing to and dissociating at
a heterojunction interface over total excitons generated. Exciton decay length, Lx , is defined
as the lattice constant times the ratio of the exciton hop rate to the exciton decay rate.

163

Figure 6.16 Charge collection efficiency for the as cast (maroon up-pointing triangles),
annealed (black circles), c-P3HT (light blue diamond) models. Solid lines represent electron
collection efficiency, while dashed lines represent hole collection. The enhanced mobility hole
transport is also shown (dark blue right-pointing triangles). Charge collection efficiency, χq ,
is defined as the ratio of charges successfully collected at their respective contacts over total
charges generated. Data presented as a function of charge recombination length, Lq .

164

Chapter 7
Conclusions
We have developed a method to efficiently generate three dimensional, interpenetrating,
percolative nanostructures which are consistent with experimental data. In this thesis, we
utilized neutron reflectometry and small angle neutron scattering data to model the bulk heterojunction nanostructure morphology as is found in many thin-film polymer based organic
photovoltaic devices, specifically the P3HT/PCBM devices generated by Jon Kiel[8, 25] of
Micheal Mackay’s group at the University of Delaware. The development of higher efficiency
thin film OPV devices is currently limited by our understanding of the influence nanoscale
morphology has on overall device performance. Though several models have been proposed
in order to explain this relationship[2, 27, 141, 136, 65, 4], none so far have extracted morphological characteristics directly from experimental data, but rather used a simple assumed
morphology. As we have shown, small changes in model morphology produce measurable
effects on calculated transport dynamics. Because of the strong influence morphology has
in these systems it is important to use the most realistic model morphology possible, as this
will lead to more accurate models and a better understanding of the system as a whole.
Previous fitting methods, using standard polydisperse sphere models[25] provide information on aggregate PCBM cluster size changes due to post-spin casting annealement of the
device, but grant no insight into the percolative interpenetrating structure which defines this
bulk heterojunction morphology and dictates many transport properties. Complex nanoscale
morphologies exhibiting such a structure are not only important in OPV, but many other
165

modern energy applications such as LEDs, batteries, fuel cells, and sensors. Each of these
fields can benefit in their modeling efforts through the use of morphologies extracted directly from experimental data, and we feel the method presented here offers a framework to
accomplish this task.
The methodology of incorporating experimentally observed features through the use of
an Ising model is quite similar to the reverse Monte Carlo technique, but distinct as we
modify the Hamiltonian of the Ising model in order to facilitate refinement, rather than
comparing directly to the scattering data. This general procedure can easily be extended to
morphological models beyond the lattice Ising model, such as a Cahn-Hillard model, phasefield models, interpenetrating sphere models, and so on. The choice of experimental data
to refine to is also much broader than the SANS and neutron reflectometry used here, as
one could easily incorporate many other forms of data such as x-ray scattering data, AFM,
STM, SEM, etc.
In this thesis, using morphologies generated with the modified Ising model, we looked
at the effects of incorporating SANS and reflectometry data into a number of different
nanoarchitecture systems, including the simple two-phase system, the sequestered PCBM
system, and a model incorporating crystallized P3HT fibers. With our modeling efforts,
we confirmed a number of well known previous postulates as well as discovering several
previously undiscussed effects of nanoscale morphology on transport in these BHJ OPV
devices. First, we confirmed the findings of previous studies[2, 18] which showed larger
characteristic feature size in nanoscale morphology lead to greater charge collection efficiency,
but also lower exciton dissociation. Second, we showed with our two-phase model that
incorporation of experimental data, here in the form of neutron reflectometry and SANS
data, has a significant effect on the results of DMC transport calculations. We also find that
166

for a specific characteristic feature size, density fluctuations throughout the morphology are
detrimental to exciton dissociation efficiency. Through our sequestered PCBM models, we
found that uniformly dropping the density of a charge carrier density (in this case, the pure
PCBM phase) will eventually lead to an abrupt drop in percolation that will effectively break
the OPV device, regardless of the effects of annealing. Finally, with the use of our c-P3HT
fiber template model we find that from a purely morphological standpoint, the introduction of
a crystal fiber network of P3HT will enhance both exciton dissociation efficiency and charge
collection efficiency. We also find that if hole mobility inside these crystal fiber networks
is enhanced only in-plane, this will dramatically increase hole transport throughout the
morphology.

167

APPENDICES

168

in figure A.1. The charge, e, is a distance z away from the interface. Inside material 1 , the
potential at a point P will be due to the charge at a distance r and the image charge, e ,
which is a distance r from point P . The potential due to this charge inside material 1 is

φ1 =

e
e
+
.
1r
1r

(A.2)

The field in material 2 will be due only to the effective charge, e , which is the charge e
screened by the dielectric boundary, such that it is in the same location but with a different
magnitude. We write the potential in material 2 as thus

φ2 =

e
2r

.

(A.3)

At the boundary between the materials, in order to fulfill Maxwell’s equations[142], the
field and it’s derivative scaled by the relative dielectric must be constant such that when
r=r

φ1 = φ2

(A.4)

and
1

∂φ1
∂φ
= 2 2,
∂n
∂n

(A.5)

where the derivatives, ∂φ/∂n are taken normal to the boundary interface. These two
boundary conditions mean that at r = r we may write

170

φ1 (r) = φ2 (r)
e
e
e
+
=
1r
1r
2r
e+e
e
=
1

(A.6)

2

and

∂φ1 (r)
∂φ (r)
= 2 2
∂n
∂n
e
∂r
e
e
∂
+
= 2
∂r
1r
1 r ∂n
2r
1

1

∂
∂r
1

−e
e
+
2
2
1r
1r

= 2

e−e =e ,

∂r
∂n

−e
2
2r
(A.7)

noting that e is on the opposite side of the boundary as e, and thus it’s derivative with
respect to the normal direction will be opposite in sign. Using the results in equation A.6
and equation A.7 we may solve for e and e to find

e =

e ( 1 − 2)
( 1 + 2)

e =

2 2e
.
1+ 2

(A.8)
(A.9)

As we are interested in a system which contains a pair of charges, we may add an

171

additional charge, q, which exists in material 2 , as is illustrated in figure A.2 The potential
on this charge would be the superposition of the field due to charge e in material 1 as is
provided in equation A.3, and the charges own image, which is similar to the second term of
equation A.2, seen as

φq (r, z) = φ2 +

q
.
2 (2z)

(A.10)

By corollary, we know that the image charge, q , must have the value

q =

q ( 2 − 1)
.
( 1 + 2)

(A.11)

Making the final assumption that we are dealing with an electron and hole, and thus
q = −e, we may write the potential experienced by charge q as

φq (r, z) =

2q
( 1 + 2)

1
r

−

( 2 − 1) q
q+ 2 2

1
2z

,

(A.12)

where we may generalize that 2 is the material the charge of interest resides in and 1
is the opposite material. The first term is clearly the attraction due to opposite charges and
the second is the image charge effect. Notice that if 1 = 2 , this reduces to the familiar
form given in equation A.1, which is what we would predict in uniform material, such as two
charges in a vacuum.
The electrostatic potential for a linear dielectric system, W , can be expressed as

W =

1
2

ρ(r)φ(r) dV

172

(A.13)

where the integral is taken over all space, ρ(r) is the charge distribution and φ(r) is the
electrostatic field. Using this, we may define the electrostatic potential of an arrangement
of two charges of opposite charge, on either side of an infinite dielectric barrier as

W =

−2q 2
q2 1 − 2
q2 2 − 1
+
+
.
( 1 + 2 )r12 4z1 1 1 + 2 4z2 2 1 + 2

(A.14)

Note that the first term of equation A.14 is the interaction energy between the two
opposite charges while the second and third terms are the interaction energy of the charges
with their respective image charges. As with the potential, we see that in the case of
1 = 2 = 1 (vacuum conditions), the static potential energy reduces to the familiar solution

of W = z1 z2 /r12

Figure A.2 Geometry of a pair of charges, with the potential on charge q given by equation
A.12

173

A.2
A.2.1

Derivation of Basic Scattering Principles
Derivation of Single Slit Scattering

We begin by deriving the results of single slit diffraction, as seen in figure A.3. Assume a
plane wave with wavelength θ is incident on a barrier, with a slit of width a. Every point
along this slit, da, will act as a point source for a spherical wave propagating in all directions,
with the assumption that all waves start out in phase from the slit. Thus, at an arbitrary
point along the detector, the contributions from the top of the slit will have to travel a
distance r1 , and the bottom of the slit will have to travel r2 . The difference of these lengths
∆r, is therefore

Figure A.3 Single slit example

174

Figure A.4 Single slit example data

∆r = r2 − r1 = a sin θ.

(A.15)

If a << θ, the slit will act as a point source, and a spherical wave will result. If λ ≈ a,
an interference pattern will occur, known as a “Fresnel diffraction” pattern if the detector is
close to the slit or a “Fraunhofer” pattern if it is far away.
The contribution to a field, E, from an arbitrary point P due to a small element of the
slit of width dy, where the edge of slit distance to P is r1 and r2 , assuming a constant source
strength at the slit of εl (units of intensity per length) is given:

ε
dE = l dy exp i (kr − ωt)
r

(A.16)

Now, consider a tiny fraction of the entire slit, dy. As seen in figure A.3, r2 is slightly

175

longer than r1 by a length of b = y sin θ, and at this point line r2 is a distance away from r1
2
of c = y cos θ. We thus can express r1 as

2
r1 = (r2 − b)2 + c2

= (r2 − y sin θ)2 + (y cos θ)2
2
= r2 + y 2 sin2 θ − 2r2 y sin θ + y 2 cos2 θ
2
= r2 + y 2 − 2r2 y sin θ
2
= r2 1 −

2y
sin θ +
r2

y2
.
2
r2
(A.17)

The limit y

r1 /2, is the Fraunhofer regime and we may drop final term of equation

A.17. We rewrite and expand this definition of r1 as

1/2
2y
sin θ
r1 = r2 1 −
r2
y
= r2 1 − sin θ + ...
r2

= r2 + y sin θ.

(A.18)

For simplicity we here define R = r1 . Using this result in equation A.16 and taking
r → R we find

176

ε
dE = l ei(kr−ωt) dy
r
ε
= l ei(k(R−y sin θ)−ωt) dy
R
εl i(kR−ωt) −i(ky sin θ)
= e
e
dy.
R

(A.19)

We can then use A.19 to find the effect of all such small sections of the total slit, a, by
integrating from −a/2 to a/2 to find E as

a/2
ε
E = l ei(kR−ωt)
e−i(ky sin θ) dy
R
−a/2
a/2

εl i(kR−ωt) e−iky sin θ
= e
R
−ik sin θ
−a/2


−i ka sin θ
i ka sin θ
ε
e 2
−e 2

= l ei(kR−ωt) 
R
−ik sin θ
ka
εl i(kR−ωt) −2i sin 2 sin θ
= e
R
−ik sin θ
ka
εl a i(kR−ωt) sin 2 sin θ
E=
e
.
ka sin θ
R

(A.20)

2

For convenience, we define β = ka sin θ, so that equation A.20 has the form
2

ε a sin (β) i(kR−ωt)
E= l
e
R
β
εa
= l sinc (β) ei(kR−ωt) .
R

177

(A.21)

We can define relative intensity of the delivered wave, Irel , as

Irel = |E|2 = I0 sinc2 β,
where I0 is the maximum intensity given by

(A.22)

εl a 2
. The sinc function behaves similarly
R

to a sin function, but decays inversely with |β|, such that Irel has an absolute maximum
at β = θ = 0. The descending maxima will occur at points satisfying β = tan β and the
minima will occur at points satisfying β = nπ for all non-zero integer values of n.
At small angles of θ, distances between maxima becomes δθ = λ . If one increases the slit
a
width until a

λ, the distance between peaks at small θ tends toward 0, and a shadowing

effect becomes dominant, rather than diffraction.
The general trends, over the range of θ from −π/2 to π/2 are that zeros only appear
once ka ≥ 2π, the peak at θ = 0 is always I0 , and as ka increases, the width of this peak
decreases. as new maxima appear with separation ∆θ = λ/a.
An alternative, and arguably simpler calculation of this scattering pattern is through the
use of the Fourier transform. If one defines the single wide slit as an aperture, A(x) as


 1
A(x)wide slit =

 0

if |x| ≤ a/2
otherwise,

and after a quick change of variables where q = 2π sin θ/λ = 2πk sin θ we find

178

(A.23)

∞

E(x) = E0
= E0

A(x)eiqx

−∞
a/2

−a/2

eiqx

E0 iqa/2
e
− e−iqa/2
iq
E
qa
= 0 2 sin
iq
2
qa
2E0 sin 2
=
,
q
=

(A.24)

where we have once again recovered the sinc behavior seen in equation A.21.

A.2.2

Derivation of Double Slit Scattering

Consider the case of two slits, spaced a distance d apart, with slit width of a, with a plane
wave incident on one side as demonstrated in figure A.5. The field propagated to the other
side of the screen will be the sum of the two individual slits, as derived previously

Figure A.5 Double slit example

179

εa
εa
E = l sinc (β) ei(kR1 −ωt) + l sinc (β) ei(kR2 −ωt) ,
R1
R2

(A.25)

where R1 and R2 are the distances from the two slits to any point on the detector. For
convenience, we define the paths from the two slits to the point P in terms of their separation,
d, and the angle θ, as R1 = R and R2 = R1 − d sin θ = R − d sin θ. As such, we may now
rewrite A.25 as

εa
εl a
E = l sinc (β) ei(kR−ωt) +
sinc (β) ei(kR−kd sin θ−ωt) .
R
R − d sin θ
First, note that we assume R

(A.26)

d, as we are considering the situation where distance

to the detector is much larger than the distance between the slits. As such, we reduce
the denominator of the second term from (R − d sin θ) to (R). Next, we define a variable
similar to β = ka sin θ from the single slit, but incorporating the distance between slits as
2
α = kd sin θ . We incorporate these changes into equation A.26 and reduce it to the form
2

εa
εa
E = l sinc (β) ei(kR−ωt) + l sinc (β) ei(kR−2α−ωt)
R
R
εa
= l sinc (β) ei(kR−ωt) 1 + e−i2α
R
recalling that: eiα + e−iα = 2 cos α
εa
= l sinc (β) ei(kR−ωt) 1 + e−i2α
R
εa
= l sinc (β) (2 cos α) ei(kR−ωt)
R

re−iα eiα + e−iα
eiα + e−iα

εa
E = l sinc (β) (2 cos α) ei(kR−α−ωt) .
R

180

2 cos α
eiα + e−iα

(A.27)

Notice that this final form, equation A.27, is similar to the single slit diffraction pattern,
equation A.21, but has an additional factor of (2 cos α), and a phase shift in the harmonic
function. Expressing relative intensity again through the use of A.22, we find the double slit
intensity to be

Ids = 4I0 cos2 (α) sinc2 (β) .

(A.28)

Notice that the intensity of the double slits is the product of the intensity of the infinitely
thin pair of slits (which is proportional to cos2 α) and the intensity of a single wide slit (which
goes as sinc2 β), a direct consequence of the Fourier relationship properties of a convolution
as described previously in equation 4.5 as

f (x) = g(x) ⊗ h(x) ⇔ F (k) =

√

2π G(k) × H(k),

where F (k), G(k), and H(k) are the resultant Fourier transformed functions of f (x), g(x),
and h(x) respectively. As the convolution of a pair of infinitely thin slits and a single finite
“thick” slit results in two finite “thick” slits, the resultant scattering of two finite slits is the
product of the scattering from an an infinitely thin pair and a single finite slit.
As the slits are brought closer, d → 0, the factor α → 0, and the single slit case is
recovered. As the slits become more narrow, a → 0, the factor β → 0, and the behavior of
the traditional Young’s double experiment emerges, in which the distance between maximas
become δθ = nθ, with n = 0, 1, 2, 3, ... and so on.
As with the single wide slit case, we may alternatively compute the scattering by defining
an aperture function with two points distance d apart

181

A(x) = δx−d/2 + δx+d/2

(A.29)

the scattering function is then a Fourier transforms

∞

E(q) = Eo

−∞

δx−d/2 + δx+d/2 eiqx dx

= Eo eikd/2 + e−ikd/2
= Eo 2 cos

dq
2

,

(A.30)

which clearly reconstructs the cosine behavior in the wave, and thus cosine squared
behavior in the intensity, as seen in equation A.28.

A.2.3

Origin of Nyquist Limits

Before discussing the various Nyquist limits placed on a Fourier transformation based scattering calculation, it is worth going through a simple example at the limits of analysis, a
one dimensional delta function in real space f(x), and it’s compliment in momentum space,
F(k), and how this relates to the well known limits of the Heisenberg uncertainty principle.
We begin with a straightforward Fourier transform of a delta function

182

∞

F (k) = F.T. [f (x)] = (2π)−1/2

−∞
∞

= (2π)−1/2

−∞

dx e−ikx f (x)
dx e−ikx δ(x)

= (2π)−1/2 ,

(A.31)

which is the expected result of a constant at all points in k-space, showing that a signal
infinitely thin is composed of an infinite number of different Fourier contributions, seen as
an infinitely wide signal in Fourier space. To get around this impasse in analysis, we model
the delta function as a Gaussian, P (x), defined

2
2
P (x) = Ax e−x /(2σx ) ,

(A.32)

where the width is σx , we derive the normalization constant, Ax as

∞

1=

−∞

∞

dx P (x) =

−∞

2
2
dx Ax e−x /(2σx )

√
= Ax σx 2π,

clearly showing that in order to normalize the probability, the prefactor must be

Ax =

1
√

σx 2π

.

Recalling that a wavefunction is related to a probability distribution as

183

(A.33)

Px = |ψ(x)|2 ,

(A.34)

we may thus define ψ(x) as a square root of the described normal distribution

ψ(x) =

2
2
Ax e−x /4σx .

(Px ) =

(A.35)

Just as this wavefunction is described in real space as ψ(x), so can it be described in
momentum space as φ(k). Making the assumption that just as the position can be described
as a normal distribution, P (x), so can the momentum be assumed to fit such a distribution
in k-space

φ(k) =

(P (k)) =

Ak e

2
−k 2 /4σk

(A.36)

where
Ak =

1
√

σk 2π

.

(A.37)

As ψ(x) and φ(k) are Fourier conjugates of one another and both are even functions we
may relate though through the cosine Fourier transform as

ψ(x) =

2 ∞
dk φ(k) cos (kx)
π 0

(A.38)

with the reverse transform being
φ(k) =

2 ∞
dx ψ(x) cos (kx) .
π 0

(A.39)

Using the definition of φ(k) from equation A.36 and ψ(x) from equation A.35, we can
184

rewrite equation A.38 as

2 ∞
dk
π 0

2
2
Ax e−x /4σx =

Ak e

2
−k 2 /4σk

cos (kx) .

(A.40)

We can solve the integral on the right with the use of the solution

∞

2 2
dx e−a x cos (bx) =

0

√

π −b2 /4a2
e
,
2a

(A.41)

such that equation A.40 becomes

√
2
2
π
−x2 σk )
Ak
e
π
2(1/2σk )
√
2
−x2 σk
= σk 2 Ak e
.

2
2
Ax e−x /4σx =

(A.42)

In order for this equality to hold, one finds that the relationship between σx and σk must
be

2 2
σx σk =

1
4

(A.43)

or more simply
1
σx σk = .
2

(A.44)

Note this equality exists only because this distribution of momentum and energy minimizes the uncertainty in both, as any function whose Fourier transform did not recover
itself would inherently have a different width (σ), and the product for a general function

185

could be larger than 1/2. The value of 1/2 on the RHS is a direct consequence of our use
of the square-root of the probability, as this is how wavefunctions are described in quantum
mechanics. If one were to instead perform an identical calculation using P (x) and P (k),
one would find the product of the widths of those functions to be 1. Finally, through this
discussion we have referred to σ as the width of a Gaussian distribution, but it could equally
be considered an uncertainty of a measurement, where the center of the distribution is the
measurement. Taking both these into account, we can write the relation

1
∆x∆k ≥ .
2
Using the relationship between momentum and wavenumber, p =

(A.45)
k, where

is the

reduced Planck’s constant, one recovers the familiar form of the Heisenberg uncertainty
principle

∆x∆p ≥

2

.

(A.46)

Although it is entertaining to see that the Heisenberg uncertainty principle can be derived
with a few simple assumptions and the Fourier transform, a more general concept is that
the uncertainty (width) of one variable is inversely related to any the uncertainty in that
function’s Fourier conjugates. This is the origin of the concept of a Nyquist interval, drN ,
which states

drN =

1
2π
=
,
2fB
2ωB

(A.47)

where B is the highest possible frequency of a system, and we have utilized the relation-

186

ship between ordinary frequency (f ) and angular frequency (ω) of f = 2π/ω. In the case
of the relationship between a PDF(r) and the scattering, I(q), which are Fourier conjugates
of one another in real (r) and wavevector (k) space. For example, the uncertainty of the
pair-distribution function[143] is inherent in the resolution in real space, here the Nyquist
interval drN , but the limiting frequency, B, is more easily understood if we consider the
entire scattering profile to be symmetric about q = 0, such that it then resembles a probability distribution with a width of qmax . As such, we recover the commonly used relationship
between resolution in the PDF(r) and qmax .

dr =

A.3

π
qmax

.

(A.48)

Detailed examination of the First Reaction Method
and Diffusion Rates

In order to verify the validity of the Walker/Greenham approach to charge transport in
Polymer-Fullerene based photovoltaics via the first reaction method (FRM) in this appendix
we examine the action of charge movement. Using FRM, particle movement is a process in
which all potential actions for a particle have an associated time, τ , generated for them:

τi = −

1
log (X)
ωi

(A.49)

Where ωi is the characteristic rate for an event and X is a random number between 0
and 1. As this transport model occurs on a cubic lattice, there can be as many as 6 different
“hop” events for a particle; potentially less if a specific geometry bars transport in some of

187

the directions prior to calculation of a hop time in that direction.
In previous works[2, 18] the calculation of the characteristic hop rates associated with a
specific mobility, µ, derive from the Einstein relation (equation A.50), and the relationship
between hop frequency and diffusion constants (equation A.51) which we will derive later.
These two equations are:

q
µ
kb T

(A.50)

a2 = 6D/ωi

(A.51)

D=

where D is the diffusion constant, q is the charge, kB T is the thermal energy, µ is
mobility, and a is the lattice constant (or characteristic hop length, which can only be the
lattice constant in this model). By combining these equations, one generates a hop rate
based on a given mobility:

ωi =

A.3.1

6kB T µ
qa2

(A.52)

Simple Diffusion

Before even approaching the validity of a charge hop rate between sites with varying energies
(something integral to both Walker’s and Greenham’s transport model) let us look at the
behavior of a simple diffusion in continuum compared to cubic lattice space. To study this
diffusion behavior, we employed a simulation in which a diffusing particle is generated at
origin, Ro = (0, 0, 0), and then hops about the system randomly. By comparison, in the
case of a cubic lattice, a hop is randomly chosen which moves the particle in one of the six
188

possible directions (±x, ±y, ±z). For continuous diffusion, the diffuser will always hop with
a distance of 1, but is in no way confined to cubic lattice directions, i.e. we allowed the
diffuser to exist at non-integer coordinates and move in any direction. To facilitate uniform
sampling on a sphere I select two random numbers: u ∈ [−1, 1] and φ ∈ [0, 2π]. These can
then be translated to Cartesian coordinates as

x =

1 − u2 cos φ

y =

1 − u2 sin φ

z = u.

(A.53)

To select a time for a hop we specify a general hop rate, ωhop , which we then either explicitly use for every hop in diffusion (producing constant hop times) or employ in our FRM
calculation through the use of equation (A.49).

A.3.2

Measuring the Effective Diffusion Constant

As we are interested in measuring the effective mobility, µeff , we here derive the relationship
between µeff and the average distance a diffusing particle travels from the origin over time.
The average displacement of a 3-dimensional cubic lattice bound random walker, after
a single hop that is equally likely to move in any of the 6 directions goes as the average of
all 6 possible lattice moves. As such, the average radial position after a single hop will be
2
r1 = 0 and the average variance after one hop will be r1 = a2 , where a is the lattice

constant. After N -hops, we can clearly see that the average displacement is going to remain
at the origin, but the average variance will grow with the number of hops such that
189

rN

2
rN

= N ri

2
= N ri

= 0

(A.54)

= N a2 .

(A.55)

As the number of hops grows large, a site’s occupation probability becomes a Gaussian
distribution about origin, with a variance given above (and normalized) as

(x − x0 )2
Px (x) = √
exp −
.
2σ 2
2πσ 2
1

(A.56)

Putting together similar expressions for Py and Pz , we can define P (r):

Pr (r) = Px (x) × Py (y) × Pz (z)

Pr (x, y, z) =
=

3· 3
2
2a3 N 3/2 π 3/2
2πN a2
3

exp −

−3/2

exp −

(A.57)

3 x2 + y 2 + z 2
2a2 N

3 r2
2a2 N

(A.58)

Note that the probability distribution is isotropic (having no angular dependence), as is
expected for a diffusing particle with no directional bias. The result of integration for the
displacement is

190

2π

r

π

=
0

∞
−∞

0

r

2πN a2
3

−3/2

exp −

3 r2
r2 sin θ dr dθ dφ
2a2 N

= 0

(A.59)

and the variance is

r2

2π

∞

π

=
0

=

0

2πN a2
3

r2

0

−3/2

3 r2
exp − 2
r2 sin θ dr dθ dφ
2a N

N a2 .

(A.60)

So we see that the variance of a diffusing particle goes as the product of the number of
hops and the average hop length squared. If particle hopping occurs at a constant rate of ω,
then after a time t has passed, we can expect the number of hops taken to be N = ωt. We
express the variance of such a particle at time t to be:

r2

= N a2 = ωta2

or if solving the 1D form :
x2

= N a2 = ωta2

(A.61)

Next we relate the diffusion constant itself, D, to the simulation measured variance, r2 .
First we derive the relationship between the diffusion constant and the hopping rate. We
start with the diffusion equation itself in 1-D:

191

∂Px
∂ 2 Px
= D
∂t
∂x2

(A.62)

If we plug the results of equ.(A.61) into equ.(A.56), and solve for D in equ.(A.62), one
finds that in 1-D:

D=

a2 ω
⇒ or equivalently ⇒ a2 ω = 2D
2

(A.63)

such that when we then plug these results back into the 1D results of equ(A.61), we find:

x2 = 2Dt

(A.64)

To get higher dimensional forms of this result, we must simply add together the variances
such that the result in 3D is

r2 = 6Dt

(A.65)

Now, we will use the equations derived and simple diffusion simulation described to
examine the results of varying an explicitly given hop rate, ωhop . We measure the variance
squared over time, knowing from equ.(A.65) that this should produce a value of 6D, the
results of the calculation are shown in figure A.6. As we see, for both continuous diffusing
(in all directions) and cubic diffusing (on a cubic lattice) the diffusion constant can easily
and accurately be extracted. It is noteworthy that both the cubic and continuous diffusion
average generate nearly identical values, illustrating that cubic hopping has no noteworthy
effect on the results. Thus, the results of diffusion transport calculations in a cubic lattice
model can be correctly assumed to be similar to the results for a real space model.
192

Figure A.6 A simple random walker simulation with hops either confined to a cubic lattice
or allowed in all directions but with fixed hop distance. Averaged over 107 walkers. On the
LHS is the average modulus of the radius traveled vs. time. On the RHS is the average
variance over time vs time, which according to equ.(A.65), will give 6D.

A.3.3

Correcting for attempt number

One problem we find with the FRM method as it has been employed previously is that
no consideration has been taken for how the number of hop-actions possible will influence
transport. For example, a diffuser in a cubic lattice which is in a corner will only have 3
potential movement actions considered in the FRM method, where as a diffuser in the bulk of
the material will have 6 potential hop actions. Based on the semi-random nature of the FRM
method, the more chances an action has to occur (in this case diffuser hopping), the faster
an averaged selected event will be. We find this can be corrected by inversely scaling the
average hop rate used to generate an event time (equation A.49) by the number of attempts
a current action is allowed (i.e. 1/6 in the bulk and 1/3 in a corner), with the form

ωeff =

ωhop
,
Ntries

193

(A.66)

where ωeff is the effective hop rate to use in the equation A.49 to generate an average
hop rate of ωhop when used over Ntries number of attempts.
To illustrate this problem and correction, we measure the value of 6D (being proportional
to the variance over time) for a series of diffusers with hop rates generated using hop rates of
ω = 1 and ω = 1/6 with the FRM method being employed in the standard way (6 attempts
on a cubic lattice). We also show the results in which only ω = 1 is employed, but only a
single time is generated with the hop vector in the cubic lattice being randomly selected.
The results are shown in figure A.7.

Figure A.7 A demonstration of how desired hop rate is affected by attempt rate in the
FRM method.

As one can see, the measured diffusion constant is a factor of 6 too high for the standard
FRM method implementation, resulting in a value of 6D=6 rather than the desired 6D=1
(seen in using A.51 with a ωi = 1 and a = 1.). The alternative method of randomly
selecting a hop vector and calculating a time recovers the desired results, as does our proposed
correction. This result shows that the previously uncorrected FRM method will on average
194

select a faster hop rate than desired.

A.4

Relationships between various correlation functions
and the scattering profiles

Here, we outline the relationships between various frequently used functions defining structure and their scattering corollaries. Common variables used here are scattering power for a
scatterer, fi , scalar distance between two scatterers (i and j), rij , total number of scatterers
in a sample, N , number density of scatterers, ρ0 . Note that although scattering power inherently has a q-dependence, in the small angle region scattering power may be considered
constant[111].

A.4.1

Structural correlation functions definitions

Radial distribution function, R(r), is given

R(r) =

1
f

fi∗ fj δ(r − rij ).
i

(A.67)

j

Real space pair density, ρ(r), is given

ρ(r) =

1

fi fj δ(r − rij ).
4πr2 N f 2 i=j

(A.68)

Atomic pair distribution function (PDF), g(r), is given

g(r) = 4πr [ρ(r) − ρ0 ] =

1

fi fj δ(r − rij ) − 4πrρ0 .
rN f 2 i=j
195

(A.69)

Reduced pair distribution function, G(r), is given

G(r) = 4πrρ0 [g(r) − 1] .

A.4.2

(A.70)

Scattering profile function definitions

The full scattering, I(q), is given

I(q) = Icoh (q) + Iincoh (q),

(A.71)

where the incoherent scattering, Iincoh , is defined as,

fi∗ fi = N f 2 ,

(A.72)

fi∗ fj exp i q · rij .

(A.73)

Iincoh =
i

and the coherent scattering, Icoh is,

Icoh =
i,j

The discrete scattering, Id , which is the coherent component minus the incoherent, may
be defined

fi∗ fj exp i q · rij .

Id = Icoh − Iincho =

(A.74)

i=j

The total scattering structure function, S(q), is defined as

S(q) =

Icoh

−
N f 2

196

(f − f )2
f 2

,

(A.75)

which if we take into account angle-averaging, may be written as

S(q) − 1 =

1
N f

fi∗ fj exp i q · rij .

2

(A.76)

i=j

The reduced total scattering structure function, F (q), being [S(q) − 1] rescaled by q and
averaged over all angles, may be written

F (q) = q [S(q) − 1] =

1
N f

fi∗ fj

2
i=j

sin q rij
.
rij

(A.77)

The sine-transform of F (q) generates a function, k(r), which is defined

k(r) =

=

sin q rij
1
2 ∞
fi∗ fj
sin(qr)dq
π 0 N f 2
rij
i=j
1
rN f

fi∗ fj δ(r − rij ) − δ(r + rij ) ,

2
i=j

limiting the calculation to values of r > 0 only, one finds

k(r) =

A.4.3

1
rN f

fi∗ fj δ(r − rij ).

2

(A.78)

i=j

Correlation relationships

The relationships between the total scattering structure function S(q), total scattering intensity I(q), reduced total scattering structure function F (q), and the coherent scattering
intensity, Icoh (q), are

197

S(q) =

I(q)

=
f 2

(f − f )2
F (q)
I
+ 1 = coh 2 −
.
q
N f
f 2

(A.79)

The relationships between the real space pair density, ρ(r), the radial distribution function R(r), and Fourier transform of the reduced total scattering structure function, k(r), the
atomic pair distribution function, g(r), and the reduced pair distribution function G(r), are

ρ(r) =

A.4.4

k(r)
R(r)
g(r)
G(r) + 4πrρ0 (1 + 4πrρ0 )
=
=
+ ρ0 =
.
2
4πr
4πr
4πr
16π 2 r2 ρ0

(A.80)

Relationship to the weighted pair density function

Our definition of the weighted pair density function, w(r), which is only sampling a random
number of pairs, Np , is
Np

w(r) =

fi∗ fj δ(r − rij ),

(A.81)

i=j

which, as the number of pairs sampled tends toward the total number of pairs in the
system, asymptotically approaches a rescaled ρ(r) as

lim [w(r)] = 4πr2 N f 2 ρ(r).
Np →all

(A.82)

The effective scattering profile we calculate, I a (q), is

Rmax

I a (q) = C

w(r)
Rmin

sin (q r)
,
qr

(A.83)

where the sum is taken over r from the smallest to largest lengths used in w(r), and C
is a normalization constant used to fit simulations to experiment. As we have calculated
198

I a (q), it’s closest scattering corollary is the coherent scattering, Icoh (q), which has been
angle averaged over the scattering interaction q · rij , though as one can see in equation A.79,
this quantity is closely related to many others.
Note also that we do not use actual scattering powers, f , as these are typically measured
in units of 10−6 ˚
A

−2

, instead using scaled scattering powers, b, with typical values of 0, 1, 5,

etc. If direct comparison to experiment is desired, correct unit scaling may be accounted for
in the normalization constant, C.

199

BIBLIOGRAPHY

200

BIBLIOGRAPHY
[1] F. Yang and S. R. Forrest, “Photocurrent generation in nanostructured organic solar
cells,” ACS Nano, vol. 2, no. 5, pp. 1022–1032, 2008.
[2] P. K. Watkins, A. B. Walker, and G. L. B. Verschoor, “Dynamical monte carlo modelling of organic solar cells: The dependence of internal quantum efficiency on morphology,” Nano Letters, vol. 5, no. 9, pp. 1814–1818, 2005.
[3] W. Chen, T. Xu, F. He, W. Wang, C. Wang, J. Strzalka, Y. Liu, J. Wen, D. J. Miller,
J. Chen, K. Hong, L. Yu, and S. B. Darling, “Hierarchical nanomorphologies promote
exciton dissociation in polymer/fullerene bulk heterojunction solar cells,” Nano Letters,
vol. 11, no. 9, pp. 3707–3713, 2011.
[4] W. Yin and M. Dadmun, “A new model for the morphology of P3HT/PCBM organic
photovoltaics from small-angle neutron scattering: Rivers and Streams,” ACS Nano,
vol. 5, no. 6, pp. 4756–4768, 2011.
[5] S. van Bavel, E. Sourty, S. Veenstra, J. Loos, et al., “Three-dimensional nanoscale
organization of polymer solar cells,” J. Mater. Chem., vol. 19, no. 30, pp. 5388–5393,
2009.
[6] S. D. Oosterhout, M. M. Wienk, S. S. van Bavel, R. Thiedmann, L. J. A. Koster,
J. Gilot, J. Loos, V. Schmidt, and R. A. J. Janssen, “The effect of three-dimensional
morphology on the efficiency of hybrid polymer solar cells,” Nature Materials, vol. 8,
no. 10, pp. 818–824, 2009.
[7] D. R. Kozub, K. Vakhshouri, L. M. Orme, C. Wang, A. Hexemer, and E. D. Gomez,
“Polymer crystallization of partially miscible polythiophene/fullerene mixtures controls morphology,” Macromolecules, vol. 44, no. 14, pp. 5722–5726, 2011.
[8] J. W. Kiel, B. J. Kirby, C. F. Majkrzak, B. B. Maranville, and M. E. Mackay, “Nanoparticle concentration profile in polymer-based solar cells,” Soft Matter, vol. 6, no. 3,
pp. 641–646, 2010.
[9] J. W. Kiel, NANOPARTICLE ASSEMBLY IN POLYMER BASED SOLAR CELLS.
PhD thesis, Michigan State University, 2010.

201

[10] D. Anderson, Clean Electricity from Photovoltaics. London: Imperial College Press,
2001.
[11] R. Gottschalg, The Solar Resource and the Fundamentals of Radiation for Renewable
Energy Systems. Sci-Notes, Oxford, 2001.
[12] S. H. Park, A. Roy, S. Beaupre, S. Cho, N. Coates, J. S. Moon, D. Moses, M. Leclerc,
K. Lee, and A. J. Heeger, “Bulk heterojunction solar cells with internal quantum
efficiency approaching 100vol. 3, pp. 297–302, May 2009.
[13] S. Forrest, “The limits to organic photovoltaic cell efficiency,” MRS Bulletin, vol. 30,
pp. 28–32, 2005.
[14] M. A. Green, K. Emery, Y. Hishikawa, W. Warta, and E. D. Dunlop, “Solar cell
efficiency tables (version 39),” Progress in Photovoltaics, vol. 20, pp. 12–20, JAN 2012.
[15] S. R. Forrest, “The path to ubiquitous and low-cost organic electronic appliances on
plastic,” Nature, vol. 428, pp. 911–918, Apr. 2004.
[16] R. Steim, S. A. Choulis, P. Schilinsky, and C. J. Brabec, “Interface modification for
highly efficient organic photovoltaics,” Applied Physics Letters, vol. 92, no. 9, p. 093303,
2008.
[17] R. Shikler, M. Chiesa, and R. H. Friend, “Photovoltaic performance and morphology
of polyfluorene blends:the influence of phase separation evolution,” Macromolecules,
vol. 39, no. 16, pp. 5393–5399, 2006.
[18] N. C. Greenham and M. Gratzel, “Nanostructured solar cells,” Nanotechnology, vol. 19,
p. 420201, Oct. 2008.
[19] M. A. Brady, G. M. Su, and M. L. Chabinyc, “Recent progress in the morphology of
bulk heterojunction photovoltaics,” Soft Matter, vol. 7, no. 23, pp. 11065–11077, 2011.
[20] D. C. Coffey, O. G. Reid, D. B. Rodovsky, G. P. Bartholomew, and D. S. Ginger,
“Mapping local photocurrents in polymer/fullerene solar cells with photoconductive
atomic force microscopy,” Nano Letters, vol. 7, no. 3, pp. 738–744, 2007.
[21] O. G. Reid, K. Munechika, and D. S. Ginger, “Space charge limited current measurements on conjugated polymer films using conductive atomic force microscopy,” Nano
Letters, vol. 8, no. 6, pp. 1602–1609, 2008.

202

[22] K. Maturova, R. A. J. Janssen, and M. Kemerink, “Connecting scanning tunneling
spectroscopy to device performance for polymer: Fullerene organic solar cells,” ACS
Nano, vol. 4, no. 3, pp. 1385–1392, 2010.
[23] L.-M. Chen, Z. Hong, G. Li, and Y. Yang, “Recent progress in polymer solar cells:
Manipulation of polymer:fullerene morphology and the formation of efficient inverted
polymer solar cells,” Advanced Materials, vol. 21, no. 14-15, pp. 1434–1449, 2009.
[24] E. J. Spadafora, R. Demadrille, B. Ratier, and B. Grevin, “Imaging the carrier photogeneration in nanoscale phase segregated organic heterojunctions by kelvin probe force
microscopy,” Nano Letters, vol. 10, no. 9, pp. 3337–3342, 2010.
[25] J. W. Kiel, A. P. R. Eberle, and M. E. Mackay, “Nanoparticle agglomeration in
polymer-based solar cells,” Physical Review Letters, vol. 105, no. 16, p. 168701, 2010.
[26] K. Schmidt-Rohr, “Simulation of small-angle scattering curves by numerical fourier
transformation,” J. Appl. Cryst., vol. 40, no. 1, pp. 16–25, 2007.
[27] C. Groves, L. J. A. Koster, and N. C. Greenham, “The effect of morphology upon
mobility: Implications for bulk heterojunction solar cells with nonuniform blend morphology,” Journal of Applied Physics, vol. 105, no. 9, p. 094510, 2009.
[28] G. A. Buxton and N. Clarke, “Computer simulation of polymer solar cells,” Modelling
and Simulation in Materials Science and Engineering, vol. 15, no. 2, pp. 13–26, 2007.
[29] O. Stenzel, L. J. A. Koster, R. Thiedmann, S. D. Oosterhout, R. A. J. Janssen, and
V. Schmidt, “A new approach to model-based simulation of disordered polymer blend
solar cells,” Advanced Functional Materials, vol. 22, no. 6, pp. 1236–1244, 2012.
[30] M. R. Hammond, R. J. Kline, A. A. Herzing, L. J. Richter, D. S. Germack, H. W. Ro,
C. L. Soles, D. A. Fischer, T. Xu, L. P. Yu, M. F. Toney, and D. M. DeLongchamp,
“Molecular order in high-efficiency polymer/fullerene bulk heterojunction solar cells,”
ACS Nano, vol. 5, no. 10, pp. 8248–8257, 2011.
[31] D. Olds, P. Duxbury, J. Kiel, and M. Mackay, “Percolating bulk heterostructures from
neutron reflectometry and small-angle scattering data,” Physical Review E, vol. 86,
no. 6, p. 061803, 2012.
[32] P. W¨rfel, Physics of Solar Cells: From Principles to New Concepts. Wiley-VCH,
u
2005.

203

[33] L. C. Tom Markvart, Practical Handbook of Photovoltaics: Fundamentals and Applications. Elsevier Science Inc, 2003.
[34] S. Fonash, Solar Cell Device Physics. New York: Academic, 1980.
[35] J. Nelson, The Physics of Solar Cells. Imperial College Press, 2003.
[36] D. M. S. C. D. F. D. J. Ellison, J. Y. Kim, “Determination of quase-fermi levels across
illuminated organic donor/acceptor heterojunctions by kelvin prove force microscopy,”
Journal of the American Chemical Society, vol. 133, pp. 13802–13805, 2011.
[37] U. Wolf, V. I. Arkhipov, and H. B¨ssler, “Current injection from a metal to a disordered
a
hopping system. i. monte carlo simulation,” Phys. Rev. B, vol. 59, pp. 7507–7513, Mar
1999.
[38] J. Bardeen, “Surface states,” Phys. Rev., vol. 71, p. 717, 1947.
[39] K. N. S.M. Sze, Physics of Semiconductor Devices. Wiley-Interscience, 2006.
[40] P. Landsberg, Recombination in Semiconductors. New York: Cambridge University
Press, 1991.
[41] C. R. McNeill and N. C. Greenham, “Charge transport dynamics of polymer solar cells
under operating conditions: Influence of trap filling,” Applied Physics Letters, vol. 93,
p. 203310, Nov. 2008.
[42] C. M. Ramsdale, J. A. Barker, A. C. Arias, J. D. MacKenzie, R. H. Friend, and N. C.
Greenham, “The origin of the open-circuit voltage in polyfluorene-based photovoltaic
devices,” Journal of Applied Physics, vol. 92, pp. 4266–4270, Oct. 2002.
[43] J. Nelson, “Diffusion-limited recombination in polymer-fullerene blends and its influence on photocurrent collection,” Phys. Rev. B, vol. 67, p. 155209, 2003.
[44] R. F. J.J.M. Halls, Organic photovoltaic devices. Imperial College Press, 2001.
[45] S. G¨es, H. Neugebauer, and N. S. Sariciftci, “Conjugated polymer-based organic solar
n
cells,” Chemical Reviews, vol. 107, no. 4, pp. 1324–1338, 2007. PMID: 17428026.
[46] P. Peumans, S. Uchida, and S. R. Forrest, “Efficient bulk heterojunction photovoltaic
cells using small-molecular-weight organic thin films,” Nature, vol. 425, pp. 158–162,
Sept. 2003.
204

[47] C. Groves, R. A. Marsh, and N. C. Greenham, “Monte carlo modeling of geminate
recombination in polymer-polymer photovoltaic devices,” Journal of Chemical Physics,
vol. 129, p. 114903, Sept. 2008.
[48] C. Y. Yang, J. G. Hu, and A. J. Heeger, “Molecular structure and dynamics at the
interfaces within bulk heterojunction materials for solar cells,” Journal of the American
Chemical Society, vol. 128, no. 36, pp. 12007–12013, 2006.
[49] S. Athanasopoulos, E. Hennebicq, D. Beljonne, and A. B. Walker, “Trap limited exciton transport in conjugated polymers,” Journal of Physical Chemistry C, vol. 112,
pp. 11532–11538, July 2008.
[50] D. S. Ginger and N. C. Greenham, “Charge separation in conjugatedpolymer/nanocrystal blends,” Synthetic Metals, vol. 101, pp. 425–428, May 1999.
[51] R. A. Marsh, C. Groves, and N. C. Greenham, “A microscopic model for the behavior
of nanostructured organic photovoltaic devices,” Journal of Applied Physics, vol. 101,
no. 8, p. 083509, 2007.
[52] R. A. Marcus, “Electron transfer reactions in chemistry - theory and experiment,”
Journal of Electroanalytical Chemistry, vol. 438, pp. 251–259, 1997.
[53] B. Ray, P. R. Nair, E. R. Garcia, and M. A. Alam, “Modeling and optimization of
polymer based bulk heterojunction(BH) solar cell,” IEEE Electron Devices Meeting,
vol. , no. , pp. 1–4, 2009.
[54] A. B. Walker, “Multiscale modeling of charge and energy transport in organic lightemitting diodes and photovoltaics,” Proceedings of the IEEE, vol. 97, pp. 1587–1596,
Sept. 2009.
[55] E. Hennebicq, G. Pourtois, G. D. Scholes, L. M. Herz, D. M. Russell, C. Silva, S. Setayesh, A. C. Grimsdale, K. Mllen, J.-L. Brdas, and D. Beljonne, “Exciton migration
in rigid-rod conjugated polymers: an improved frster model,” Journal of the American
Chemical Society, vol. 127, no. 13, pp. 4744–4762, 2005. PMID: 15796541.
[56] I. G. Scheblykin, A. Yartsev, T. Pullerits, V. Gulbinas, and V. Sundstrm, “Excited
state and charge photogeneration dynamics in conjugated polymers,” The Journal of
Physical Chemistry B, vol. 111, no. 23, pp. 6303–6321, 2007.
[57] Y. Kanai and J. C. Grossman, “Insights on interfacial charge transfer across
p3ht/fullerene photovoltaic heterojunction from ab initio calculations,” Nano Letters,
vol. 7, no. 7, pp. 1967–1972, 2007.
205

[58] Z. Li, X. Zhang, and G. Lu, “Dipole-assisted charge separation in organicinorganic
hybrid photovoltaic heterojunctions: Insight from first-principles simulations,” The
Journal of Physical Chemistry C, vol. 116, no. 18, pp. 9845–9851, 2012.
[59] T. W. Holcombe, J. E. Norton, J. Rivnay, C. H. Woo, L. Goris, C. Piliego, G. Griffini,
A. Sellinger, J.-L. Bredas, A. Salleo, and J. M. J. Frechet, “Steric control of the
donor/acceptor interface: Implications in organic photovoltaic charge generation,”
Journal of the American Chemical Society, vol. 133, no. 31, pp. 12106–12114, 2011.
[60] T. Liu and A. Troisi, “Absolute rate of charge separation and recombination in a
molecular model of the p3ht/pcbm interface,” The Journal of Physical Chemistry C,
vol. 115, no. 5, pp. 2406–2415, 2011.
[61] T. Liu, D. L. Cheung, and A. Troisi, “Structural variability and dynamics of the
p3ht/pcbm interface and its effects on the electronic structure and the charge-transfer
rates in solar cells,” Phys. Chem. Chem. Phys., vol. 13, pp. 21461–21470, 2011.
[62] F. J. Martinez-Veracoechea and F. A. Escobedo, “Bicontinuous phases in diblock
copolymer/homopolymer blends: Simulation and self-consistent field theory,” Macromolecules, vol. 42, no. 5, pp. 1775–1784, 2009.
[63] R. Hamilton, C. G. Shuttle, B. ORegan, T. C. Hammant, J. Nelson, and J. R. Durrant, “Recombination in annealed and nonannealed polythiophene/fullerene solar cells:
Transient photovoltage studies versus numerical modeling,” The Journal of Physical
Chemistry Letters, vol. 1, no. 9, pp. 1432–1436, 2010.
[64] 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,” Phys. Rev. B,
vol. 72, p. 085205, Aug 2005.
[65] G. A. Buxton and N. Clarke, “Predicting structure and property relations in polymeric
photovoltaic devices,” Physical Review B, vol. 74, p. 085207, Aug. 2006.
[66] I. Hwang, C. R. McNeill, and N. C. Greenham, “Drift-diffusion modeling of photocurrent transients in bulk heterojunction solar cells,” Journal of Applied Physics, vol. 106,
p. 094506, Nov. 2009.
[67] Y. M. Nam, J. Huh, and W. H. Jo, “Optimization of thickness and morphology of
active layer for high performance of bulk-heterojunction organic solar cells,” Solar
Energy Materials and Solar Cells, vol. 94, no. 6, pp. 1118 – 1124, 2010.

206

[68] N. Christ, S. Kettlitz, S. Zufle, S. Valouch, and U. Lemmer, “Trap states limited
nanosecond response of organic solar cells,” in Numerical Simulation of Optoelectronic
Devices (NUSOD), 2010 10th International Conference on, pp. 69–70, IEEE, 2010.
[69] C. R. McNeill, I. Hwang, and N. C. Greenham, “Photocurrent transients in all-polymer
solar cells: Trapping and detrapping effects,” Journal of Applied Physics, vol. 106,
p. 024507, July 2009.
[70] T. Kirchartz, B. E. Pieters, K. Taretto, and U. Rau, “Mobility dependent efficiencies
of organic bulk heterojunction solar cells: Surface recombination and charge transfer
state distribution,” Phys. Rev. B, vol. 80, p. 035334, Jul 2009.
[71] R. Li, Y. Peng, C. Ma, R. Wang, Y. Wang, H. Xie, T. Yang, J. Xie, S. Yan, and
J. Zhang, “Effect of mixture ratio on the performance of mdmo-ppv:pcbm bulk heterojunction solar cells: A numerical study,” Materials Science and Engineering: B,
vol. 172, no. 3, pp. 305 – 310, 2010.
[72] T. Kirchartz, W. Gong, S. Hawks, T. Agostinelli, R. MacKenzie, Y. Yang, and J. Nelson, “Sensitivity of the mott–schottky analysis in organic solar cells,” Journal of Physical Chemistry C, vol. 116, no. 14, p. 7672, 2012.
[73] C. G. Shuttle, R. Hamilton, B. C. ORegan, J. Nelson, and J. R. Durrant, “Chargedensity-based analysis of the currentvoltage response of polythiophene/fullerene photovoltaic devices,” Proceedings of the National Academy of Sciences, vol. 107, no. 38,
pp. 16448–16452, 2010.
[74] T. Kirchartz, J. Mattheis, and U. Rau, “Detailed balance theory of excitonic and bulk
heterojunction solar cells,” Phys. Rev. B, vol. 78, p. 235320, Dec 2008.
[75] L. Y. Meng, Y. Shang, Q. K. Li, Y. F. Li, X. W. Zhan, Z. G. Shuai, R. G. E. Kimber,
and A. B. Walker, “Dynamic monte carlo simulation for highly efficient polymer blend
photovoltaics,” Journal of Physical Chemistry B, vol. 114, no. 1, pp. 36–41, 2010.
[76] L. Y. Meng, D. Wang, Q. K. Li, Y. P. Yi, J. L. Bredas, and Z. G. Shuai, “An improved
dynamic monte carlo model coupled with poisson equation to simulate the performance of organic photovoltaic devices,” Journal of Chemical Physics, vol. 134, no. 12,
p. 124102, 2011.
[77] J. Anta, “Random walk numerical simulation for solar cell applications,” Energy &
Environmental Science, vol. 2, no. 4, pp. 387–392, 2009.

207

[78] J. C. Blakesley and D. Neher, “Relationship between energetic disorder and opencircuit voltage in bulk heterojunction organic solar cells,” Phys. Rev. B, vol. 84,
p. 075210, Aug 2011.
[79] C. Groves, R. G. E. Kimber, and A. B. Walker, “Simulation of loss mechanisms in
organic solar cells: A description of the mesoscopic monte carlo technique and an
evaluation of the first reaction method,” Journal of Chemical Physics, vol. 133, no. 14,
p. 144110, 2010.
[80] J. Nelson, J. J. Kwiatkowski, J. Kirkpatrick, and J. M. Frost, “Modeling charge transport in organic photovoltaic materials,” Accounts of Chemical Research, vol. 42, no. 11,
pp. 1768–1778, 2009.
[81] T. Offermans, S. Meskers, and R. Janssen, “Monte-carlo simulations of geminate
electron-hole pair dissociation in a molecular heterojunction: a two-step dissociation
mechanism,” Chemical Physics, vol. 308, pp. 125–133, JAN 10 2005.
[82] C. Groves and N. C. Greenham, “Bimolecular recombination in polymer electronic
devices,” Physical Review B, vol. 78, no. 15, p. 155205, 2008.
[83] B. P. Lyons, N. Clarke, and C. Groves, “The Quantitative Effect of Surface Wetting
Layers on the Performance of Organic Bulk Heterojunction Photovoltaic Devices,”
Journal of Physical Chemistry C, vol. 115, pp. 22572–22577, NOV 17 2011.
[84] C. Groves, J. C. Blakesley, and N. C. Greenham, “Effect of charge trapping on geminate recombination and polymer solar cell performance,” Nano Letters, vol. 10, no. 3,
pp. 1063–1069, 2010.
[85] B. Zacher and N. R. Armstrong, “Modeling the effects of molecular length scale electrode heterogeneity in organic solar cells,” The Journal of Physical Chemistry C,
vol. 115, no. 51, pp. 25496–25507, 2011.
[86] H. P. Yan, S. Swaraj, C. Wang, I. Hwang, N. C. Greenham, C. Groves, H. Ade, and
C. R. McNeill, “Influence of annealing and interfacial roughness on the performance of
bilayer donor/acceptor polymer photovoltaic devices,” Advanced Functional Materials,
vol. 20, no. 24, pp. 4329–4337, 2010.
[87] C. R. McNeill, S. Westenhoff, C. Groves, R. H. Friend, and N. C. Greenham, “Influence
of nanoscale phase separation on the charge generation dynamics and photovoltaic performance of conjugated polymer blends: Balancing charge generation and separation,”
Journal of Physical Chemistry C, vol. 111, no. 51, pp. 19153–19160, 2007.

208

[88] M. Casalegno, G. Raos, and R. Po, “Methodological assessment of kinetic monte carlo
simulations of organic photovoltaic devices: The treatment of electrostatic interactions,” Journal of Chemical Physics, vol. 132, no. 9, 2010.
[89] J. Williams and A. B. Walker, “Two-dimensional simulations of bulk heterojunction
solar cell characteristics,” Nanotechnology, vol. 19, no. 42, 2008.
[90] D. T. Gillespie, “Exact stochastic simulation of coupled chemical reactions,” The Journal of Physical Chemistry, vol. 81, no. 25, pp. 2340–2361, 1977.
[91] S. Westenhoff, W. J. D. Beenken, A. Yartsev, and N. C. Greenham, “Conformational
disorder of conjugated polymers,” Journal of Chemical Physics, vol. 125, p. 154903,
Oct. 2006.
[92] S. Athanasopoulos, J. Kirkpatrick, D. Martinez, J. M. Frost, C. M. Foden, A. B.
Walker, and J. Nelson, “Predictive study of charge transport in disordered semiconducting polymers,” Nano Letters, vol. 7, pp. 1785–1788, JUN 2007.
[93] S. Athanasopoulos, E. V. Emelianova, A. B. Walker, and D. Beljonne, “Exciton diffusion in energetically disordered organic materials,” Phys. Rev. B, vol. 80, p. 195209,
Nov. 2009.
[94] W. B. Russel, D. A. Saville, and W. R. Schowalter, Colloidal Dispersions. New York:
Cambridge University Press, 1989.
[95] N. Sutin, “Theory of electron transfer reactions: insights and hindsights,” Progress in
inorganic chemistry, vol. 30, pp. 441–498, 1983.
[96] J. R. Miller, L. T. Calcaterra, and G. L. Closs, “Intramolecular long-distance electron
transfer in radical anions. the effects of free energy and solvent on the reaction rates,”
Journal of the American Chemical Society, vol. 106, no. 10, pp. 3047–3049, 1984.
[97] J. A. Barker, C. M. Ramsdale, and N. C. Greenham, “Modeling the current-voltage
characteristics of bilayer polymer photovoltaic devices,” Physical Review B, vol. 67,
p. 075205, Feb. 2003.
[98] H. J. Snaith, N. C. Greenham, and R. H. Friend, “The origin of collected charge and
open-circuit voltage in blended polyfluorene photovoltaic devices,” Advanced Materials,
vol. 16, pp. 1640–+, Sept. 2004.
[99] M. Dresch, R. Isidoro, M. Linardi, J. Rey, F. Fonseca, and E. Santiago, “Influence of
sol-gel media on the properties of nafion-sio¡ sub¿ 2¡/sub¿ hybrid electrolytes for high
209

performance proton exchange membrane fuel cells operating at high temperature and
low humidity,” Electrochimica Acta, 2012.
[100] E. J. Roche, M. Pineri, R. Duplessix, and A. M. Levelut, “Small-angle scattering
studies of nafion membranes,” Journal of Polymer Science: Polymer Physics Edition,
vol. 19, no. 1, pp. 1–11, 1981.
[101] G. Gebel and J. Lambard, “Small-angle scattering study of water-swollen perfluorinated ionomer membranes,” Macromolecules, vol. 30, no. 25, pp. 7914–7920, 1997.
[102] D. Shah, P. Maiti, D. Jiang, C. Batt, and E. Giannelis, “Effect of nanoparticle mobility
on toughness of polymer nanocomposites,” Advanced Materials, vol. 17, no. 5, pp. 525–
528, 2005.
[103] A. Huq, R. Welberry, and E. Bozin, “Advances in structural studies of materials using
scattering probes,” MRS bulletin, vol. 35, no. 07, pp. 520–530, 2010.
[104] D. Rangappa, K. Murukanahally, T. Tomai, A. Unemoto, and I. Honma, “Ultrathin
nanosheets of li2msio4 (m= fe, mn) as high-capacity li-ion battery electrode,” Nano
letters, vol. 12, no. 3, pp. 1146–1151, 2012.
[105] D. Spagnoli and J. Gale, “Atomistic theory and simulation of the morphology and
structure of ionic nanoparticles,” Nanoscale, vol. 4, no. 4, pp. 1051–1067, 2012.
[106] E. Lorenz, C. Keil, and D. Schlettwein, “Structure and morphology in thin films of
perfluorinated copper phthalocyanine grown on alkali halide surfaces (001),” Journal
of Porphyrins and Phthalocyanines, vol. 16, no. 07n08, pp. 977–984, 2012.
[107] A. McGlashon, W. Zhang, D. Smilgies, M. Shkunov, K. Genevicius, K. Whitehead,
A. Amassian, G. Malliaras, D. Bradley, M. Heeney, et al., “Spectroscopic and morphological investigation of conjugated photopolymerisable quinquethiophene liquid crystals,” Current Applied Physics, vol. 12, pp. e59–e66, 2012.
[108] H. Wang, H. Jeong, M. Imura, L. Wang, L. Radhakrishnan, N. Fujita, T. Castle,
O. Terasaki, and Y. Yamauchi, “Shape-and size-controlled synthesis in hard templates:
Sophisticated chemical reduction for mesoporous monocrystalline platinum nanoparticles,” Journal of the American Chemical Society, vol. 133, no. 37, pp. 14526–14529,
2011.
[109] D. Williams and C. Barry, Transmission Electron Microscopy: A Textbook for Materials
Science. Springer, 2009.
210

[110] D. S. Sivia, Elementary Scattering Theory for X-ray and Neutron Users. Oxford University Press, 2011.
[111] A. Guinier, G. Fournet, C. Walker, and K. Yudowitch, Small-angle scattering of X-rays,
vol. 14. Wiley New York, 1955.
[112] J. O’Donnell, X. Zuo, A. Goshe, L. Sarkisov, R. Snurr, J. Hupp, and D. Tiede,
“Solution-phase structural characterization of supramolecular assemblies by molecular
diffraction,” Journal of the American Chemical Society, vol. 129, no. 6, pp. 1578–1585,
2007.
[113] D. Schneidman-Duhovny, M. Hammel, and A. Sali, “Foxs: a web server for rapid
computation and fitting of saxs profiles,” Nucleic acids research, vol. 38, no. suppl 2,
pp. W540–W544, 2010.
[114] F. Poitevin, H. Orland, S. Doniach, P. Koehl, and M. Delarue, “Aquasaxs: a web server
for computation and fitting of saxs profiles with non-uniformally hydrated atomic
models,” Nucleic acids research, vol. 39, no. suppl 2, pp. W184–W189, 2011.
[115] D. Tiede, K. Mardis, and X. Zuo, “X-ray scattering combined with coordinate-based
analyses for applications in natural and artificial photosynthesis,” Photosynthesis research, vol. 102, no. 2, pp. 267–279, 2009.
[116] H.-C. Liao, C.-S. Tsao, T.-H. Lin, C.-M. Chuang, C.-Y. Chen, U.-S. Jeng, C.-H. Su,
Y.-F. Chen, and W.-F. Su, “Quantitative nanoorganized structural evolution for a high
efficiency bulk heterojunction polymer solar cell,” Journal of the American Chemical
Society, vol. 133, no. 33, pp. 13064–13073, 2011.
[117] N. A. Gumerov, K. Berlin, D. Fushman, and R. Duraiswami, “A hierarchical algorithm
for fast debye summation with applications to small angle scattering,” Journal of
Computational Chemistry, vol. 33, no. 25, pp. 1981–1996, 2012.
[118] S. R. Kline, “Reduction and analysis of sans and usans data using igor pro,” Journal
of Applied Crystallography, vol. 39, no. 6, pp. 895–900, 2006.
[119] J. P´rez and Y. Nishino, “Advances in x-ray scattering: from solution saxs to achievee
ments with coherent beams,” Current Opinion in Structural Biology, 2012.
[120] R. Jones, S. Kumar, D. Ho, R. Briber, and T. Russell, “Chain conformation in ultrathin
polymer films using small-angle neutron scattering,” Macromolecules, vol. 34, no. 3,
pp. 559–567, 2001.
211

[121] H. C. B. J. S. Higgins, Polymers and Neutron Scattering. Oxford Science Publications,
1996.
[122] B. H. Zimm, “Apparatus and methods for measurement and interpretation of the
angular variation of light scattering; preliminary results on polystyrene solutions,”
The Journal of Chemical Physics, vol. 16, no. 12, pp. 1099–1116, 1948.
[123] H. Hoppe and N. S. Sariciftci, “Morphology of polymer/fullerene bulk heterojunction
solar cells,” J. Mater. Chem., vol. 16, no. 1, pp. 45–61, 2006.
[124] H. Hoppe, M. Niggemann, C. Winder, J. Kraut, R. Hiesgen, A. Hinsch, D. Meissner,
and N. Sariciftci, “Nanoscale morphology of conjugated polymer/fullerene-based bulkheterojunction solar cells,” Advanced Functional Materials, vol. 14, no. 10, pp. 1005–
1011, 2004.
[125] J. Jo, S.-I. Na, S.-S. Kim, T.-W. Lee, Y. Chung, S.-J. Kang, D. Vak, and D.-Y.
Kim, “Three-dimensional bulk heterojunction morphology for achieving high internal
quantum efficiency in polymer solar cells,” Advanced Functional Materials, vol. 19,
no. 15, pp. 2398–2406, 2009.
[126] P. Vanlaeke, A. Swinnen, I. Haeldermans, G. Vanhoyland, T. Aernouts, D. Cheyns,
C. Deibel, J. DHaen, P. Heremans, J. Poortmans, and J. Manca, “P3ht/pcbm bulk
heterojunction solar cells: Relation between morphology and electro-optical characteristics,” Solar Energy Materials and Solar Cells, vol. 90, no. 14, pp. 2150 – 2158,
2006.
[127] M. S. White, D. C. Olson, S. E. Shaheen, N. Kopidakis, and D. S. Ginley, “Inverted
bulk-heterojunction organic photovoltaic device using a solution-derived zno underlayer,” Applied Physics Letters, vol. 89, no. 14, p. 143517, 2006.
[128] Z. Xu, L.-M. Chen, G. Yang, C.-H. Huang, J. Hou, Y. Wu, G. Li, C.-S. Hsu, and
Y. Yang, “Vertical phase separation in poly(3-hexylthiophene): Fullerene derivative
blends and its advantage for inverted structure solar cells,” Advanced Functional Materials, vol. 19, no. 8, pp. 1227–1234, 2009.
[129] R. G. E. Kimber, A. B. Walker, G. E. Schroder-Turk, and D. J. Cleaver, “Bicontinuous
minimal surface nanostructures for polymer blend solar cells,” Physical Chemistry
Chemical Physics, vol. 12, no. 4, pp. 844–851, 2010.
[130] C. F. Majkrzak and N. F. Berk, “Exact determination of the phase in neutron reflectometry,” Phys. Rev. B, vol. 52, pp. 10827–10830, Oct 1995.

212

[131] E. Teller, N. Metropolis, and A. Rosenbluth, “Equation of state calculations by fast
computing machines,” J. Chem. Phys, vol. 21, no. 13, pp. 1087–1092, 1953.
[132] R. Glauber, “Time-dependent statistics of the ising model,” Journal of mathematical
physics, vol. 4, p. 294, 1963.
[133] K. Kawasaki, “Kinetics of ising models,” Phase transitions and critical phenomena,
vol. 2, pp. 443–501, 1972.
[134] A. Nelson, “Co-refinement of multiple-contrast neutron/X-ray reflectivity data using
M OTOFIT,” Journal of Applied Crystallography, vol. 39, pp. 273–276, Apr 2006.
[135] A. Salleo, R. J. Kline, D. M. DeLongchamp, and M. L. Chabinyc, “Microstructural
characterization and charge transport in thin films of conjugated polymers,” Advanced
Materials, vol. 22, no. 34, pp. 3812–3838, 2010.
[136] D. A. Chen, A. Nakahara, D. G. Wei, D. Nordlund, and T. P. Russell, “P3ht/pcbm
bulk heterojunction organic photovoltaics: Correlating efficiency and morphology,”
Nano Letters, vol. 11, no. 2, pp. 561–567, 2011.
[137] 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, vol. 2,
pp. 3135–3145, DEC 15 2011.
[138] C. R. McNeill, “Morphology of all-polymer solar cells,” Energy and Environmental
Science, vol. 5, pp. 5653–5667, 2012.
[139] D. DeLongchamp, R. Kline, and A. Herzing, “Nanoscale structure measurements for
polymer-fullerene photovoltaics,” Energy Environ. Sci., vol. 5, no. 3, pp. 5980–5993,
2012.
[140] W. Oosterbaan, J. Bols´e, A. Gadisa, V. Vrindts, S. Bertho, J. D’Haen, T. Cleij,
e
L. Lutsen, C. McNeill, L. Thomsen, et al., “Alkyl-chain-length-independent hole mobility via morphological control with poly (3-alkylthiophene) nanofibers,” Advanced
Functional Materials, vol. 20, no. 5, pp. 792–802, 2010.
[141] I. Hwang and N. C. Greenham, “Modeling photocurrent transients in organic solar
cells,” Nanotechnology, vol. 19, p. 424012, Oct. 2008.
[142] L. D. Landau and E. Lifshitz, Electrodynamics of Continuous Media - 2nd ed. Pergamon Press, 1984.
213

[143] C. Farrow, M. Shaw, H. Kim, P. Juh´s, and S. Billinge, “Nyquist-shannon sampling
a
theorem applied to refinements of the atomic pair distribution function,” Physical
Review B, vol. 84, no. 13, p. 134105, 2011.

214