A 2...:
5..» runs”
a. 2L .

ﬁx

 

mm at

£35

3
O

p. ‘
. 4.0.. :I
31-. L in“.

 

‘ ,_..m..uv . . . ‘ . . ‘ , ... c , ‘ ‘ basin. 472.»... ...

i-

 

 

THESIS

Z Illlllllllllllllllllllllllnlllllllllllllllllllllllllll

193 01688 7519

This is to certify that the

dissertation entitled
1‘ g‘ley 0f ConcludjVT-b/ ‘m a Comp/9X 960144557.

presented by

36mg“ HyKAVL

has been accepted towards fulﬁllment
of the requirements for

Pk ~ D. degree in WMW/

141%

Major professor?

 

’_—

Date 33/), 5/ 0

MSU is an Afﬁrmative Action/Equal Opportunity Institution 0-12771

 

 

 

LIBRARY

Michigan State

University

 

 

PLACE lN RETURN BOX
to remove this checkout from your record.
TO AVOID FINES return on or before date due.

 

DATE DUE

MTE DUE

DATE DUE

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

use wmmn

A STUDY OF CONDUCTIVITY IN A
COMPLEX GEOMETRY

By

Sangil Hyun

A DISSERTATION

Submitted to
Michigan State University
In partial fulﬁllment of the requirements
For the degree of

DOCTOR OF PHILOSOPHY
Department of Physics and Astronomy

1998

ABSTRACT
A STUDY OF CONDUCTIVITY IN A COMPLEX GEOMETRY
By

Sangil Hyun

I present a study of the conductivity1 of a material with a complex geometry.
There are three main parts in this dissertation, the determination of the isotropic
conductivity, anisotropic conductivity, and effective conductivity. In the ﬁrst part, a
method using a multi-probe measurement together with a numerical simulation is
introduced for the determination of the isotropic conductivity in a complex geometry.
This method has been demonstrated on micron-sized diamond crystallites and diamond
homoepitaxial ﬁlms. The 3d-computer image for the ﬁnite element method (FEM) has
been reconstructed from the scanning electron microscope (SEM) photos. For a sample
with four probes, a 6 x 6 geometrical factor matrix is generated from the FEM analysis. In
addition, the experimental measurement on the sample provides another 6 x 6 resistance
matrix. The sample conductivity and contact resistances can be determined by the least-
square fitting for the geometrical factors and the resistances. In the second part, this
method is generalized for the anisotropic conductivity that is represented by a 3 x 3
tensor. It has been validated on computer models and real material (bismuth). This
modiﬁed numerical scheme using the FEM analysis and the iterative linearization
technique has identiﬁed the anisotropic conductivity tensor within acceptable errors in

most cases. Some conditions for the application of this technique have been suggested

from the analysis of the results. In the last part, I describe the effective conductivity of a
composite material. Much research has been done to investigate the electric, thermal, and
elastic properties of composite materials. However, only a few problems with a simple
geometry have been solved analytically. For instance, the dielectric properties of a
composite that includes two perfect circular conductors have been determined analytically
by a multipole expansion method. However, for complex composites, the analytic
solutions cannot be obtained in most cases. Therefore, a numerical technique using the
ﬁnite element method has been introduced for problems where the analytic approach
fails. The FEM analysis is known as the most effective scheme to solve a boundary value
problem in a complicated geometry. The simulation using this method is a primary step to
understanding the effective medium theory for a complex composite material. To
demonstrate this technique in the study of composite materials, numerical simulations
have been performed and the results have been compared with the analytical solutions for

a simple composite material.

To all whom I am indebted for this,

especially North Koreans in the horrible famine.

ACKNOWLEDGMENTS

I would like to thank Professor M. F. Thorpe for his considerate and patient
guidance over the years. His insights into physics and passions on the study have made a
deep impression on me. Also, I like to appreciate Professors S. D. Mahanti, N. O. Birge,
D. Stump, and S. Hawley for being my guidance committee members. The collaboration
with Professor B. Golding, Professor A. R. Day, and Dr. M. D. Jaeger has been a great
pleasure and experience to me. I want to acknowledge Ms. J. King for her help.

I appreciate many people who have shared their lives in this small building,
especially, H. Seong, K. Chun, J. S. Moon, B. Djordjevic, and J. U. Kim. They have made
my life here always delightful. I appreciate also L. Malete for his helping on my
manuscript. I know there have been invaluable and uncountable help from many people
who have not been named here. They should be acknowledged, too.

My study would not be accomplished without the prayer, support, and patience of
my dear parents, M. G. Hyun and G. J. Kim, and whole family including a newly born
nephew, SeongWon, whom I have not yet meet. My deep thanks should go to them.

I want to represent a special gratitude to my beloved Myeoung H. Oh who has
shared the valuable and memorable times with me. I thank you more than I can express.

Finally, my thanks would go to Him who have created the world and me. In my

whole life, He has lead me to the truth with His boundless love.

TABLE OF CONTENTS

 

 

LIST OF TABLES ................................................................................ viii
LIST OF FIGURES ................................................................................ X
CHAPTER 1 INTRODUCTION 1

CHAPTER 2 CONDUCTIVITY MEASUREMENT OF MICRON-SIZED DIAMOND

CRYSTALLITES 4
2.1 INTRODUCTION .................................................................................................................. 4
2.1.1 motivation ............................................................................................................... 4
2.1.2 van der Pauw formula for 2d ﬁlm .......................................................................... 8
2.1.3 inverse problem .................................................................................................... 1 1
2.1.4 Finite Element Method (FEM) ............................................................................. 15
2.2 THEORY ........................................................................................................................... 23
2.2.1 conductivity and geometrical factors ................................................................... 23
2.2.2 resistance table and its properties ....................................................................... 25
2.2.3 contact resistances ............................................................................................... 31
2.3 COMPUTER SIMULATION .................................................................................................. 34
2.3.1 reconstruction of 3d computer image ................................................................... 34
2.3.2 modeling and analysis using FEM ....................................................................... 39
2.3.3 least square ﬁtting ................................................................................................ 47
2.4 EXPERIMENT .................................................................................................................... 50
2.4.1 sample preparation ............................................................................................... 50
2.4.2 four-probe measurement ...................................................................................... 53
2.5 ANALYSIS ......................................................................................................................... 55
2.5.1 diamond ﬁlm ......................................................................................................... 55
2.5.2 sample I (Ohmic contacts) ................................................................................... 60
2.5.3 sample 2 (non-Ohmic contacts) ............................................................................ 62
2.6 CONCLUSION .................................................................................................................... 68

CHAPTER 3 ANISOTROPIC CONDUCTIVITY MEASUREMENT FOR COMPLEX

 

GEOMETRY 70
3.] INTRODUCTION ............................................................................................................... 70
3.2 THEORY .......................................................................................................................... 72
3.2.1 anisotropic conductivity tensor ............................................................................ 72
3.2.2 iterative linearization technique for nonlinear inverse problem ......................... 74
3.2.3 orthogonal transformation ................................................................................... 80
3.3 COMPUTER SIMULATION ................................................................................................ 81
3.3.1 sample geometry ................................................................................................... 81
3.3.2 simulated experimental data using FEM ............................................................. 83
3.3.3 iterative linearization ﬁtting using FEM ............................................................. 84
3.4 EXPERIMENT .................................................................................................................. 85
3.4.1 sample preparation (Bi cylinder) ......................................................................... 85
3.4.2 six-probe measurement ......................................................................................... 88
3.5 ANALYSIS ....................................................................................................................... 89
3.5.1 computer generated sample (2d) .......................................................................... 89

vi

3.5.2 computer generated sample (3d) ................................................................................ 92

3.5.3 real sample (Bi cylinder) ............................................................................................ 96
3.5.4 stability test ............................................................................................................... 101
3.6 CONCLUSION ....................................................................................................................... 106

CHAPTER 4 EFFECTIVE CONDUCTIVITY OF A MEDIUM WITH CIRCULAR

 

 

INCLUSIONS 108

4.1 INTRODUCIION ............................................................................................................. 108

4.2 A BRIEF HISTORY OF THE STUDY ON A COMPOSITE MATERIAL .................................... 109

4.3 THEORY ........................................................................................................................ 1 13
4.3.1 4-point resistance (R4) ........................................................................................ 114

4.3.2 Z-point resistance (R2) ........................................................................................ 118

4.3.3 The ratio of R4 and R2 ..................................................................................... 121

4.4 COMPUTER SIMULATION .............................................................................................. 122
4.4.] FEM analysis for 2d composite .......................................................................... 122

4.4.2 FEM analysis for 3d composite .......................................................................... 124

4.5 ANALYSIS ..................................................................................................................... 126
4.5.1 2d composite ....................................................................................................... 126

4.5.2 3d composite ....................................................................................................... 128

4.5.3 The ratio of 4-point and 2-point resistances ...................................................... 130

4.6 CONCLUSION ................................................................................................................ 130
APPENDICES 133
A.I RECIPROCITY THEOREM .............................................................................................. I34
A.2 COUNTING OF INDEPENDENT 4-TERMINAL TERMS ....................................................... 138

A.3 RESISTANCES OF TWO PERFECT HYPERSPHERICAL CONDUCTORS ............................. I40
REFERENCES 145

 

vii

Table 2.1

Table 2.2

Table 2.3

Table 2.4

Table 2.5

Table 2.6

Table 2.7

Table 2.8

Table 2.9

Table 2.10

Table 2.11

Table 2.12

Table 2.13

Table 2.14

Table 2.15

Table 2.16

Table 2.17

Table 2.18

Table 2.19

Table 2.20

Table 2.21

LIST OF TABLES

 

Selected properties of diamond structure semiconductors 4

A resistance table from all possible combinations of voltage and current
measurements 27

 

How to generate the triangular matrix from six diagonal terms so

 

 

 

 

 

that the triangular matrix produces a whole table 30
The number of independent elements in the resistance table 30
Resistance table with two contact resistances for current flows 32
The characteristics of the meshings for each sample 42
The Ohmic resistance table obtained from Table 2.5 (symmetric) 48
The non-Ohmic resistance table obtained from Table 2.5 (non-symmetric) ----- 48
The characteristics of CVD diamond samples 52

 

Resistances measured for homoepitaxial ﬁlm by 4—probe experiment (k9) ----- 56

 

Geometrical factor table obtained by FEM (mm") 56

The resistivity and contact resistances of diamond ﬁlm (Film 1) given by

 

 

 

van der Pauw formula (VDP) and our method. 56
Resistances measured by 4—probe experiment (1(9) 58
Another resistance table for the opposite current ﬂows of Table 13 (k9)-------- 58
Geometrical factors calculated from FEM (mm") 58

The resistivity and contact resistances of Film 2 given by van der Pauw
formula (VDP) and our method. 58

 

Stability test of 4—terminal resistances on the position of one contact -- -------- 59
Resistances obtained from the experimental measurement for Sample 1 (kﬂ)-—— 61
Geometrical factors from the ﬁnite element method for Sample 1 (um")-------- 61

Resistivity and contact resistances for Sample 1, assuming that all the contact
resistances were Ohmic. 61

 

Resistances obtained from the experimental measurement for Sample 2 (kQ)--- 64

viii

Table 2.22 Geometrical factors from the ﬁnite element method for Sample 2 (um’l )------ 64

Table 2.23 Resistivity and contact resistances for Sample 2 64

 

Table 2.24 Fit result of conductivity model with two activation energies for data in Figure
2.27 66

 

Table 3.1 Electrical resistivity of metalic crystals. Values of principal resistivities at
T=20° C. (in unit 1045 (2 cm) 71

 

Table 3.2 The comparison of the ﬁtted results of the conductivity tensor elements with the
prescribed values for the computer model in 2d 91

 

Table 3.3 The comparison of the ﬁtted results of the conductivity tensor elements with the

 

 

 

prescribed values for the computer model in 3d 93
Table 3.4 Principal values of the resistivity of Bi. (104’ 9 cm) 99
Table 3.5 Different meshings for the stability test of a cylindrical sample 102
Table 3.6 Different meshings for the stability test of Cylinder2 106

 

Table A.l The number of independent 4-terminal terms in a resistance matrix for n-probe
measurement 1 39

 

ix

Figure 2.1

Figure 2.2

Figure 2.3

Figure 2.4

Figure 2.5

Figure 2.6

Figure 2.7

Figure 2.8

Figure 2.9

Figure 2.10
Figure 2.11

HngH

Figure 2.13

Figure 2.14

Figure 2.15

Figure 2.16

LIST OF FIGURES

A polycrystalline diamond ﬁlm (a) and a bicrystal (b) that was formed by two
single crystallites 5

 

A typical shape of CVD diamonds used in this measurement. It is a micron-sized
sample grown on a Si substrate. 6

 

The conductivity measurement of a simple shaped sample is usually
straightforward. This cylindrical geometry forms constant electric current pattern
inside the sample 7

 

The conductivity measurement of 2D ﬁlm in arbitrary geometry is done by four-
probe measurement. 9

 

A clover-shaped ﬁlm (a) and a bridge-shaped ﬁlm (b) are used for four-probe
measurement. 10

 

A schematic of the inverse and direct problem. Those two problems can be
characterized by what is known and what is determined in the problem.----- 12

Heat conduction problem on a simple bar 13

 

A boundary value problem is solved in this irregular geometry. Boundary
conditions on each surface are needed. 16

 

Some typical element shapes are presented for FEM. Each element has nodes on
the vertices or some extra nodes on the edges to enhance the accuracy. 17

 

Highly distorted meshes that are not usually acceptable in mesh generating ---- 18
A sample of meshed geometry with square-shaped elements in 2d ----------- 19

By the FEM analysis, a voltage distribution on a sample is obtained. The
interpolation technique is used to obtain the distribution inside the elements.--- 22

For 2-terminal measurement, the voltage is calculated by a line integral along the
path I and the current is obtained by a surface integral on the cross section S.-- 23

Showing a conductor that has four probes that are attached on the surface mm 24

Micron-sized diamond crystal grown by a chemical vapor deposition technique
before Au leads were attached. Five titanium contacts were formed on the

surface and four of them (indexed) were used for the measurement. --- 26

 

Contact resistances have two different values (I; , r, ') for non-Ohmic contacts. A

resistance measurement consists of a sample resistance Rm“, and contact
resistances added in series. 32

 

Figure 2.17

Figure 2.18

Figure 2.19

Figure 2.20

Figure 2.21

Figure 2.22

Figure 2.23

Figure 2.24

Figure 2.25

Figure 2.26

Figure 2.27

Figure 3.1

Each SEM photo is given in an orientation for viewing angle (9, 4)) as shown in
this ﬁgure. 34

 

SEM photos of a sample (Sample 2) in various orientations. (a) 9, ¢ = 0° , 0°
(b) 6, ¢ = 50° , 0° (c) 9, ¢ = 50° , 90° ((1) 6, ¢ = 50° , 270° 35

 

Computer images reconstructed from SEM photos by image processing. Five small
pyramids on the surface denote ﬁve contacts. Compare these two images with the
SEM photos in the Figure 2.18. 38

 

Triangular planes generated from three keypoints were combined into a large plane
as shown in (a). A volume could be constructed by connecting these large planes
into a volume unit. The volume elements for the diamond crystallites were
presented in (b) and (c). 40

 

The elements used for the meshing of 2d ﬁlm (square) and 3d crystallites
(tetrahedra). Sometimes elements that have only four nodes on the comers were
used. 42

 

Meshed structures of the diamond ﬁlm(a) and crystallites (Sample 1(b), Sample 2
(c)). Figure 2.20(b) and (c). Each model has square-shaped contacts on the surface.
43

 

Equipotential lines on the surface of the samples given by ABAQUS. (a) is for
homoepitaxial ﬁlm. (b), (c) are for Sample 1. (d),(e) are for Sample 2. The ﬁgures
of 3d samples were obtained in different orientations. 46

 

A normal picture of the 2d sample (diamond homoepitaxial ﬁlm) is shown in (a).
Two schematics ((b), (c)) of the sample denote the positions of the contacts. Film 1
(b) has four contacts and ﬁlm 2 (c) has six contacts. Gray-colored squares in the
ﬁgures represent the contacts actually used for the measurement. 51

 

Two micron-sized CV D diamond crystallites, Sample 1 (a) and Sample 2 (b).
Sample 1 has Au leads attached on the Ti contacts and Sample 2 is a shape before
the lead attachment 51

 

Temperature dependence of a resistance of diamond crystallites (Sample 1 (solid
circles) and Sample 2 (open circles)). To investigate the temperature dependence

of the resistivity of these samples, a four-terminal resistance (R5) were measured

 

over a wide range of temperatures. 54

Temperature dependence of the resistivities of the diamond samples. Two
microcrystallites (MC-A, MC-B) correspond to Sample 1 and Sample 2. PF 1 and
PF2 are polycrystalline ﬁlms. HF 1 and HF2 are homoepitaxial ﬁlms. The
resistivity of these samples was obtained by rescaling the resistance measurement
in Figure 2.26 using the ﬁtted resistivity. 67

 

A diamond crystallite has a local defect (twin boundary) denoted by an arrow. This
defect was neglected in the isotropic conductivity determination. 70

 

xi

Figure 3.2

Figure 3.3

Figure 3.4

Figure 3.5

Figure 3.6

Figure 3.7

Figure 3.8

Figure 3.9

Figure 3.10

Figure 3.11

Figure 3.12

Figure 3.13

Iterative linearization technique of 12 minimization in 1d (Newton method). This

technique is useful when the function form of the 12 is not known. Each step, the

function form is approximated to a linearized form by Taylor expansion. After
many iterations, it reaches the minimum point. 76

 

Many probes are attached on the surface of samples in 2d (a) and 3d (b). More
contacts are needed for anisotropic conductivity tensor compared to isotropic
conductivity tensor. 82

 

Schematics of a real sample (Bi cylinder) with six circular Au electrodes. (a) It is a
relatively thin cylinder. Three contact pads were formed on each circular plane of
the sample. (b) Three pads on a plane were formed to have an angle from three
pads on the other plane. 86

 

The normal photos of the Bi sample for the real measurement.2 It has six contacts
that are formed on both sides of the circular planes. (a) and (b) show the formation
of contacts on these two sides. The ruler shows the size of the sample. 87

 

A computer generated sample for the determination of anisotropic conductivity.
This 2d computer model has six circular probes. The contour map (b) shows the
equipotential lines. The meshing was done by ANSYS and the analysis was done
by ABAQUS. 89

 

These ﬁgures show the ﬁtting procedures of the tensor in 2d. The ﬁttings were
performed using only 4-point resistances (a) or whole resistances in the table (b).
Three data points are equivalent to three conductivity tensor terms.--—------ 90

Two principal directions were drawn on the surface of 2d computer model. Three
directions (one real and two ﬁtted) represent a remarkable coincidence.—- 91

 

The meshing (a) and equipotential lines (b) of a 3d computer generated sample.
Six brick-shaped probes were formed on six facets of the sample. 94

 

These ﬁgures show the ﬁtting procedures of the tensor in 3d. The ﬁttings were
performed using only 4-point resistances (a) or whole resistances in the table (b).
Six parameters in a tensor were obtained each step. The iterations in 3d needed a
few more steps to the optimum than in 2d. 95

 

Sensitivity test of the conductivities ﬁtted by 4-point terms and whole table has
been done to check the relation between the result and the Gaussian noise.---- 97

Equipotential lines of the real sample (Bi) obtained by FEM analysis. -------- -- 98

Fitting procedures of a resistivity tensor of the real sample. Two meshings were
used to test our method on this sample. Only three principal components of the
tensor were shown in these pictures. Two planar resistivity components (solid
circle and square) are observed stable around the real value, but the perpendicular
component (triangle) shows a substantial variation, especially in (b). —-—----- 100

xii

Figure 3.14

Figure 3.15

Figure 3.16

Figure 3.17

Figure 3.18

Figure 4.1

Figure 4.2

Figure 4.3

Figure 4.4

Figure 4.5

Figure 4.6

Resistivity ﬁttings were done for the Bi sample using three different meshings.
BilO—l is a ﬁne meshing and Bi4-2 is a rough meshing. It turned out that the
perpendicular resistivity component was so sensitive on the meshing while the
other components were stable. 103

 

Resistivity ﬁttings were performed on some modiﬁed structures. The
modiﬁcations were done by rotating the contact pads. The dependence of the result
on the angles was shown in the ﬁgure. Huge variations were observed for all

components. Interestingly, two planar resistivity components ( p,r , py) crossed

over somewhere between 2° and 5° . 103

 

Equipotential lines of thicker cylindrical sample (Cylinder 2). This model has more
layers of elements along the perpendicular direction. That may enhance the ﬁtting
result for that direction. 104

 

This graph was obtained from another computer test for the same geometry of the
Bi cylinder (Cylinder l). The simulated resistances were generated by the meshing
Bi10-l. The ﬁtted results from different meshings were presented. The
perpendicular component (triangle) shows still a huge variation for this computer
model. 105

 

Another computer test for a thick sample (Cylinder 2) as shown in Figure 3.16.
The simulated resistances were generated by the meshing BulklO—l. The ﬁtting
from other meshings could reproduce the all resistivity components nicely
including the perpendicular component. 105

 

A composite material consists of two different dielectric constant regions (t:0 , 8,.)
in dilute limit. 110

 

Elastic energy distribution in a composite material with many circular inclusions
under hydrostatic pressure. The contour map was given by the elasticity analysis
using ABAQUS. Elastic energy turned out to be distributed mostly near the necks
between two circles. 112

 

4-point resistance measurement of two circular inclusions. The voltage between
two inclusions is calculated in terms of the geometrical parameters when a
constant electric ﬁeld is applied. 115

 

Normalized current 10 is deﬁned by the total current flowing through the cross
section of the inclusions. 116

 

One dimensional composite material with two perfect conductor inclusions.--- 117

2-point resistance measurement of two circular inclusions. A resistance is obtained
when the current ﬂows between these two inclusions. 119

 

xiii

Figure 4.7

Figure 4.8

Figure 4.9

Figure 4.10

Figure 4.11

Figure 4.12

Figure 4.13

Figure 4.14

Figure A.1

Figure A2

The geometry of 2d composite material with circular inclusions. Huge plates were
attached to the host to form a constant electric ﬁeld in 4-point resistance
measurement. Those plates were removed for 2-point resistance measurement - 123

Contour graphs of voltage distribution near the neck in 2d. (a) is given from 4-
point measurement and (b) is given from 2-point measurement. 123

 

Three—dimensional composite material with two spherical inclusions. Similar large
bricks were attached to the host for 4-point measurement 125

 

Contour graph of voltage distribution near the neck in 3d. (a) is given from 4-point
measurement and (b) is given from 2-point measurement 125

 

The meshed structure of the spherical inclusions. For a rough meshing, as these
spheres approach together, the shape of the neck is getting less realistic because of
the sharp ends of the spheres 126

 

Numerical solutions of the 4-point and 2-point resistances were obtained for 2d
composite using FEM analysis. They agree well in all regions with the analytical
solutions. (b) shows the detail in the close neck limit 127

 

Numerical solutions of the 4-point and 2-point resistances were obtained for 3d
composite using FEM analysis. They agree nicely in most regions with the
analytical solutions, but show some disagreements in an asymptotic regime (close
limit). (b) shows the detail in that limit 129

 

The ratio between 4-point resistance and 2-point resistance in 2d and 3d. Each
ratio was observed to approach a constant in the close neck limit ( w/ a << 1 ). The
limiting value is given by Itin 2d and n2 / 3 in 3d 131

 

A sample with four probes and its equivalent circuit. The electric circuit is
composed of many resistors that are fully connected 135

 

Two circular perfect conductors in a host material in the close neck limit ---—--- 141

xiv

Chapter 1 INTRODUCTION

The determination of the physical properties of a material has been one of the
main topics in experimental and theoretical research over the years. It has been a primary
step to characterize a speciﬁc material in most experiments. Moreover, the study has been
considered as an important approach to investigate the structural dependence of material
properties for some structured materials such as glasses, amorphous solids, and
composites. The effective material properties have been considered for a composite that
consists of more than one material. Some studies have been conducted to ﬁnd an
optimum structure associated with these physical properties (electric, mechanical, and
thermal). The determination of these properties is usually straightforward in a simple
geometry. However, it is not possible for complex structures such as an arbitrary-shaped
diamond crystallite. The trouble in the determination usually comes from the complicated
geometry of the material. This is one of main motivations of this study. The ﬁnite element
method (FEM) has been adapted to deal with this trouble numerically. This technique is a
powerful tool to solve a boundary value problem especially in a complex geometry.
Sometimes, another computer simulation technique could be used to enhance the
performance of this method. This numerical technique in the determination of the
material properties provides a way to overcome a major limitation that results from
complex geometry.

In this dissertation, I describe the numerical determination of material properties.
Throughout my work, only the electric properties (conductivity, resistance) have been

considered. However, the methods suggested in this study can be applied to the

determination of other material properties (thermal, elastic, etc.), too. In Chapter 2, I
present the numerical method for the determination of the conductivity of micron-sized
diamond crystallites. For this diamond sample, the conductivity is isotropic due to its
symmetric lattice structure. The method is demonstrated on a diamond ﬁlm for the
comparison with the well-known van der Pauw method. Two main procedures, the multi-
probe measurement and the FEM analysis are required to identify the sample
conductivity. This work could be considered as the ﬁrst measurement of a diamond
crystallite in micron scale. In addition to conductivity, contact resistances are determined
also quantitatively. The conductivity measurement of a single diamond crystal provides a
primary step to the study of the transport properties in polycrystalline ﬁlms. In Chapter 3,
I introduce a more generalized technique for the evaluation of anisotropic conductivity
that involves the 3 x 3 tensor. Aided by another computer simulation (iterative
linearization technique), the method using FEM analysis determines the anisotropic
tensor in a complex geometry. This generalized technique has been tested for computer
models and for a real material (bismuth). The conditions for the application of this
technique are obtained from the analysis of the result. In Chapter 4, this numerical
technique (FEM) is introduced to study the electrical properties of a composite. A simple
composite in 2d, which consists of a host and two circular inclusions has been selected
for the demonstration of the method because some analytic solutions of the electric
conduction problem are available in this geometry. Two resistances (4-point and 2-point)
in the composite are deﬁned and obtained numerically and analytically over a wide range
of the neck distances between the inclusions. The numerical method turns out to

reproduce remarkably well the analytic solutions. In addition, a simple relation between

the resistances has been observed in the close neck limit. This relation may be useful in
the study of the effective conductivity of a composite in random close packing. The
numerical simulation using FEM could be considered a primary step in the study of the

effective medium theory for a composite material, especially in a complex geometry.

Chapter 2 CONDUCTIVITY MEASUREMENT OF
MICRON-SIZED DIAMOND CRYSTALLITES

2.1 Introduction

2.1.1 motivation

Because of its unique electrical and mechanical properties, diamond attracts much
interest these days. It has very wide band gap, high thermal conductivity, and extreme
strength as in Table 2.1. Motivated from these concerns, much research has been done on
the determination of the electrical properties of this material. Although it behaves
similarly in many respects to Si and Ge located in the same column of the periodic table,
it has some unique properties compared to these semiconductor materials. However, most
studies have been conducted on two-dimensional ﬁlms because they can be fabricated
easily and many important devices are in this form. A diamond ﬁlm grown on a Si
substrate (polycrystalline ﬁlm) is a mixture of many isolated single crystallites (see Figure

2.1(a)), it includes many local defects (grain boundaries, twins, and dislocations)

Table 2.1 Selected properties of diamond structure semiconductors.

 

 

 

 

Band Ga Lattice Parameter Bond Energy Thermal Conductivity
at 300K at 300K3 _1 4 at 300K
(CV) (A) (kJ-mol ) (W-cm'l-C'l)
Diamond 5.47 3.56683 347 24 ~ 255
Si 1.12 5.43095 196 1.53
Ge 0.66 5.64613 163 0.63

 

 

 

 

 

 

 

that contribute on the transport property of this ﬁlm. One of these defects in the ﬁlm is in
the form of a bicrystal in Figure 2.1(b). The electric transport in a Si polycrystalline ﬁlm
has been investigated and modeled by the Trap Transistor Model.6'7 It turns out that the
transport in a polycrystalline ﬁlm is dominated by grain boundaries. The role of these
local defects on transport in the diamond film has not been studied rigorously before,
because some electric properties of a single crystal are not measured simply due to the
complex geometry. Figure 2.2 shows an electron micrograph of typical single crystal that
is a well faceted micron-sized diamond.

Usually the determination of the conductivity of a regularly shaped macroscopic
sample is straightforward, both in principle and in practice. For a suitable choice of the
sample geometry, there is a region where the current ﬂow is uniform and, if the voltage is
sampled across this region the measurements are independent of the contact resistances [a
so-called 4-terminal measurement] and the resistivity is easily obtained. An example of

this is a cylinder in Figure 2.3 with current contacts near both ends and voltage probes

 

(a) (b)

Figure 2.1 A polycrystalline diamond ﬁlm (a) and a bicrystal (b) that was formed by two
single crystallites.

near the center of the bar. Since this cylindrical conductor has constant area through the
current direction, we have a simple relation for resistance R from the geometrical
parameters (L, A) by

L
R: _ .
pA, (21)

where L is length and A is area. If I neglect the contact resistances coming from the 2-
terrninal measurement, I need to use two extra contacts. The voltage between these two
leads is free of a voltage drop across the electrical contacts. However, this approach is not
possible for many samples for which regular shapes are difﬁcult to obtain or fabricate,
especially in the micron scale. There is a well-known 4-probe method to determine the
conductivity in complex geometry, as ﬁrst suggested by van der Pauw.8 However, this
method requires some restrictions on the geometry of the sample and contact pads. The
sample should be a quasi two-dimensional sample. The contacts should be on the
circumference of the sample. Therefore, the van der Pauw method cannot be used for the

3d diamond crystallites.

 

g); ZSKU 102.809 IP11 11019

Figure 2.2 A typical shape of CVD diamonds used in this measurement. It is a micron—
sized sample grown on a Si substrate.

 

 

Figure 2.3 The conductivity measurement of a simple shaped sample is usually
straightforward. This cylindrical geometry forms constant electric current
pattern inside the sample.

I suggest a technique based on the multi-probe measurement of I-V characteristics
and computer simulation using the ﬁnite element method (FEM). This method can be
applied to an arbitrary sample geometry under a restriction that the sample is Ohmic.
However, the contacts can be non-Ohmic. The diamond samples in this study have been
grown by a chemical vapor deposition (CVD) technique on a silicon substrate with boron
doping. This technique has been tested on a diamond ﬁlm, and the conductivity of
diamond crystallites has been determined. The conductivity of 3d microcrystallites has
been measured and compared with that of a 2d ﬁlm (especially polycrystalline ﬁlm).
Throughout this chapter, the conductivity of diamond is taken to be isotropic because of
its lattice symmetry. Although the determination by this technique has been done on
electrically isotropic diamond crystallites, it is a quite general scheme that can be applied
to other materials.

In Section 2.1.2, the van der Pauw method is brieﬂy reviewed and compared with
our technique. The general idea of the inverse problem and the ﬁnite element method

(FEM) will be described in Section 2.1.3 and Section 2.1.4. In the following Section 2.2, I

will introduce the resistance table from multi-probe measurements (Section 2.2.2) and the
geometrical factor will be deﬁned in Section 2.2.1 and the way to consider the contact
resistances will be described in Section 2.2.3. Section 2.3 is devoted to computer

simulations of image processing (Section 2.3.1), mesh generating and analysis using FEM
(Section 2.3.2) and the x2 ﬁtting technique to extract the material properties will be

described in Section 2.3.3. A brief description of the experimental measurements will
follow in Section 2.4. In Section 2.5, we tested this method with the van der Pauw
method on a 2d ﬁlm (Section 2.5.1) and then apply this technique to the data taken from
two diamond crystallites and discuss the results (Section 2.5.2). Finally, I summarize the

application of the method and give some suggestions in Section 2.6.

2.1.2 van der Pauw formula for 2d film

The determination of the resistivity of a regular shaped conductor is simple
because it provides a uniform current inside the sample so that a simple relation between
the resistivity and even the Hall coefﬁcient and parameters for the sample shape can be
given. However, this is not possible in case where a sample is grown and shaped
arbitrarily as in Figure 2.4. A method based on the well-known technique of conformal
mapping of two—dimensional ﬁelds9 to measure the resistivity and the Hall coefﬁcient of
2-dimensional ﬁlms was developed by van der Pauw.8 This method is very useful to
derive those quantities from simple I-V measurements and without considering contact
resistances between the sample and contacts by using multi-probe measurements.

Suppose we have a ﬂat sample as in Figure 2.4 and the current I ,2 is injected

 

Figure 2.4 The conductivity measurement of 2D ﬁlm in arbitrary geometry is done by
four—probe measurement.

from lead 1 to 2 and the voltage V34 between lead 3 and 4 is measured. Then a resistance
R13; is deﬁned by V34 / I12 and it is independent of the contact resistances because there is
no current flowing through the leads where the voltage is measured. The resistivity p of

the sample and 4- -point resistances (R12 , R3) are related by

 

crux-”R 34d)+ex x-p( ”R3331 —=) l (2.2)

where R5“ = V," / I u and d is the thickness of the sample. The resistivity p is given by

w R34+ R4l R34
p= 1112Mﬂ R13.) (2.3)

where f is a function of the ratio R1324 IR;I only and has a relation

R1324 —R;3] = fcosh"{exP(ln2/f)} ' (24)
12:; mg; 2

The Hall mobility can be obtained by the change of R1334 when a magnetic ﬁeld is applied

perpendicular to the plane. A resistance R1234 is given by a cross measurement of voltage

and current as done in the usual Hall geometry. If the magnetic induction is B and the

resistivity of the sample is known by the previous procedure,

_1AR13§‘

— B p 9 (2'5)

#H

where ARE: indicates a change of the resistance. Usually a clover-shaped sample , shown

in Figure 2.5(a), is used to eliminate the inﬂuence of the contacts. There are some
advantages with this shape compared to the bridge-shaped sample (Figure 2.5(b)). For
low electronic mobility, it provides more Hall current at same amount of heat dissipation.
Moreover, it has mechanical strength and can be used for small sample. However, this
method works only for a situation that satisﬁes the following conditions. First, the
contacts should be located on the circumference of the sample. Although there are some
marginally safe regions inside the sample, for instance near the comers of a square
sample, generally this condition should be satisﬁed.10 Secondly, the contacts should be

small compared to the sample size. This second condition raises trouble in micron-sized

   

(a) (b)

Figure 2.5 A clover-shaped ﬁlm (a) and a bridge-shaped ﬁlm (b) are used for four-probe
measurement.

samples. Moreover, the van der Pauw method requires homogeneous thickness of the
sample. Following a 2-dimensional approximation, it does not consider any current
variation along the perpendicular direction to the plane. Finally, the sample should have a
singly connected surface. There should not be any hole inside the sample, which
generates a separated boundary from the outer surface. With these restrictions, it does not

help us to calculate the conductivity of three—dimensional shaped microcrystallites.

2.1.3 inverse problem

All physical problems might be classiﬁed into one of two categories, the direct
problem and the inverse problem depending on what is known and what is to be
identiﬁed in the problem. The term inverse problem can be deﬁned simply as the
determination of the causes from the results. It is the opposite to the direct problem that
determines the results from the given causes (see Figure 2.6). The direct problem is the
determination of some physical quantities (potential, free energy, temperature, current
etc.) from the given conditions (geometry, boundary condition, initial condition, and
material property). Most analytic problems are usually direct problems. To solve a
problem analytically, all necessary conditions that cause some physical phenomena
should be deﬁned. Here is a simple example of direct problem. Suppose we have an
electrostatic problem to obtain voltage distribution in a spherical conductor that has
constant charge distribution. This can be done by solving the Poisson’s equation with the
prescribed conditions (geometry, charge distribution, etc.). Even some numerical

simulations are also direct problems. One of these examples is molecular dynamics.

11

 

 

CAUSES RESULTS

 

- Geometry Direct Problem - Potential

- Operator - Current

- Material Properties > - Displacement
(electric conductivity, - Temperature

 

A
‘

elastic constants,
thermal conductivity.)

- Boundary Conditions Inverse Problem

 

 

 

 

 

 

Figure 2.6 A schematic of the inverse and direct problem. Those two problems can be
characterized by what is known and what is determined in the problem.

Consider a system consisting of many particles in a ﬁxed volume. This system at a later
time can be described by the initial conditions (potential, initial position and momentum
of each particle, etc.). All positions and momenta of the particles at this later time can be
given by the simulation. However, some problems are quite different from the direct
problems. Suppose we determine some physical quantities from experimental data by
some ﬁtting procedure. The results (experimental data) are used to ﬁnd the causes
(physical properties). It can be called an inverse problem. The ﬁtting is usually performed
to identify the unknown causes from the measured results. The optimization technique

can be classiﬁed by an inverse problem. It determines some parameters that minimize the
measured quantities (energy, power, 12 etc.).

Suppose we have a heat conduction problem as in Figure 2.7. Given a temperature
distribution at an initial time t = 0 , ﬁnd the temperature distribution on the bar at later

time t. This problem is a direct problem because the solution of the

12

 

Figure 2.7 Heat conduction problem on a simple bar.

temperature distribution is determined by the initially given conditions. And this problem
is a well-posed direct problem because the temperature distribution at time t can be
determined uniquely from the given conditions (temperature distribution at t = 0 , heat

capacity cv of bar, and boundary conditions on the surface) and the heat conduction

equation governing the analysis.ll However, if we solve this problem in a different way
so that the initial temperature distribution on the bar should be found from the
temperature distribution at later time t, it is an inverse problem. This problem is an ill-
posed problem that does not have a unique solution because the temperature distribution

function T(x) in all regions of the bar cannot be determined from a finite number of

measurements of the temperature distribution at time t. Tikhonov developed a
regularization technique to solve a problem that cannot be solved by normal methods like
the inversion of linear equations.12 However, this regularization is not needed for our
problem that is a well-posed problem. Enough information is available to extract the
unknown parameters in our problem. We will determine the conductivity of the sample
and the contact resistances where the unknown parameters are fewer than the measured

(known) resistances provided by the multi-probe method. In our problem, the I-V

characteristics of the crystallite can be viewed as results and we identify the conductivity
and contact resistances as causes.

Suppose o is a ﬁeld variable representing a physical state of the system, which is a
response to an applied force f. Then a governing equation can be written as Eq. (2.6)

where L(k) is an operator and k denotes material properties of the system.

L(k)¢ = f . (2.6)
The inverse problem is, given the response 41, and the applied force f, to determine any
missing information involving geometry, operator form, boundary conditions and
material properties. Our problem can be more speciﬁcally classiﬁed as a material
property inverse problem following the classiﬁcation of Kubo.l3 He suggested that
inverse problems could be classiﬁed depending on what quantities to be determined.
There are ﬁve main inverse problems; domain/boundary problem (that includes electric
potential computer tomography (CT) method, inversion scheme for electric potential CT
method, crack identiﬁcation), governing equation inverse problem, initial/boundary value
inverse problem, force inverse problem, and material property inverse problem. In our
problem for the determination of the conductivity of a material with a complex geometry
by multi probe I-V measurement, the applied force f is the current and the response is the
potential. We know the form of the operator L deﬁned by the analysis type. The unknown

quantity that we need to identify in this problem is the conductivity as a material property.

14

2.1.4 Finite Element Method (FEM)

Initially the ﬁnite element method has been developed for the analysis of
structural mechanics because of its excellent applicability to complex geometry analysis.
It is a very efﬁcient way to analyze the mechanical property of an object with irregular
shape. However, it has been recognized it could be used to solve other problems like
electromagnetic, heat conduction, ﬂuid mechanics, acoustics, and piezoelectric studies.
This method can be extended to dynamic problems as well as static. There are numerical
methods such as ﬁnite difference method (FDM), Monte Carlo, spectral, variational
method for solving a partial differential equation. However, the ﬁnite element method is
considered the most efﬁcient for a study in a complicated and highly irregular geometry.
Only a few geometries of boundary value problem can be solved analytically even for
simple partial differential equations like Laplace’s equation. Such geometries might be
line, square, circle, sphere, and cylinder. However, most real problems generally have
very complicated shapes compared to these simple and ideal structures. Some problems
of irregular shapes might be solved by using the perturbation method. Otherwise, it
should be approached by successive numerical approximations.

Suppose we solve a boundary value problem as in Figure 2.8 and the inside
domain is governed by a partial differential equation

AI]! = f , (2.7)
where A is an operator that may represent Laplace’s equation, Poisson’s equation ,or heat
conduction problem depending on the type of analysis. f can be any type of source as
charge, heat source, and so on. The solution 11!, which represents voltage

(electromagnetic) or temperature (heat conduction), can be solved by well-deﬁned

15

    

SI

Figure 2.8 A boundary value problem is solved in this irregular geometry. Boundary
conditions on each surface are needed.

boundary conditions on the surface of the domain as

1;! = We : Dirichlet condition (2-8)
811,; = V’o' : Neumann condition , (29)
n

where div is a normal derivative of up. Since our problem is electrodynamic analysis of
n

a diamond crystallite, Laplace’s equation should be satisﬁed inside the region when there
is no static electric charge in the sample. We need Dirichlet conditions as boundary
conditions on the leads (SI, 82) where speciﬁc voltages are applied. Neumann condition

on the other surface S3 with I/lo'= 0 is given as no current ﬂows through the surface.

These boundary conditions can be represented by

VZV =0 (2.10)
V: V0 onSl,V=0 MS; (2.11)
ri-E=0 ons3, (2.12)

2.1.4.1 preprocessor

There are three main procedures in FEM, the preprocessor, analysis,
postprocessor. The preprocessor generates geometry and initial conditions of the
problem. Initial conditions consist of the construction of the geometry, meshes (nodes and
elements), boundary conditions, and material properties. In the procedure of geometry
reconstruction, we need to have the coordinates of vertices, faces and their connectivity.
We developed an algorithm to reconstruct a 3d sample from 2d SEM pictures.14 The
coordinates of the sample vertices obtained by this technique are used to generate the
geometry of the computer model in the FEM analysis. This process to build the computer
geometry is one of main steps in the preprocessor. Some commercial FEM packages
including AN SYS provide a CAD type utility for this process. The computer model
consists of lines in 1d, areas in 2d and volumes in 3d. The whole domain of the model
should be divided into subdomains if different materials are involved in the analysis.
Each material is assigned to a given subdomain.

After the construction of the geometry, this whole domain should be meshed into

/[\.

Figure 2.9 Some typical element shapes are presented for FEM. Each element has nodes
on the vertices or some extra nodes on the edges to enhance the accuracy.

 

 

 

 

l7

the elements. Each domain of the volume should be meshed separately. These elements
are the most important objects in FEM, where the analysis is performed inside. Each
element is deﬁned by nodes, lines, and facets as in Figure 2.9. Elements can be lines in
1d, triangles, squares, honeycombs in 2d, and tetrahedras, cubes in 3d. However, highly
distorted elements as in Figure 2.10 are not recommended because they may cause serious
numerical errors in simulation. Following the deﬁned shapes and sizes of elements, the
whole domain can be covered by the elements as shown in Figure 2.11. Apparently ﬁner
meshes can generate more realistic shape for a complex structure and can enhance the
accuracy of the simulation. However, many meshes increase degrees of freedom so that
the performance of the simulation would be limited by the computing power. A spherical
object is one of the most complicated structures to mesh well. All boundaries of elements
are deﬁned by a line segment in 2d, plane in 3d. A curve in 2d or a curved plane in 3d
cannot be used for the boundaries of the elements. A circle in 2d can be meshed only into
a polygon. It may cause a serious problem in the numerical analysis of two adjacent

circles. Except for these extreme structures, the meshing itself is not so

. « '5"

Figure 2.10 Highly distorted elements are not usually acceptable in mesh generating.

 

l8

 

Figure 2.11 A sample of meshed geometry with square-shaped elements in 2d.

important in the numerical simulations. The meshing generates the coordinates of all
nodes and deﬁnes the topology of each element deﬁned by these nodes. This information

is moved into the next stage of analysis.

2.1.4.2 analysis

In the analysis procedure, a function is deﬁned inside each element. As mentioned
earlier, since the physical quantities in the problem are represented only on the nodal
points at the boundaries of elements, they are not deﬁned yet inside the elements. We
express the physical quantities (voltage, displacements, etc.) in terms of the functions
whose magnitudes are deﬁned on the nodal points. This function is interpolated to
generate physical quantities inside the elements. This property is the main difference with
the finite difference method (FDM), which is very similar to the FEM in many aspects.
All variables, including ﬁrst and second derivatives of the physical quantities, can be

obtained by the interpolation of this function inside the elements.

The next equation is one example of this function deﬁned by the basis function

¢i~

q
‘l’ = X‘Pﬂli (2.13)
'=1

where (1),- is a known basis function and u, is an unknown coefﬁcient (Ritz
coefﬁcient).'5’16 We have q unknown coefﬁcients here. In FEM, the basis functions are

selected to calculate the physical values inside the elements. Then all values in the whole
domain can be represented by the coefﬁcients and the physical quantities on each node.
The global potential energy is represented by the sum of the local potential energy deﬁned
for each element. The variational principle is applied to ﬁnd out the values of won the
nodes. By solving for the coefﬁcients in the above equation, we can construct the
functional wand obtain the physical quantities in the whole region. One choice for the
basis functions is orthogonal polynomials. The degree of the polynomials is an option to
enhance the efﬁciency of the analysis. However, the increase of the degree of the
polynomials means the increase of degrees of freedom so that it causes a memory
problem in computing. We have a trade off here between the accuracy of the simulation
and the limit of computing power. This trouble happens similarly in the mesh generation
as mentioned in the previous section. The ﬁner mesh gives a better result, but it increases
the number of nodes where each functional tilis deﬁned. After setting up the functionals
of the physical variable, we solve the boundary value problem inside the elements. That
means we have the equations of the unknown coefﬁcients by solving the differential

equation and applying the boundary conditions between the elements and their adjacent

20

elements. The boundary conditions should be the continuity relations of the physical
quantities around the boundary of the elements. The ﬁrst derivative of the variable also
should be continuous. For instance, the voltage and electric ﬁeld on the boundary
between two ﬁnite elements should be continuous. Then, the potential energy of each
element is represented by the unknown coefﬁcients and combined together for the global

potential energy II as in
IT = J¢du 2*- 2(1)].(uij) . (2.14)
j=l

The variational principle should be applied to minimize the potential energy functional II.
The unknowns are the nodal point variables that are related to the unknown coefﬁcients
in Eq. (2.13). The global function 11 would be the electric power in electrodynamics,
potential energy in electrostatic, heat dissipation in heat conduction, and elastic potential
energy in the structural mechanics. In the variational method, the functional should satisfy
the minimization condition of

L“ = 0, i=l,..,n (2.15)
(911,.

where the approximate solution is given by

u" =ﬁf,a,. (2.16)

i=1
where j". is a Ritz function equivalent to a basis function and a, is a Ritz parameter. In
this procedure, the correct boundary conditions of u" should be applied on the surfaces

and the differential equation should be satisﬁed inside the complete domain. By solving

the minimization method, the physical quantities on all nodes and the unknown

21

coefﬁcients (Ritz parameters) are determined. Those values generate a complete solution

for the whole region.

2.1.4.3 postprocessor

In last step of postprocessor, a complete solution has been obtained for the whole
region. Since each variable on the node is given by the analysis, values inside the
elements can be determined by the interpolation of the functional in Eq. (2.13). From the
given quantities on the nodes, we can calculate other information. For instance, we solve
the electric conduction problem using only one variable that is voltage. However, after
the voltage on each node is obtained from the analysis, the current ﬂow can be calculated
from the electric ﬁeld and the conductivity of the material. The postprocessor provides
the visualization of the results by the interpolation technique; the contour map of voltage,
current, stress and strain distributions as presented in Figure 2.12. This is a very useful

and fast way to check if the simulation works well. More detail about FEM can be

 

Figure 2.12 By the FEM analysis, a voltage distribution on a sample is obtained. The
interpolation technique is used to obtain the distribution inside the
elements.

22

found in other references.” 18' '9

2.2 Theory

2.2.1 conductivity and geometrical factors

For a complicated geometry, the conductivity cannot be determined directly from
the measurement of the resistances since we measure the voltage and current that depend
on the geometry. We cannot separate the conductivity (geometry independent) from a
resistance (geometry dependent). Thus, we deﬁne a geometrical factor that only depends
on the geometry to treat the conductivity separately from the resistance. At this moment,
the contact resistances are neglected. The way to treat contact resistances (Ohmic and
non-Ohmic) will be described later in Section 2.2.3.

When we measure the voltage V and current I in a sample like Figure 2.13, a

resistance is deﬁned by

R=_=_,__'_,=ir=pr , (2.17)
. . S 0'

 

Figure 2.13 For 2-terminal measurement, the voltage is calculated by a line integral
along the path I and the current is obtained by a surface integral on the cross
section S.

23

where we assumed the conductivity tensor is isotropic; otherwise acan not be taken out
the integral. A newly deﬁned quantity I‘= I E - dl / IE ~ d8 in this relation depends on
the geometry only, not on the absolute value of the conductivity.

Suppose we have a sample with four probes as in Figure 2.14, 36 resistances can
be deﬁned from six measurements of voltage and six measurements of current. We can

represent each resistance with a corresponding quantity I‘ deﬁned in Eq. (2.17). From a
voltage between two leads I, m and a current through i, j, a resistance R5" is deﬁned in

the same manner.

 

jiiai jaw”
R5" =¥1=§ =—'j——=pr;" (2.18)
pas dyad“

We call the quantity 1‘5.“ deﬁned in Eq. (2.18) a geometrical factor. By deﬁning this

 

Figure 2.14 Showing a conductor that has four probes that are attached on the surface.

24

geometrical factor, the conductivity (which is a constant as an intrinsic material property)
can be separated from the resistance. This geometrical factor F is given by L / A for a
simple geometry of cylindrical shape, where L is a length and A is a constant area.
However, in most cases it can not be obtained by this simple form because of complicated
geometry that induces a complicated non-uniform current pattern. Following this
deﬁnition in Eq. (2.18), each resistance has its corresponding geometrical factor with the
conductivity as a scale factor. These geometrical factors can be given by a numerical
simulation described in Section 2.1.4 using unit conductivity. These calculated
geometrical factors will be used later to determine the conductivity by ﬁtting with

resistances given by the experimental measurement.

2.2.2 resistance table and its properties

In a multi-probe measurement, we can deﬁne many different conﬁgurations for
the resistance measurement depending on how to choose the leads for voltage and current
measurement. All possible measurements provide a set of resistances that forms a matrix,
a so called resistance table.

One of the diamond crystals (Sample 2) with contact pads prior to lithographic
lead attachment is presented in Figure 2.15. In a normal multi-probe measurement,
contacts are attached on the surface of a sample and the I-V characteristics are measured.
At least four contacts are needed to obtain contact-free resistances. Although we have ﬁve
pads on the surface in this sample, four of the ﬁve pads were used for the measurement,

with the furthest to the left having no lead attached. This pad was at a

25

 

Figure 2.15 Micron-sized diamond crystal grown by a chemical vapor deposition
technique before Au leads were attached. Five titanium contacts were
formed on the surface and four of them (indexed) were used for the
measurement.

constant (but initially unknown) potential that was determined by the ﬁnite element code.
With the four probes for the I-V measurement, a set of resistances can be given by all
possible voltage and current measurements. This set is called a resistance table that has
all necessary information of the sample conductivity and contact resistances. A current I 1.].
ﬂows through the terminal i, j and a potential difference Vb" is measured across terminal

I, m. We can deﬁne a (signed) resistance R5" = V," / I. which reduces to the standard

,1.
deﬁnition of resistance (2—terminal resistance) when i, j = l,m . With this deﬁnition of
resistance, we can generate a resistance matrix as Table 2.2 formed by V,"I column and I ,1.
row for all possible combinations. This set of measurement forms a "C2 (= n(n — l)/2) by
"C2 resistance matrix of all possible pairs of current terminals and voltage terminals but

there are only "C2 independent elements in this matrix because of the following

constraints. For an Ohmic sample, the matrix is symmetric under the exchange of current

and voltage terminals, Le.

26

Table 2.2 A resistance table from all possible combinations of voltage and current

 

 

 

 

 

 

 

 

 

 

 

measurements.
R g" Vl2 Vl3 V14 V23 V24 V3,,
1'2 R1122 R1123 R112
1‘3 R133 R11;
1” R11:
’23 R333
1,, R23:
1,, R33:

 

 

 

 

 

 

 

27

 

i’ I" m V
Rd. =73 = R;- = J- (2.19)

m
9
1m [if

where V0. is the voltage between the leads i, j and I ,m is the current between the leads I,

m. This is a reciprocity theorem,20 which is proved in the Appendix A.1 for the present

geometry. Along any row of the resistance matrix, the current I ,m is a constant, and the
voltage is a scalar potential that satisﬁes
%=%+%. em

Thus, the resistance elements have the additive property

R1} = vi)“ = Vr'k +ij
Im II

I = Riii. + Riff. . (2.21)

 

lm m

Using these two rules, we can show that there are only ”C2 independent elements in the
resistance matrix by the following construction:
Consider the n —1 elements (R13; , i = 2,..,n) in the ﬁrst row of the resistance

matrix. These terms can generate all the other elements in this row using the addition rule
Rg=Rg—Rg, am)

where i < j = 2,..,n. Note that changing the order of the indices changes the sign of the

resistance. In the second row, only n — 2 elements of R1]; (i = 3,.., n) are needed to

construct the remaining elements in that row because, by the reciprocity theorem,

R}; = R3 and R3 was given in the ﬁrst row. Continuing in this manner we need n —1

terms in the ﬁrst row, n - 2 terms in the second row, for a total of

28

n(n — 1)

 

(n-1)+(n—2) +....+2+1 = ="C2 (2.23)

independent elements. There are many possible choices of the "C2 independent elements,

but they cannot be selected arbitrarily. By the above construction, the (n — 1) by (n — 1)

upper diagonal block forms one independent set. Surprisingly, the "C2 diagonal elements
(the set of 2 terminal elements) also form an independent set, related to the off diagonal
elements by

jl im _ jm _ 1'!
R jl + Rim R jm Ril

lm _
RU. _ 2

 

(2.24)

This relation means that knowledge of only the 2-terminal measurements can be
used to generate the complete resistance matrix, which can be used as a useful self
consistency-check on the experimental results. Table 2.3 shows how the diagonal terms
reconstruct the whole resistance table. Three off-diagonal terms in the triangular matrix in
the bold box can be obtained from the six diagonal terms by Eq. (2.24). Six terms in the
bold box was already shown to generate the whole matrix in Table 2.3. In Table 2.4, I
summarized how the number of elements of the resistance matrix, and the number of
independent elements depend on the number of terminals. We need at least four terminals
to extract the resistivity p and the contact resistances because otherwise there are more
unknowns than independent terms.

As an example, consider the case for four terminals (n = 4) shown in Table 2.4.
There are six ways to select two of four terminals and so there is a 6 by 6 resistance

matrix with six independent elements. The six diagonal elements form the most

29

Table 2.3 How to generate the triangular matrix from six diagonal terms so that the
triangular matrix produces a whole table

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

R3" VIZ VIS V14 V23 V24 V34
1.. 11:; R3 +14: we 11:: +14: 42::
2 2
[13 R11; Rii 1' R1]: " R3:
2
Ii. R13:
123 R 2233
12. R22:
1,, R332
Table 2.4 The number of independent elements in a resistance table.
Number of Dimension of Elements in Independent Unknowns
Terminal Resistance Table Resistance Elements ( p , contact
Table resistances )
2 l 1 1 3
3 3 9 3 4
4 6 36 6 5
5 10 100 10 6
6 15 225 15 7
" "c2 =n(n—1)/2 ["c2]2 "C, =n(n—l)/2 n+1

 

 

 

 

 

3O

 

 

important set of independent elements from which we can use Eq. (2.24) to obtain the
three off-diagonal terms, R1123 , R11; , R11; as shown in Table 2.3. The 3 x 3 upper

diagonal block (shown inside the bold line in Table 2.2 and Table 2.3) is an alternative
independent set that can be used to construct the entire matrix using the procedure
outlined above. However, it is interesting to note that the set of 4-terminal measurements
is not all independent. Note that by Eq. (2.24), only two of the 4-terminal elements are
independent with
R3; = R3: + RS; . (2.25)

As a result, a set of six terms including three 4-terminal measurements cannot be used to
reconstruct the whole resistance table.

It is important to note that the above analysis is only true for samples where the I-
V curve is independent of the current direction (otherwise the reciprocity theorem fails).20'
2‘ Some of the experimental data analyzed in Section 2.5.3 has non-Ohmic characteristics
indicated by asymmetry of a table, and so the resistance matrix has more than 6
independent elements. The necessary modiﬁcations to the above analysis are model

dependent, and will be discussed for the speciﬁc situation of the experiment.

2.2.3 contact resistances

In two previous sections, it has been found that the geometrical factors separate
the unknown conductivity from a resistance. However, the contact resistances should be
considered in 2-terminal and 3-terminal resistances. The measured resistances include the

contact resistances, so the representation of the resistances should be modiﬁed. We

31

 

Figure 2.16 Contact resistances have two different values (r, , r,. ') for non-Ohmic

contacts. A resistance measurement consists of a sample resistance Rum
and contact resistances added in series.

Table 2.5 Resistance table with two contact resistances for current ﬂows.

 

 

 

 

 

 

 

 

Riljm Viz VI} VI4 V23 V24 V34

112 pl“),2 + rl + rz' pl‘1‘23 + rl p131;t + r1 P1753 - 72' pl": ” ’2' P 1324

1,3 p11; + ’1 pi]? + r1 + r3' p11: + rl p113,3 + r3' pI‘lza4 p17,? — r3'
1,, p113 + r. pFi'i +ri PW: + n + '4' p11? p113.“ + r.‘ pF.3.“ + r.‘
123 prz‘i - r. p1“); + r,‘ tori; PF2? + '2 + '2' 10132: + 4 p133: - r.‘
124 pm? - r2 pFi‘i‘ pFlZ + r.‘ p133} + r. PF2? + '2 + '4' pl"; + 0'
I34 p133} p133,3 — r3 p1“;3 + r,‘ pl‘ff — r3 p133: + r4' p133,4 + r3 + r,‘

 

 

 

 

 

 

 

32

 

assume these contact resistances can be added in series to the sample resistance because
the contact area is assumed an equipotential due to the high conductivity of the Au/T i

contact pads. This leads to

Rim =priim+’i'(5u"5im)+rj(5jm—6ﬂ) , (2.26)
where 1;. is the Ohmic contact resistance of terminal i and 5,]. is Kronecker delta function.

However, contact resistances are generally non-Ohmic which is indicated by non-

symmetric three terminal elements in the observed resistance table. We assume that the
sample itself is Ohmic, which is justiﬁed if the four terminal elements are symmetric.
Since the I-V characteristics are non-linear and not symmetric for non-Ohmic behavior of
the contacts, each contact has two independent resistances for the two different current
directions [i.e. current reversal with the magnitude of the current held constant]. To
investigate such non-Ohmic I-V characteristics of a contact, we used two ﬁtting

parameters, 4, ’1' corresponding to the two contact resistances for constant current ﬂow
as in Figure 2.16. If r, is a contact resistance when a current is ﬂowing into the sample
and '1' is for a current ﬂowing out of the sample, then Eq. (2.26) is modiﬁed to become

R5?" = p13.“ + r, (6,, — 6m) + 5.16),, — 6],) (2.27)
to be a more general expression for non-Ohmic contacts. The whole expression of the

resistance table following the relation of Eq. (2.27) is given in Table 2.5.

33

2.3 Computer Simulation

2.3.1 reconstruction of 3d computer image

In FEM, all the information of the sample geometry is needed for the
preprocessor that constructs a computer model for the analysis. It may include the
coordinates and topology of all vertices, lines, planes, and volumes associated with the
sample shape. However, since our samples are micron-sized crystals, it is not
straightforward to measure the shape of such small objects. One way is using electron
lithography to get the SEM (Scanning Electron Microscope) photos. An inverse ray
tracing technique is used to reconstruct the three dimensional object that will be used in
simulations. Ray tracing is a computer process used to create a two dimensional
representation of a three dimensional object by tracing the path of the ray of light between
the light source, objects, and an observer.22 We apply this technique in reverse to
construct a three dimensional image from the two-dimensional SEM photos. Since these

photos were only two dimensional, we needed several pictures taken along different

viewing angle

 

sample y

 

X

 

 

Figure 2.17 Each SEM photo is given in an orientation for viewing angle (6, ¢) as
shown in this ﬁgure.

34

 

 

(C) (d)

Figure 2.18 SEM photos of a sample (Sample 2) in various orientations. (a) 6, ¢ = 0° ,
0° (b) 0, ¢ = 50° , 0° (c) 0, ¢ = 50° , 90° (d) 9, ¢ = 50° , 270°

35

orientations deﬁned by viewing angles in Figure 2.17. Figure 2.18(a)~(d) present four
different pictures projected in each direction before gold leads are attached on the ﬁve
contacts formed on the surface. In this sample, contacts were located on the upper surface
since it has been formed by the electron beam along the vertical direction. Therefore, the
picture in Figure 2.l8(a) taken along the same direction shows most of the detailed

geometry. However, some contacts were observed smeared down the surface (see Figure in

2.18(c)). This contact formed a wide contact area that could affect the I-V measurements.
For a more complicated contact formation or some serious deformations in the lower part,

it is necessary to select all viewing angles carefully to take the whole geometry of the

 

‘1‘ -..- m.‘ m

sample. After taking these SEM pictures, the coordinates of all vertices of the sample and
contacts were digitized on each photo.

The algorithm treats the crystal as an irregular polyhedron and attempts to assign
coordinates to the vertices in a way that is consistent with the SEM images. From each
picture of the crystal we can measure the projected length of an edge of the diamond
crystal, and compare it with the projected length calculated from the assumed vertex

coordinates. The vertex coordinates are adjusted to minimize

 

, 1.4.55“ '2
x, =2[' 6? ]. (2.28)

where ll. is a calculated projected length, If” is the projected length from the SEM and

o". is an estimate of the experimental error in measuring [1.35M . The sum is over all edges

visible in each picture, and over all pictures. The positions of the contact pads are
determined by including the coordinates of the comers of the contact pads, and the

lengths from the comers to some of the crystal vertices. Taking pictures from enough

36

different orientations to ensure that every edge was visible in at least two images was
adequate to ensure that we could construct a representation of the sample that matched
the SEM images very well. Because we knew that the sample was a single crystal, we
added some additional constraints to the minimization process. Speciﬁcally we required
the facets to be planar, and the facet angles to conform to the known facet angles of the

diamond structure, e. g. 70.5° or 109.4° between different (111) facets. The ﬂatness

constraint was achieved by adding to 12 the term

22? (r —-—c-‘)2], (2.29)

where if). denotes j—th vertex coordinate in i-th facet and E. is a constant vector that is

normal to the i-th plane. The facet angles are constrained by adding to 12 a term
6 _0 theory 2
1,3 = [.-_.-__:| , (2.30)

where the facet angles 9,.”"0'3’ are those appropriate for a diamond structure. The sum is

over pairs of facets. The 0",.” in Eq. (2.29) and of in Eq. (2.30) are standard deviations

denoting the weighting factors for each constraint. Although it is arbitrary how these
factors are assigned, they should be larger than the standard deviation associated with the

edge lengths I, because these lengths could be determined very easily and directly from
the SEM photographs.

12 =x,2+xp2+xf2. (2.31)
The total 12 , the sum of the three quantities in Eq. (2.28)~ (2.30) is minimized by

adjusting the vertex coordinates using a standard non-linear minimization routine.23 The

37

ﬁnal coordinates associated with the minimized 152 were used as the input for the current
ﬂow analysis that was done using the ﬁnite element method.

After this optimization procedure, we obtained a computer model of the sample. It
provided all information about the geometry deﬁned by the coordinates of vertices and
the connectivity of planes and volumes. Nearly 80 vertices were needed in our models
and most vertices were for the contacts. The contact areas have been constructed very
carefully since the transport in the sample mostly depends on the area. By initially
assigning high conductivity to the contacts, the contribution of the contacts to the
resistance measurements was determined only by the effective area for the current ﬂow.
We attached ﬁve contact pads in the shape of small pyramids on the sample as in Figure
2.19. Although only four contacts of ﬁve were used for the real measurement, we made
the ﬁfth contact also to generate the equipotential area below the contact. Figure 2.19(a)

and (b) shows two views of the computer model (Sample 2), which correspond to the

 

(a) (b)

Figure 2.19 Computer images reconstructed from SEM photos by image processing.
Five small pyramids on the surface denote ﬁve contacts. Compare these two
images with the SEM photos in the Figure 2.18.

38

SEM image in Figure 2.18(a) and ((1) respectively. It was not possible to observe the
bottom of a sample because of the Si substrate. Some diamond crystallites may have an
unusual structure of local defect below that cannot be identiﬁed in the SEM photos.
However, usually it does not contribute signiﬁcantly to the transport as much as the top
part of the sample does. On the surface near the contact pads, a large variation of current
ﬂow was observed. However, it was not so on the surface below where no contact pads

were formed.

2.3.2 modeling and analysis using FEM

2.3.2.1 modeling (preprocessor : ANSYS)

After the geometry reconstruction, I did the simulation using the ﬁnite element
method to solve the electrical conduction problem. All FEM tasks in this problem were
performed by the commercial packages AN SYS24 and ABAQUS”. ANSYS 6.1 was used
for a preprocessor (mesh generating) and the other processes (numerical analysis and post
processor) were done by ABAQUS 5.5. In modeling by AN SYS, I generated volumes that
correspond to the sample and contacts. Since our sample had highly irregular geometry, it
was not straightforward to generate the structure. Some simple structures (squares,
circles, bricks, spheres, etc.) can be constructed by giving a few coordinates and
geometrical parameters. First, it begins with the formation of keypoints for all vertices of
the sample and contacts. Planes for the volume construction are made by linking the
keypoints that are associated with a plane. Some vertices of a plane may not be located on
the same plane even after the geometry optimization with the constraint of a ﬂat plane.

AN SYS does not allow making a plane by linking these vertices. All keypoints should be

39

 

(b) (C)

Figure 2.20 Triangular planes generated from three keypoints were combined into a
large plane as shown in (a). A volume could be constructed by connecting
these large planes into a volume unit. The volume elements for the diamond
crystallites were presented in (b) and (c).

40

on the same plane to construct a plane. Thus, a triangle plane was used to overcome this
trouble. Suppose we need to make a plane structure by four vertices (1, 2, 3, 4) in Figure
2.20(a) and the vertex 4 is not on the plane deﬁned by the other three vertices (1 , 2, 3).
Then we make two triangular planes and connect them to form a folded plane made of
four vertices. Most planes in our samples were constructed in this way as shown in Figure
2.20(b) and (c). A volume was constructed by combining these planes. Figure 2.20(b) and
(c) present the structured volumes of two diamond crystallites before the meshing. Both
of our samples consisted of six separated volumes (one huge volume for the sample and
ﬁve small pyramids for the contact pads). Different material properties could be assigned
to the separated volumes in the FEM analysis.

After constructing these volumes, we need to select the type of element for the
meshing. Element types are selected depending on the dimension and the analysis type.
Figure 2.21 shows the element shapes used in the meshing of the diamond ﬁlm (2d) and
the crystallites (3d). For the analysis of electrical conduction in 2d, AN SYS provides
triangular elements with three or six nodes and squares with four nodes or eight nodes. In
3d, we can use cubic elements with 8 or 20 nodes and tetrahedral elements with 4 or 10
nodes. Tetrahedral elements were usually used to mesh an irregular shape in 3d where it
was not done by cubic elements. Tetrahedral elements with four nodes were used for the
diamond ﬁlms and one diamond crystallite (Sample 1) and tetrahedral elements with 10
nodes were used for the other diamond crystal (Sample 2). The diamond ﬁlm was meshed
by square elements (2d) as well as tetrahedral elements (3d) to test whether the 2d
approximation of the ﬁlm was valid in the simulation. Figure 2.22 shows the meshed

structures of all samples; a diamond ﬁlm(a) and crystallites ((b) and (c)).

41

Figure 2.21 The elements used for the meshing of 2d ﬁlm (square) and 3d crystallites
(tetrahedra). Sometimes elements that have only four nodes on the comers
were used.

Table 2.6 The characteristics of the meshings for each sample.

 

 

 

 

 

 

Homoepitaxial ﬁlms Sample 1 Sample 2
Element type DC3D4 DC3D4 DC3D 10
Nodes per elements 4 4 10
Total Nodes 800 1100 4700
Number of elements 2300 5500 3300

 

 

 

 

42

 

 

(b)

 

(C)

Figure 2.22 Meshed structures of the diamond film(a) and crystallites (Sample 1(b),
Sample 2 (c)). Figure 2.20(b) and (c). Each model has square-shaped
contacts on the surface.

2.3.2.2 analysis (processor : ABAQUS)

The meshed structures given by ANSYS were delivered to ABAQUS for the
analysis. All coordinates of nodes and the topology of elements were used to generate an
ABAQUS input ﬁle that includes all necessary information for the analysis; material
property, boundary conditions, and output format as well as the geometry. Since
ABAQUS does not provide the electromagnetic analysis, the heat conduction analysis
was used. Both analyses follow the same differential equation, Laplace’s equation with no
charge and heat source. The temperature in the heat conduction analysis is equivalent to
the voltage in the electric conduction analysis. The heat ﬂux and the thermal conductivity
correspond to the current and the electric conductivity respectively. By these
correspondences, the electric conduction problem could be solved by the heat conduction
analysis using ABAQUS. In the heat conduction analysis, the temperature is the only
degree of freedom that is deﬁned on each node. By minimizing the heat dissipation
deﬁned by the temperature, the solution on each node could be obtained under the
subscribed boundary conditions and the material properties. After the temperatures on
every node were obtained, the heat ﬂux could be determined from the derivatives of the
temperature functions and the thermal conductivity. An electric resistance could be
deﬁned by the temperature (voltage) and the heat ﬂux (current). Some details of these
procedures were described earlier in Section 2.1.4.2 and the rest can be found in some
other references or ABAQUS manuals.25

Our main object in this analysis is to obtain the geometrical factors, not
resistances. They represent the shape factors that are not related to the intrinsic material

property (conductivity). To obtain these factors in this analysis, unit conductivity is given

44

to the conductor. Then, the resistances calculated from the voltage and current given are

equivalent to the geometrical factors Fiﬁ.” . In the absence of any contact resistances, these

factors are related to the resistances R5" by a simple linear relation of

lm _ lm
RU — p11,. . (2.32)
All possible combinations of the current and voltage measurements generate a 6 by 6

geometrical factor table. These geometrical factors include any information of complex

current paths across the arbitrarily shaped sample for each combination.

2.3.2.3 contour map (postprocessor : ABAQUS)

From the analysis in the previous section, the voltages only on the nodes were
obtained. Inside the elements, the equipotential lines were calculated from the
interpolated functions of the voltage. The contour maps of the voltage distributions on the
sample surfaces are shown in Figure 2.23 for a selected conﬁguration of the current and
voltage measurement. The current ﬂows from the red pad (high potential) to the blue pad
(low potential). The voltage distributions are indicated by the contour map. Figure 2.23(a)
presents the potential distribution on the surface of a ﬁlm. Figure 2.23(b) and (c) present
the potential distribution of a crystal (Sample I) viewed by different orientations. Figure
2.23(d) and (e) show the similar contour maps of another crystal (Sample 2). Figure
2.23(b) and (d) were viewed by the vertical direction (view point =(0, 0, 1)) and Figure
2.23(c) and (d) were rotated and shown by a different View point (1, -1, 1). By assigning

high conductivity ( > 106 ) to the contacts and unit conductivity to the sample, the contour

45

 

(d) (6)

Figure 2.23 Equipotential lines on the surface of the samples given by ABAQUS. (a) is
for homoepitaxial ﬁlm. (b) (c) are for Sample 1. (d),(e) are for Sample 2.
The ﬁgures of 3d samples were obtained in different orientations.

46

maps show the constant potential distributions on all contact areas including the lead-free

contact as expected.

2.3.3 least square fitting
Table 2.5 in Section 2.2.3 is a modiﬁed resistance table for a 4-terminal problem

showing each element in terms of the sample resistivity p, the geometric factors 1“,?" and

6 of the 8 contact resistances (rl , r2 , r2 ' , r3 , r3' , r4 ' ).26 Another table for reverse current
ﬂows is necessary to extract all the contact resistances for non-Ohmic contacts. From

these two tables given by opposite current directions, a resistance can be deﬁned as

R1?" = P135" +7.-(5.-i JANE-(5).. —5,-i)+Ari(5.-i *5im)-Ar,-(5,-m -5,-i) (2.33)

where
_ l ,
and
1 .
An=-2-n—n). (2.35)

A similar resistance table with both current and voltage leads reversed enables us

to construct an Ohmic resistance table RE“ = (R53 + Rf,“ )/ 2 which is shown in Table 2.7,

and the non-Ohmic resistance table M3?” = (R? — R7,.“ )/ 2 , shown in Table 2.8. Note that

the Ohmic resistance table displays all the symmetries of the analysis in Section 2.5.2,
and in particular, only the six diagonal (2-terminal) elements are needed to generate the

entire table. For the non-Ohmic resistance table, the six diagonal elements only

47

Table 2.7 The Ohmic resistance table obtained from Table 2.5. (symmetric)

 

 

 

 

 

 

 

 

 

 

 

 

 

 

f?- Vrz VI3 Vi4 V23 V24 V34
112 Frizz + 7'1 + 72 131-1'23 '3' 71 [71.1124 + 7"] 13131223 - 7'2 Frizz4 _ 72 1,11324
1.. 1111?”. pr}: +4 +4 pri‘: +F. 10113.3 + F. p113.“ p113: -F.
I” p F1142 + Fl PIT: +71 Fri: + F1 '1’ F4 pr}? [713114 + 74 [31—1344 + F4

12 — 13 — 14 23 - — 24 - 34 -
123 Przs “’2 pF23 +"3 Prza 131-‘23 +"2 +"3 przs + '2 przs ’ ’3
I 2. p133 - F. pri'i p1“): + r, prjf + 7, m3: + F. + F. pl‘i’: + F.

12 13 — 14 — 23 — 24 — 34 - -

134 Pr34 Pry —r3 pram +"4 p134 -r3 pr” +"4 p134 "H3 +r‘

 

Table 2.8 The non-Ohmic resistance table obtained from Table 2.5. (non-symmetric)

 

 

 

 

 

 

 

 

 

 

 

 

 

 

AR Viz Vrs Vr4 V23 V24 V34

I ,2 Ar. - Ar; Ar; Ar; Ar; Ar; 0

I ,3 Ar] Ar] - Ar3 Ar] - An 0 Ar;

I 1 4 Ar] Ar] Ar] - Am 0 - Ar4 - Am
123 - Arz - An; 0 Ar; - Ar3 Ar; Am

I 24 - Arz 0 - Ar4 Ar; Ar; - Ar4 - Am

I 34 0 - Ar3 - Am An - Am An - Ari;

 

48

 

 

determine the off diagonal elements up to an arbitrary constant. To determine all the Ar; it
is necessary to have at least one of the 3-terminal elements. For this model, the non-
Ohmic resistance table is independent of the resistivity, because we have assumed that the
material itself is Ohmic.

Because the geometric factors are known from the FEM, we can in principle
obtain the resistivity and the contact resistances by doing a least squares ﬁt for the

expressions intheoretical resistances (Table 2.7 and Table 2.8) and the measured

resistances; i.e. we minimize a x2 deﬁned by

 

2
X2=2|:Ri —Ri'(p9r):l , (236)

i of

where R f are the experimental measurements and R} (p, r) are the theoretical

resistances deﬁned by Eq. (2.33) as a function of the resistivity p and contact resistances

r. The sum is over the elements of the resistance matrix used for ﬁtting and o". is a

standard deviation of each experimental measurement, which includes both errors in the
resistance measurements. The errors in the geometric factors might be from the incorrect
computer modeling of the geometry as well as the small numerical error in the FEM. In
practice, the least squares ﬁt is best done separately for the Ohmic and the non-Ohmic
parts (Table 2.7 and Table 2.8). The practical details of this ﬁtting will be discussed in the

next section, when it is applied to our two samples.

49

2.4 Experiment

2.4.1 sample preparation

In this section, several experimental processes of sample preparation, contact
formation, and the I-V measurement by multi probes will be brieﬂy described. Details of
this work can be found in J aeger et al.27 The experimental measurements have been done
on homoepitaxial diamond ﬁlms and microcrystallites. The ﬁlms were grown on the
diamond substrates where as the microcrystals were grown on Si substrate. Since the
conductivity of diamond crystallites has not been measured before in highly irregular
geometry, our technique should be compared with a well-known method. We chose the
van der Pauw method on the conductivity measurement of a 2d ﬁlm. The conductivity of
a 2d ﬁlm can be given by the van der Pauw formula in Eq. (2.1), thus our result can be
compared with the result from this formula.

The size of a homoepitaxial ﬁlm (Figure 2.24) is 3 x 3 mm2 and its thickness is
1.86 um. It is a good two-dimensional object. The I-V measurement was done twice on
the same film with different contact formations. Once it was done with a ﬁlm with four
contacts as in Figure 2.24(b). Later two more contacts were attached on the same sample
as in Figure 2.24(c) and the measurement was done again. However, only four contacts
(gray pads in the ﬁgures) were used for the measurement in both cases. Figure 2.25
presents two SEM photos of the CVD diamond crystallites whose conductivities were
measured by our method. Although both samples had ﬁve contacts, only four contacts
were used for the resistance measurement. Two crystallites were grown in the same

condition. However, a twin boundary between the contacts was observed in one of the

50

 

 

 

 

 

 

 

1:1 1:1 1: E1
2 3 2 3
1:11 4121 1:1.1 4-
(a) (b) (C)

Figure 2.24 A normal picture of the 2d sample (diamond homoepitaxial ﬁlm) is shown
in (a). Two schematics ((b), (c)) of the sample denote the positions of the
contacts. Film 1 (b) has four contacts and ﬁlm 2 (c) has six contacts. Gray—

colored squares in the ﬁgures represent the contacts actually used for the
measurement.

 

(a) (b)

Figure 2.25 Two micron-sized CVD diamond crystallites, Sample 1 (a) and Sample 2

(b). Sample 1 has Au leads attached on the Ti contacts and Sample 2 is a
shape before the lead attachment.

51

Table 2.9 The characteristics of CVD diamond samples.

 

 

 

 

 

Sample Description
Diamond ﬁlm Homoepitaxial ﬁlm on single crystal (100).
(homoepitaxial) Natural diamond substrate.
Dimension : 3.0mm x 3.0mm x 1.86pm
Contact size : 0.3mm x 0.3mm
Sample 1 Isolated microcrystal with no visible defects.
(microcrystal) Ohmic contacts.
8102 substrate.
Diameter : ~ 4 mm , Height : ~ 2 pm
Sample 2 Isolated microcrystal with a twin defect.
(microcrystal) Non-Ohmic contacts.

 

SiOz substrate.

Diameter : ~ 4 um , Height : ~ 2 um

 

52

 

sample (b). Although no signiﬁcant difference of the contacts was observed in the photos,
it turned out that the contacts of one sample (a) were Ohmic where as the other contacts
in (b) were non-Ohmic. Our method could be tested more generally for these two
different samples whether it determines the conductivity, which is independent of the
contact properties. The characteristics of all the samples in this work are listed in Table
2.9. All diamond samples were prepared by microwave plasma enhanced CVD at Kobe

Steel USA, Inc.

2.4.2 four-probe measurement

With four probes formed on the samples, the I-V characteristics were measured by
DC voltage as a function of temperature in a vacuum level less than 1045 torr. To
investigate the I-V characteristics of the contacts, a ﬁxed amount of current was used for
each current ﬂow. All possible combinations of the voltage and current measurement
generated a 6 by 6 resistance table that was described in Section 2.2.2. A total of 36
elements in the resistance table consist of six 2-terminal terms, six 4-terminal terms, and
twenty-four 3-terminal terms. While these measurements were done through a broad
temperature regime (-100 ~ 170° C) as in Figure 2.26, a complete measurement of a

resistance table was done only for a speciﬁc temperature for each sample. At the other
temperatures, only one 4-terminal resistance (R;4 = V14 / I 23 ) was measured. Since the

geometrical factor is temperature independent, the conductivity at the other temperatures
can be obtained by dividing the 4-point resistances with this geometrical factor. The

investigation of the temperature dependence of the conductivity provides the activation

53

 

 

 

 

 

 

 

 

 

1000001 ' 3 F ' Y ' C ’ r ' ' 1

1 1

J o d

1 ° 3,, 3

° 1

j <3, . MC-A (lavu) ﬂ

oo o MC-B (“by“)

10000 1 o 1

1 o 1

1 .

“ 1

1000 -. g

1 l

A . 4

G 4 4
x

V g :1
O‘.’

100 1 1

= i

I 4

d J

l l

10 1 "‘

1 i

2 J

1 1
V I j ' ﬁ' 1 V r ﬂ T I I '..__.
-150 -100 -50 0 50 100 150 200
T (°C)

Figure 2.26 Temperature dependence of a resistance of diamond crystallites (Sample 1
(solid circles) and Sample 2 (open circles)). To investigate the temperature
dependence of the resistivity of these samples, a four-terminal resistance

(R;3 ) were measured over a wide range of temperatures.

54

energies of the doped semiconductors.28 The conduction mechanism of each sample can
be understood and compared using these activation energies.

The resistance measurements for some samples were done with two different
current directions to determine the non-Ohmic contact resistances. Usually two diagonal
resistances were measured for these two current directions; they could provide complete
information of non-Ohmic contact resistances. In this submicron scale, the leads were
electrically fragile because static discharges could happen near the contacts. The melting
of metal leads that might cause the missing of an entire portion of the contacts sometimes
happened during the measurement. So the current level should be kept below 10 11A to
prevent the contact disconnection from the melting. With this problem, some
measurements did not provide enough information (two complete sets of resistance table

for two current ﬂows) needed for the analysis.

2.5 Analysis

2.5.1 diamond film

To test our method, we used two films with multi probes whose conductivity
could be determined by the van der Pauw method. One (Film 1) had only four leads on
the ﬁlm and the other (Film 2) had six leads. With the Film 1 with four contact pads of
Figure 2.24(b), we measured all possible combinations of current and voltage
measurement to construct a 6 by 6 resistance table of Table 2.10. It was given by constant
current (1 (IA) to investigate Ohmic behaviors of contacts at T=24.7° C. For this sample,

we measured two diagonal terms (two terminal measurements) for two different current

55

Table 2.10 Resistances measured for homoepitaxial ﬁlm by 4-probe experiment (k9).

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

R VI2 V13 Vl4 V23 V24 V34
I12 264.6(329.5) 106.65 90.95 -158. 17 -173.94 -_l_5_.'_72
I13 107.41 l98.7(l98.7) 105.988 91.239 1.43 -92.658
I14 91.62 106.097 l84.5(184.5) M 92.636 78.158
I 23 -232.6 90.23 M 322.4(248.1) 247.3 -75.834
I 24 -246.7 1.2—12 91.391 246.5 337.7(265.1) 92.612
I34 ;_1_5i9_ 93.7 14 78.268 -77.86 94.170 l72.l(l72.l)
Table 2.11 Geometrical factor table obtained by FEM (mm'l).

F Vl2 V13 V14 V23 V24 V34
I12 904.04 515.40 393.80 -388.65 -510.24 -121.59
Il3 515.32 1023.48 506.32 508.16 M -517.16
114 393.80 506.27 947.54 _l_l_2._42 553.74 441.27
I 23 -388.68 508.13 M6 896.80 501.13 -395.67
124 -510.18 M 553.81 501.14 1063.99 562.85
I 34 M -517.22 441.30 -395 .68 562.84 958.52

 

 

 

 

 

 

 

Table 2.12 The resistivity and contact resistances of diamond ﬁlm (Film 1) given by van
der Pauw formula (V DP) and our method.

 

 

 

 

P (9 cm) P (52 cm) 71 . Ar 1 (m) 72 . Arz (kg) 7'3 . A73 (k9) 7.. . Ar4 (kn)
(VDP)
12.51 12.91 37.8, -0.9 139.8, 34.1 24.55, -1.45 19.95, -l.35

 

 

 

 

 

 

56

 

 

directions to identify all non—Ohmic contact resistances. Notice that some of the_matrix
terms showed substantial asymmetry indicating that some contacts were non-Ohmic.
However, the 4-terminal measurements (underlined) showed a relatively good symmetry.
We also did the FEM analysis using ABAQUS to ﬁnd out the geometrical factor table in
Table 2.11. From the 36 elements in this table, we used only the 4-terminal terms to
determine the resistivity of this film that would be compared with the same quantity
calculated from the van der Pauw method. The contact resistances were obtained by
ﬁtting the diagonal elements in both tables. They turned out substantially non-Ohmic as
anticipated from the asymmetry of the resistance table. The results from the van der Pauw
and our method are presented in Table 2.12. The resistivity given by van der Pauw
method was 12.51 (Q cm) and the resistivity from our method was 12.91 (0 cm) which
agrees within 3% error. Thus, the ﬁtted resistivity from our method agrees well with the
resistivity from the van der Pauw formula.

Another test has been done on Film 2 that had six contact pads. At this time, I
obtained two complete resistance tables for the different current directions as shown in
Table 2.13 and Table 2.14. The minus signs in Table 2.14 came from the deﬁnition of the
resistance. Table 2.13 and Table 2.14 were given by constant currents (100 nA) at
T=22.8° C. These tables show a good symmetry in each term that indicates that the
contacts are Ohmic. The geometrical factors of the ﬁlm were given by the numerical
simulation and shown in Table 2.15. All six pads were constructed on the ﬁlm as perfect
conductors, although two pads were not used in the measurement. From these geometrical
factors and the resistance tables, the conductivity of the ﬁlm could be determined by the

van der Pauw method and our method. By the van der Pauw formula,

57

Table 2.13 Resistances measured by 4-probe experiment (k9).

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

R V. V. V... V.. V. V...
I12 149.20 76.30 53.51 -70.91 -91.58 M
I ,3 76.14 162.20 66.30 83.66 ﬂ -93.36
I ,4 55.23 66.07 141.10 _10_.9_8_ 83.69 72.65
I 23 -70.40 83.03 M 156.80 81.17 -72.20
I 24 -90.63 ﬂ 83.1 1 80.85 177.80 93.12
134 M -92.50 72.08 -72.00 92.53 168.60
Table 2.14 Another resistance table for the opposite current ﬂows of Table 13. (k9)
R V12 VI3 V14 V23 V24 V34
I21 -l49.00 -76.06 -53.50 71.03 91.73 20.20
I31 -75.89 -l62.00 -66.06 -83.66 23 93.45
I 4, -55.04 -65.86 -141.00 M -83.76 -72.77
132 70.30 -83.09 M -156.70 -81.07 72.24
I 42 90.46 _9._9_2_ -83.23 -80.78 -l77.80 -93.22
I 43 29$ 92.42 -72.28 71.96 -92.64 -168.70
Table 2.15 Geometrical factors calculated from FEM (mm'l).
l" VI2 Vis V14 V23 V24 V34
I ,2 712.19 408.59 262.44 -303.61 —449.75 -146. 14
I 3 408.63 878.20 350.58 469.57 _-__5__8_._(_)_’_5_ -527.62
I 4 262.47 350.58 779.07 M 516.60 428.49
I 23 -303.57 469.67 M 773.24 391.65 -381.59
124 -449.70 £807 516.57 391.63 966.27 574.64
I 34 -146. 19 -527.69 428.44 -381.50 574.64 956.13
Table 2.16 The resistivity and contact resistances of Film 2 given by van der Pauw
formula (VDP) and our method.
p (9 cm) p ((2 cm) 7', , Ar; (k9) 72 , Arz (k9) F3 , Ar3 (kﬂ) F4 , Ar4 (kﬂ)
(VDP)
12.62 I 13.25 17.5, 0.05 27.35, 0.05 18.15, 0.05 11.65, 0.05

 

 

 

58

 

 

 

 

the resistivity p was given by 12.62 (9 cm). We used only two 4-terminal resistances
(R324 , R1243 ) in the formula because the other term R1234 was very unstable due to the

symmetric formation of the contact pads. And the same resistivity as well as the contact

resistances were obtained by our method and compared with the result from the van der
Pauw method in Table 2.16. The same 4-terminal geometrical factors (R13;1 , R1243 ) used in

the van der Pauw method were used to extract another resistivity 13.25 (Q cm) by our
method. It agreed with the van der Pauw result (12.62 (9 cm)) within 5% error. All of the
four contact resistances in the same table show good Ohmic behaviors as expected from
measured tables. Our technique works well with a two dimensional sample. It was
veriﬁed by the comparison of the result from the van der Pauw method. Therefore, we
suggest that it can be applied to the 3d sample, too. One of the advantages of our method
is that it could determine the contact resistances besides the resistivity of a sample. The
detailed ﬁtting procedure to select data set and ﬁtting parameters will be mentioned in
next section for the analysis of crystallites.

Another test simulation was done to check how the result was changed depending

on the positions of leads. When one of the leads was moved by -0.]mm in x-direction

Table 2.17 Stability test of 4-terminal resistances on the position of one contact.

 

34 24 23 14 13 12
R12 R13 R14 R23 R24 R34

 

Before (k0) -121.59 -9.007 1 12.47 1 12.46 -9.044 -121.54

 

After (k9) — 124.84 -13.06 1 11.79 1 11.79 -13.06 -124.84

 

Error (%) 2.7 44.4 0.6 0.6 44.4 2.7

 

 

 

 

 

 

 

 

 

59

(sample size is about 3x3 mmz), six terms of 4-terminal method were shifted from the

initial values as in Table 2.17. Two of these terms (R1234 , R53) given by the crossed

formations of voltage and current measurements showed the signiﬁcant changes up to
44.4% by the shifting of one lead. However, the other terms turned out stable. From this
result, it has been realized that the position of the contacts should be reconstructed

accurately for the FEM especially for the symmetric formation.

2.5.2 sample 1 (Ohmic contacts)

As an illustration of this technique on 3d samples, we used two diamond
crystallites with four terminals (Figure 2.25). A resistance matrix for Sample 1 (Table
2.18) was completed by making all 36 resistance measurements. Diagonal terms in
boldface denote 2-terminal measurements and the underlined values are 4-terminal
measurements. A 6 x 6 matrix of geometrical factors is obtained by using the FEM and
shown in Table 2.19. Diagonal terms in boldface denote 2-terminal measurements and the
underlined values are 4-terminal measurements. Within the experimental error, Table

2.18 is symmetric which indicates that r2 and r3 are Ohmic. If the contact resistances are
Ohmic, we can determine the resistivity p and the four contact resistances (I; = r, '= T; ,
i = l,..,4 ).

As a check of the validity of our method, we first calculated the resistivity from

the 4- terminal measurements (where R5” = pI‘é’" ) with a linear ﬁt through all six data

points and the origin. We then used this p to determine the contact resistances from the 3-

terrninal measurements (each of which involves only one contact resistance) and then

60

Table 2.18 Resistances obtained from the experimental measurement for Sample 1 (k9).

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

R V. V. V. V. g V. V.

I ,2 720.0 560.0 508.0 -176.0 -231.8 :55_._5_
I 13 568.7 1039.0 534.5 497.3 43$ -531.1
I l 4 518.4 540.2 986.0 M 500.0 476.9
I 23 - 174.5 498.7 M 665.0 195.1 -477.7
I 24 -230.6 :35._2 498.5 195.4 718.0 531.0
I 34 £7 -534.2 473.7 -477.6 529.2 999.0

Table 2.19 Geometrical factors from the ﬁnite element method for Sample 1(|.I.m'l ).

r V. V. V. V. V. V.
I12 1.6609 0.9635 0.7955 -0.6974 -0.8653 M72
I 13 0.9634 1.9550 0.9568 0.9916 -0.006647 -0.9982
I , 4 0.7957 0.9569 1.5608 1.161; 0.7651 0.6039
I 23 -0.6973 0.9916 Q._1_6_1_3_ 1.6889 0.8586 -0.8303
124 -0.8652 -0.006684 0.7651 0.8585 1.6303 0.7718
I 34 M -0.9982 0.6038 -0.8303 0.7717 1.6021

 

 

 

 

 

 

 

Table 2.20 Resistivity and contact resistances for Sample 1, assuming that all the contact

resistances were Ohmic.

 

 

 

 

 

 

 

 

 

 

 

p (Q cm) rl (k0) r2 (k9) r3 (k9) r4 (k9)
(a) 23.72 316.3 6.1 276.2 329.0
(b) 22.50 314.0 15.0 285.0 336.0
(c) 23.72 303.3 5.6 274.0 329.0
(d) 23.58 317.6 7.3 277.5 330.0

 

61

 

 

 

averaged the results. The results are shown in row (a) of Table 2.20. I then did a least
squares ﬁt for ﬁve parameters to the six 2-terminal measurements. The results are shown
in row (b). In (c), the resistivity was ﬁtted by 4-terminal measurement as in (a), and four
contact resistances were obtained from six diagonal terms by a least squares ﬁt. Finally,
we ﬁtted ﬁve parameters with all 36 elements of the resistance matrix, the results of
which are shown in row ((1) of Table 2.20. It is clear that all four methods give essentially
the same results. Note that one of the contacts has a very small resistance compared to the

others.

2.5.3 sample 2 (non-Ohmic contacts)

The measured resistance matrix for Sample 2 is given in Table 2.21. All
measurements were done at constant current (100 nA) at room temperature ('1‘: 20° C).
The values in parentheses in the diagonal terms were measured for a reversed current
flow. The different magnitudes in these two measurements show the non-Ohmic behavior
of the contacts. Diagonal terms in boldface denote 2-terminal measurements and the
underlined values are 4-terminal measurements. Using Table 2.5, it is clear from the
asymmetries that contact resistance 3 was very non-Ohmic, and that contact resistance 2
had a small non-Ohmic part. As discussed in the previous section, a second set of
measurements had to be taken, with the current probes reversed. This was only done for
the six diagonal (2-terminal) elements before the leads were damaged. Nevertheless, this
is a sufﬁcient set of data to determine the resistivity and all the contact resistances. The
geometrical factors are given for the Sample 2 as in Table 2.22. We processed the data for

Sample 2 in four distinct ways to test the accuracy of the ﬁtting procedure.

62

In the ﬁrst method we ﬁrst calculated the resistivity from the 4-terminal measurements
and with this p determined six of the eight contact resistances from the 3-terminal
measurements, and the remaining two contact resistances (r,' and r4 ) from the diagonal

elements with the currents reversed. The results of this analysis are shown in row (a) of
Table 2.23. In the second method we did a least squares ﬁt to obtain the resistivity and the
Ohmic part of the contact resistances (7,, i = l,..,4 ), from the six diagonal elements of

Im

a. , and the 3—terminal elements of the original

if" . From the six diagonal elements of AR

 

resistance matrix that involve both r2 and rz', i.e. the elements (R1223 , R1224 , R; , R3 ), we
were able to fit the non—Ohmic part of the four contact resistances (Ar, , i=1,..,4). The

results are shown in row (b) of Table 2.23. There is rather good agreement between the
two methods, with only a 10% difference in the resistivity. This agreement obtains in

spite of the very large and non-Ohmic contact resistance r3 that dominates the least

squares ﬁtting procedure for the diagonal elements. In general, it is probably best to
determine the resistivity from the four terminal measurements, and so as a ﬁnal check, we

used the value of p determined from the 4-terminal measurements. And a least squares ﬁt
was done to obtain the Ohmic part of the contact resistances ( 7} , i: l,..,4 ), from the six
diagonal elements of El” . The non-Ohmic parts of the contact resistance are unchanged

by this alternative method of determining the resistivity. These results are shown in row
(0) of Table 2.23. Finally, the row (d) of Table 2.23 was given by ﬁtting all parameters by

whole table.

63

Table 2.21 Resistances obtained from the experimental measurement for Sample 2 (kn).

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

R V12 V13 V14 V23 V24 V34
I12 875(875) 386 369 -432.3 -451 -_1_8,_1_
I13 392 3771(4878) 587 3240 _2_0§ -3032
I14 369.5 389 1591(1591) 21—2 1114 917

I 23 -450 3255 2_2’_1 3774(4936) 663 -3026
I 24 -471 £2 1126 654 1705(1705) 943

I 34 1041 -4207 918.9 -4190 935 5236(4098)

Table 2.22 Geometrical factors from the ﬁnite element method for Sample 2 (}.tm'l ).
r V. V. V. V. V. V.
I ,2 1.0825 0.6336 0.6049 -0.4489 -0.4777 -0.02869
I13 0.6336 1.4178 0.9440 0.7843 QQLOA -0.4737
114 0.6049 0.9440 1.4112 _Q._3§2_1_ 0.8064 0.4671
I 23 -0.4489 0.7842 9,1321 1.2332 0.7880 -0.4450
I z4 —0.4776 _0_.3;l_(_)4 0.8064 0.7880 1.2840 0.4958
I 34 @2862 -0.4738 0.4673 -0.4452 0.4960 0.9407
Table 2.23 Resistivity and contact resistances for Sample 2.
10(5) cm) 7., Ar. (k9) 7.. Ar. (k9) 7.. Ar. (k9) 7.. Ar. (k0)

(a) 62.68 -3.0 , -7.00 162.6 , 10.7 3379.6 , 600.7 675.5 , 51.3
(b) 69.08 -26.8 , 11.6 152.9 , 4.72 3361.7 , 575.9 654.4 , 7.70
(c) 62.68 19.8 , 11.6 189.5 , 4.72 3398.1 , 575.9 692.2 , 7.70
(d) 60.78 9.30 , -7.50 173.1 , 12.5 3390.5 , 600.9 687.0 , 52.0

 

 

 

 

 

 

 

In this analysis, it should be noted that there are some possibilities of error in
modeling and ﬁtting as well as experimental measurements. First, we suspect the
effective contact area in the sample is less than that apparent in the SEM photographs.
Some part of the contact might not form a good electrical contact with the sample, and
hence decreases the real contact area. The sensitivity of diagonal and 4-terminal terms has
been tested by calculating the resistance changes by shrinking the area of all the contacts
by 1/3. The 2-terminal measurements varies 30~40 % but the 4-terminal ones were
changed by less than 10%. The resistivity ﬁtted by 4-terminal elements turned out more
stable than the contact resistances included in 2-terminal terms under this change of the
contact geometry. Hence the resistivity should be ﬁtted from the 4-terminal resistance
measurements, which are contact independent resistances, and the contacts are then
obtained subsequently by ﬁtting the whole resistance table. We note that a ﬁne meshing is
necessary for reliable results for the geometrical factors for the 2-terminal measurements,
but a coarser meshing is sufficient to get reliable geometric factors for the 4-terminal
measurements. And we could get better results by doing the measurements on many
different samples that have the same resistivity and with many contacts that give more
terms for consistency checks.

Throughout our theoretical treatment, although we have taken into account the
spatial distribution of currents within the sample, there has been no attempt to consider
the charge redistribution at surfaces and interfaces, for example due to band bending.
Contact resistances, in particular, will depend on local electric ﬁelds which, in turn, will;
depend on details of the semiconductor-metal interface. In retrospect, it is somewhat

surprising that a model invoking local contact resistances should work so well, with

65

contact pad dimensions and separations approaching a space charge depth. If the scale of
the experiments were to shrink further, these considerations will need to be treated with
more rigor. Even if the reconstructed geometry is correct, one more error is possible. The
conductivity of the sample was assumed to be isotropic due to the diamond lattice
structure. So, the conductivity has been treated as a scalar, not a tensor. However, this
assumption may not be true for a sample with a twin boundary that could be observed in
Sample 2. This twin boundary may affect the isotropy of conductivity or homogeneity of
the material property. Nevertheless the reciprocal theorem is still true for anisotropic
conductivity; a resistance table is still symmetric for anisotropic tensor. We did not
analyze this effect on Sample 2 when the conductivity was determined in this chapter
since the defect is not apparent to be considered in the SEM photo.

The temperature dependence of conductivity has been one of main concerns
because it is crucial to understand the conduction mechanism of the doped
semiconductors. The temperature dependence of the resistivity in Figure 2.27 was
obtained from the temperature dependence of a 4-terrninal resistance in Figure 2.26. It

was done by rescaling the resistances with a geometrical factor for the measurement. The

Table 2.24 Fit result of conductivity model with two activation energies for data in

 

 

 

 

Figure 2.27.29
Sample a] E] 0', E2
(102 o" cm") (eV) (102 :2" cm") (eV)
Sample 1 1 14: 4 0.347 :0.001 1711.9 :345 0.109 :0.004
Sample 2 163: 5 0.351 :0.001 21.8: 17 0.059 :0.014
Polycrystalline film 0.95 :0.12 0.199: 0.004 102: 6 0.0089 :0.009

 

 

 

 

 

 

 

66

 

 

 

 

 

 

 

100000 I ' l r l v 1 1 1 I I l I ' I d
j . HF1 1
. a HF2 ‘
' 0 91:1 ' ‘ . o ‘
10000 -. ° PF2 ' .
E o MC-A 5
: o uc-a - :
d ‘ I
1000 -: 1
a . =
. 1
al 1
A
(a). 100 1 ’ 1‘ g
I 1 . .
g 1 1
o. 2 I '
cl ‘ .1
I ‘ .
1o ‘5 I, A -:
Z . r' A :
. . y' A .
ul $~ g. 4
1 .3 A ‘
47° " ‘
1 1 ‘ 1
-I a A d
: _ an a
001 I l I l V l I ' I ' I r I I I’
2 0 2.5 3 0 3.5 4.0 ‘ 4.5 5.0 5.5 6.0
-1
1000/T (K )

Figure 2.27 Temperature dependence of the resistivities of the diamond samples. Two
microcrystallites (MC-A, MC-B) correspond to Sample 1 and Sample 2.
PF 1 and PF2 are polycrystalline ﬁlms. HF] and HF2 are homoepitaxial
ﬁlms. The resistivity of these samples was obtained by rescaling the
resistance measurement in Figure 2.26 using the ﬁtted resistivity.

67

geometrical factor for the 4-terminal measurement is a constant throughout whole
temperature regime. The temperature dependence of the conductivity is governed by

a = 0, exp(— E, /kT)+ 0, exp(— E, /kT), (2.37)
where E is the impurity ionization energy and E2 is from the nearest neighbor hopping.
0'l , 0'2 , El and E2 were obtained by ﬁtting the data with this relation. I present the ﬁtted
results given by J aeger in Table 2.24.29 It was observed that the activation energy E1 was

the same for two diamond crystallites, but the polycrystalline ﬁlm presented signiﬁcantly

different activation energy. That indicates that the conduction mechanism in the

 

polycrystalline ﬁlm is quite different from the homoepitaxial ﬁlm or microcrystallites.
The discrepancy might be from the effect of grain boundaries in the ﬁlm as mentioned in
the motivation of this work. No signiﬁcant structural defects such as grain boundaries are

formed in a homoepitaxial film since it is grown on the diamond substrate.

2.6 Conclusion

A method to measure the resistivity and contact resistances of single crystals in
complex geometry has been introduced and demonstrated on micron-sized diamond
crystallites. The method is quite general, if the material of the sample is homogeneous
and Ohmic, with contacts that can be non-Ohmic but should be equipotential surfaces.
The sample resistivity and contact resistances were obtained experimentally using 4-
probe measurements, aided by ﬁnite element calculations. This technique can be applied

to determine the sample resistivity using any shape of sample. It is possible to determine

68

the sample resistivity and the Ohmic part of the contact resistances solely from the set of
2-terminal measurements, with currents in both directions, but it is probably best to use
the 4-terminal measurements to determine the resistivity. The most accurate results will
be obtained from a complete set of resistance measurements with currents in both

directions.

V_L'

This technique was tested to determine the conductivity of some materials in
complicated geometry like CVD diamond crystallites, but it should be noted that there are 3.

some restrictions on the application of this method. For a material that has an anisotropic

 

conductivity, our method cannot be applied. There are six terms in a general conductivity
tensor because of its symmetry. In our method described in this chapter, the conductivity
has been treated as a scalar, not a tensor. Therefore, we need a different approach for the
anisotropic case. The next chapter will be devoted to a presentation of a method to
determine the anisotropic tensor using the iterative linearization technique as well as the

multi probe measurement.

69

Chapter 3 ANISOTROPIC CONDUCTIVITY
MEASUREMENT FOR COMPLEX GEOMETRY

3.1 Introduction

In this chapter, I introduce a method to determine the anisotropic conductivity of
an arbitrary shaped sample. The conductivity of our samples (CVD diamond crystallites)
was assumed isotropic throughout the previous chapter due to the symmetric lattice
structure of the samples.”31 However, one of the crystallites had a clear twin boundary
between the contact pads (see Figure 3.1). The homogeneity may not be satisﬁed in this
sample because of this defect, that is a main requirement for the application of our
method. Under the assumption that the effect of the defect could be neglected, its
contribution to the transport has not been considered rigorously before.

However, this inhomogeneity of the sample induced by some signiﬁcant defects

 

Figure 3.1 A diamond crystallite has a local defect (twin boundary) denoted by an arrow.
This defect was neglected in the isotropic conductivity determination.

7O

 

should be considered more carefully in the determination of the material properties.
Moreover, some materials (Bi, Sn) and some superconductors are known to have
anisotropic conductivity tensors from the anisotropic lattice structures (see Table 3.1).
Diamond and Aluminum have an isotropic conductivity tensor from the cubic lattice
structures, but some materials in other lattice structures (tetragonal, trigonal, and
hexagonal) usually have anisotropic conductivity tensors. Therefore, the conductivity '1"
should be determined as a tensor form that has six terms in 3d. Only three diagonal terms
can be determined when the three principal axes are known.

The main difference between the isotropic and anisotropic case is related to the

 

geometrical factor. For an anisotropic conductivity tensor, the resistances can not be
represented by the simple form of a geometrical factor times conductivity. In fact, the
geometrical factor cannot be deﬁned in the anisotropic case. Nevertheless, it may be
possible to calculate the anisotropic conductivity tensor for some simple shapes (brick,
cylinder) that build up uniform current patterns. However, it is not easy to fabricate the
samples into such shapes in micron scale, and the principal axes of conductivity are

generally unknown. To identify this conductivity tensor, I developed an algorithm using

Table 3.1 Electrical resistivity of metalic crystals.32 Values of principal resistivities at
T=20° c. (in unit 10'6 :2 cm)

 

 

 

 

 

 

 

 

 

 

Crystal System pl , p2 p3
Aluminum Cubic 2.72 2.72
Tin (Sn) Tetragonal 9.9 14.3
Bismuth (Bi) Trigonal 109.0 138.0
Zinc Hexagonal 5.91 6.13

 

71

the multi-probe measurement and computer simulations (the FEM analysis and the
iterative linearization technique). This method has been tested for the 2d and 3d
computer models and a real sample (Bi) whose anisotropic tensor components are already
known as in Table 3.1.

In the next section, a brief description of the nonlinear problem induced by
anisotropic conductivity is presented and I will describe the basic idea of a method based
on iterative linearization in Section 3.2.2. The orthogonal transformation is introduced to
ﬁnd the principal axes of a tensor in Section 3.2.3. In Section 3.3, I introduce a general
scheme of computer simulation using FEM on the computer models (Section 3.3.1) and
the interfaces between the FEM package and the minimization routine (Section 3.3.2). I
will explain how to apply this method to real material (Bi) in Section 3.4. The sample
geometry and the lead formations will be described in that part. The results of the
computer models and the real sample are presented and analyzed in Section 3.5.1 ~ 3.5.3.
The stability test on our method is discussed in Section 3.5.4. To close this chapter, I will

provide the conclusion and suggestions in Section 3.6.

3.2 Theory

3.2.1 anisotropic conductivity tensor

For an isotropic resistivity tensor, a sample resistance ’R;."' has been deﬁned by

(3.1)

i] ’

3R3?" = 1".”

Where 1‘5" is a geometrical factor that depends only on the geometry and the resistivity p

18 Single parameter. The geometrical factor does not have any orientation dependence

72

since the sample is isotropic. This simple relation in Eq. (3.1) makes it possible to

determine a resistivity by ﬁtting the geometrical factors 1‘5?" and the sample resistances

5R3". However, the conductivity should be considered as a 3 x 3 tensor in general. This

anisotropic behavior of the conductance could be observed in many materials. The
conductivity tensor is symmetric when there is no external magnetic ﬁeld.33 Then, the

number of the independent terms in the tensor is given by d (d + l)/ 2 in d—dimensions. It

becomes three in 2d and six in 3d as shown in the following equations.

on a , ,
0' =[ x’ J in 2d (3.2)
on, 0'
0,, 0",). on
a: 0,, 0'”, a), in 3d (3.3)
OX2 a}! all

If the principal axes of a conductor are known and the coordinate axes in the
simulation are set up equal to these principal axes, the conductivity tensor becomes a
diagonalized tensor. Since all off-diagonal terms are set zeroes, the only unknown terms

in a tensor are the diagonal terms. Even so we cannot use the previous method to

determine the conductivity tensor because the geometrical factor 1".” cannot be deﬁned
as in Eq. (3.1) for the anisotropic tensor (5' ), i.e.
R.’."' = -R.’"’(6) ¢ p17?" . (3.4)

So, the resistivity p cannot be obtained by dividing a resistance by a geometrical factor.
All elements in a conductivity tensor are correlated in the expression of the resistances.

These elements cannot be separated from the geometrical factors as previously done for

73

the isotropic tensor. More general relations between the experimental measurements and

the theoretical calculations should be

‘R.j-'"=‘R§?"(6)+n(5. —6,-,,.)+r,'(6,,,, —6,,). (3.5)
where 'Rfjf" is an experimental measurement and r, and r,‘ represent non-Ohmic contact
resistances. We need to extract the conductivity tensor using Eq. (3.5) aided by the multi
probe measurement and the FEM analysis.

In the presence of a magnetic ﬁeld B, the tensor becomes antisymmetric and can

be separated into a symmetric and an antisymmetric part as in Eq. (3.6).

on on, 0 on on “also,2 0
o = 0",, a” 0 = on, -RZB0'02 0),, O
0 0 oz, 0 0 oz,
2 , (3-6)
on on, O 0 R200 0
= on. 0,, 0 +3 —R,o,2 0 0
0 0 on 0 0 0
where RZ is the Hall coefﬁcient and 0'0 is given by
/2
0'0 = (03.0"... —ox}.2) . (3.7)

However, this Hall measurement will not be considered in this work.

3.2.2 iterative linearization technique for nonlinear inverse problem

The unknown conductivity tensor can be determined by ﬁtting two resistances
(experimental and theoretical) in Eq. (3.5). Since the functional forms of a resistance in
terms of the tensor elements are generally unknown, we should use the iterative

linearization method. While the tensor elements are changed, the linearized resistances

74

calculated from the FEM are fitted to the experimental measurements to extract the
optimized tensor elements. Usually the optimized solution can be achieved after many
iterations of the linearization and least square fittings. This technique is also known as the
Gauss-Newton method.34 The term Gauss is named by the least-square ﬁtting of the
resistances and the term Newton comes from the iterative linearization technique in the
optimization procedure. This technique has been widely used to solve nonlinear
Optimization problems, for instance the earth conductivity problem in geophysics.35

Here I describe the detailed procedures of this technique. Suppose we are trying to
ﬁnd an anisotropic conductivity tensor with unknown principal axes. Simply we choose
an arbitrary tensor to start the iteration, which is deﬁned by six parameters. To get the
linearized form of each resistance for this tensor, all ﬁrst derivatives of the resistances
with respect to the six tensor elements are calculated using the FEM. Then, the
expressions of the resistances are known as linearized forms, we can perform the least
square ﬁtting to extract a tensor. However, since the linearized forms are correct only
near the previous tensor, the ﬁtted tensor may not be the global solution of the
optimization. Therefore, the ﬁtted tensor substitutes for the old tensor and the same
procedures are done. This iteration will go on until the best-optimized tensor is obtained.
After the tensor is obtained, it is diagonalized by the orthogonal transformation to ﬁnd the
principal axes and values of the tensor.

In this iteration, the ﬁrst step is very important. The iteration begins from an
arbitrary selected point that might be far from the global solution. Then it is very crucial
to ﬁnd a next point close to the global solution in the ﬁrst step. We begin the ﬁrst step

with an isotropic tensor (unit conductivity) and keep the tensor isotropic throughout the

75

 

 

 

 

 

 

 

 

p‘ p'+Ap p" p2 p°p°+Ap

Figure 3.2 Iterative linearization technique of 12 minimization in 1d (Newton method).

This technique is useful when the function form of the 12 is not known.

Each step, the function form is approximated to a linearized form by Taylor
expansion. After many iterations, it reaches the minimum point.

76

step. This procedure is exactly the same as our method in the previous chapter. We can
determine the scale factor po that is nothing but an isotropic conductivity. When a tensor
is isotropic, there is only one parameter p0 ( = p” = p),y = pa ). Since the conductivity is
isotropic in this step, we can use the relation of

R = pOFO , (3.8)
where F0 is a geometrical factor given by the initial tensor that is a unit conductivity

tensor. The ﬁtted tensor that will be used for the next step has same diagonal terms of p0

as

 

p0 O O
p°= O p° 0 . (3.9)
O 0 p°

From the second step, we use the linear approximation method of the resistances.
Each tensor element is treated separately in the procedure. A resistance depending on the

six ﬁtting parameters can be represented by

RU)”,p).).,pzz,pn,pyzipxz)=R0032,pg)»pgrpgnpgz’pgz)+2Aleprj ’ (3°10)

i< j
where R0 is a resistance determined by the resistivity tensor p0 given by a previous step.

The linear expansion of the resistances can be deﬁned by

(3.11)

 

R 9+A .. —R .9
pij=p3+Apij , Al}: A5— 2 (p0 p") (p0).
Apij Apr}

Selecting Apij is somewhat arbitrary depending on the functional shape around the tensor

p3. The constant Ar; can be given by taking the numerical derivative of the resistance.

77

Then the resistances in Eq. (3.10) can be ﬁtted with the experimental values to ﬁnd a
conductivity tensor that minimizes the 12 of the resistances. By continuing this process
until the best optimum point is achieved, the anisotropic conductivity tensor and contacts
can be given. This is the main idea of this algorithm.

Next, I explain the way to save one ﬁtting procedure in the ﬁrst step. Each
iteration in this simulation requires a heavy computer job; sometimes it is very efﬁcient to
skip one iteration. Previously, I consider only an isotropic tensor in the ﬁrst step for the
rescaling of the resistivity. However, an anisotropic tensor also can be used for the
rescaling and the determination of each tensor term simultaneously in the ﬁrst step. First,
I will describe the skip for a diagonalized tensor. It will be explained for a general tensor
later. Suppose we have a diagonalized tensor. Then, we have only three unknown

variables in a tensor. A resistance can be represented by a geometrical factor F0 and some

constants Bi}. as in Eq. (3.12).

 

u + W + :2
R(pxx’pyy’pzz)=r0[p p3 p

+Bupu +Bwp”, +Bapa] . (3.12)

Since a resistance R should be equal to F0 for the unit tensor ( p“ = p”, = pZZ = 1.0 ), a
relation between three coefﬁcient B”. can be found by

BU+BW+Ba=O. (3.13)
Then Eq. (3.12) can be simpliﬁed as

R(p...p.,..p.) = 131(1)) +3....(p... -p..)+B.(p. -p..)l . (3-14)

where

78

 

p.+p.+p.
<P>= " -

3 (3.15)

One coefficient Bu has been removed from Eq. (3. 12) using Eq. (3.13). If we deﬁne new

parameters of 1,]. and Ba by

py‘pu AR 1
x.~=——— ,B,.=r —-, (3.16)
J (p) } ()[Apu'] 3

then we have a relation between a resistance and the ﬁtting parameters (p) and Zr;-

 

R(pxx’ pyy’ pzz) = I-‘0<I))[1-'- B)’)‘XX)’ + Buxxz] (3’17)
In the first step, the rescaling is done by (p) while 1,}. provides the separation of each

tensor element. If the off-diagonal terms are considered as well, this relation can be

generalized by adding some extra terms 6,]. associated with off-diagonal terms in Eq.
(3.12).
R(pxx’py_y’pzz’pxy’pyz’pxz) : r0<p)[1+Byylxy +821sz +Bxy§xy +Byz€yz +sz§xz]’ (3'18)

where new parameters of 5,]. can be deﬁned by
15,, = 3L. (3.19)

So we have six ﬁtting parameters (< p ), x”, x”, 5.. , 5 Cu) in the resistances with

yz ’

six geometrical factors (I‘0 , 8 Ba , B”, B sz ) calculated at the starting point.

However, this amended ﬁrst step has been used only for especially heavy computer jobs
in the real simulation. Usually I applied the simple algorithm for the ﬁrst step given by

Eq. (3.8).

79

In each iterative step, we minimize 12 of the experimental resistances R‘ and the

theoretical resistances R ' ( [5, rj ) . The theoretical resistances can be given by a resistivity

tensor [5 and the contact resistances rj as in Eq. (3.5).

2

R? _ R-t ”, r-

. .(p ,) , (3.20)
0 i

l
where i denotes the index of all data and 0',- is the standard deviation of each experimental

measurement. We used Levenberg-Marquardt method to minimize 12 to determine the

unknown parameters (srx terms 1n a resrstlv1ty tensor and the contact resrstances).2

3.2.3 orthogonal transformation

After the iteration, an optimized conductivity tensor is obtained. However each
element in this tensor varies under the rotation of the coordinate; it is changed for
different coordinate systems. So, the principal axes and diagonal elements should be
obtained. The sample properties are characterized by these invariant quantities. A
symmetric conductivity tensor in Eq. (3.3) is a Hermitian matrix that has three real
eigenvalues and three orthogonal eigenvectors.36 Three eigenvalues (principal values) and

eigenvectors (principal directions) of the matrices can be determined by doing the

Orthogonal transformation on the tensor.

P
011 012 013 0.11 O O
013 023 0'33 0 0 0.3,;

80

 

To identify the principal axes and values, we solve an eigenvalue problem of the tensor

0,}. .

6 - i. = 2,55,. , (3.22)

l

where A, is an eigenvalue that denotes a conductivity and If, is an eigenvector for the

eigenvalue.

3.3 Computer Simulation

3.3.1 sample geometry
First, I have counted how many leads are needed on the sample for enough

information to determine the tensor. Most ﬁttings in this chapter have been done with
only 4-terminal measurements to obtain the tensor. The contact resistances were not
considered in most cases. In the previous chapter, the ﬁtting with 4-terminal terms turned
out to give a better result than the ﬁtting with whole terms including 2-point, 3—point
terms. The isotropic tensor can be determined by four leads since the tensor has single
unknown parameter that is fewer than the independent 4-terminal terms. However, the
number of the unknowns is increased for an anisotropic tensor, so we may need more
leads. Moreover, it is hard to form many contacts on some samples. Therefore, it is
important to know what is the minimum number of the contacts necessary to reconstruct
the conductivity tensor. Since there are six unknown parameters in the tensor now, we
Should have more than six independent 4—point resistances. Following the argument in

Appendix A.2, we should have at least six leads on the sample. To identify the contact

81

 

 

(b)

Figure 3.3 Many probes are attached on the surface of samples in 2d (a) and 3d (b). More
contacts are needed for anisotropic conductivity tensor compared to isotropic
conductivity tensor.

 

resistances also, we need six leads. Now there are 6 + n unknowns (6 tensor terms and n
contact resistances) and this number should be less than "C, (= n(n — 1)/ 2) that is the

number of the independent terms in a resistance table. Therefore, we need six leads in
both cases whether the contact resistances are ﬁtted or not.

To illustrate this technique, we generated computer models in 2d and 3d as shown
in Figure 3.3. Six surface probes were formed on the sample both in 2d and 3d. By the
counting argument of the 4-terminal terms, only ﬁve probes were enough to determine
three tensor elements in 2d. Six circular contacts of different sizes were formed in the 2d
planar sample and six identical small bricks of contact pads were attached on the surface
of the large brick of 3d sample. The contacts were randomly distributed on the 2d sample,
but for the 3d sample, the six contacts were attached one on each facet. A 2d sample was
meshed by ~3000 brick elements (DC2D4 in ABAQUS) and a 3d sample was by ~5000
tetrahedral elements (DC3D10). Modeling (generating geometry and meshes) was done

again by the commercial package AN SYS.24

3.3.2 simulated experimental data using FEM

To test my method for these computer models, I generated simulated resistance
tables. An anisotropic conductivity tensor has been arbitrarily selected to generate the
simulated resistance table. By performing the FEM with this prescribed tensor, I obtained
a 15 x 15 resistance table constructed by all possible voltage and current measurements.
With this tensor, the contact resistances were added to the sample resistances. Only
Ohmic contacts were considered for this test sample. Moreover, to simulate the real

experimental measurement, the resistance table was decorated by Gaussian noise. This

83

simulated resistance table was used as experimental resistances to test how our method
could ﬁnd the prescribed tensor and contact resistances. A commercial package

ABAQUS25 was used for the analysis.

3.3.3 iterative linearization fitting using FEM

In each iteration, the resistance calculation using FEM was done by the ﬁtted
conductivity tensor obtained from the previous step. And these resistances were ﬁtted

with the experimental values again to determine another tensor for the next step. In each

step, I calculated the whole resistances R5" (p j) (225 terms) using the fitted conductivity
tensor p j given from the previous step. And some other resistances R5?" ( p j. + Apj) were

calculated to obtain the linear expansions of Ag. in Eq. (3.11) or By. in Eq. (3.16).

 

 

 

 

AR?!"
Ar = '1 , (3.23)
Apj
ARI?" 1
3:1 = 1‘0 -- , (3.24)
Apj 3
ARI?" _ Ry“); +Apj)-R,;m(pj)
Apr A” 1' (3.25)

The least square ﬁtting could be done in each step using numerical derivatives
given by the FEM. The whole table or only 4-terminal terms were used for the ﬁtting.
After approaching close to an optimum point, the iteration was stopped. The ﬁtted
conductivity tensor and contact resistances were compared with the initially given

quantities. While the ﬁtting was done iteratively, I obtained a diagonalized conductivity

84

tensor by the orthogonal transformation. The diagonalized tensor showed the principal

axes of the tensor and the principal values along the principal axes.

3.4 Experiment

3.4.1 sample preparation (Bi cylinder)

A test of this method has been performed on a bismuth cylinder. As shown in
Table 3.1, bismuth has an anisotropic resistivity tensor. It has two distinguishable

principal directions: one (p p = p1 = p2) is the isotropic on xy-plane and the other ( p,)
is perpendicular to the plane. Later, I will call the p p a planar resistivity and p3 a

perpendicular resistivity. This property comes from its trigonal lattice structure. The
planar resistivity component is known to be 109.0 (1045 Q-cm) and the perpendicular
component is 138.0 (10'6 Q-cm). The discrepancy between these two elements is about
40%; this number is larger than the possible error estimations of the measurement or the
simulation. In Figure 3.4, I display the geometry of our test sample. It is a relatively thin
cylinder that may be treated as a circular plate. Its diameter is 10.4mm and thickness is
2.0mm. Six gold electrodes were formed on both circular planes of the sample by the
evaporation technique. All the contact areas were circular (diameter : 1.0mm). Three
contacts were formed on the top, the others were formed on the bottom. Three pads on a
plane were rotated by some angle from the other three pads (see Figure 3.4(b)).
Otherwise, we would have measured extremely small 4-point resistances that might have

caused some serious numerical error. Usually the symmetric formations of contacts raise

85

 

(a)

top contact

bottom contact

 

(b)

Figure 3.4 Schematics of a real sample (Bi cylinder) with six circular Au electrodes. (a)
It is a relatively thin cylinder. Three contact pads were formed on each
circular plane of the sample. (b) Three pads on a plane were formed to have
an angle from three pads on the other plane.

86

«tartan... .1.

,.
“a" ems.
. . ‘3.

 

(b)

Figure 3.5 The normal photos of the Bi sample for the real measurement.37 It has six
contacts that are formed on both sides of the circular planes. (a) and (b) show
the formation of contacts on these two sides. The ruler shows the size of the
sample.

87

the errors in the simulation because of the sensitivity of the 4-terminal resistances that has
been observed in the isotropic case. The geometry was constructed by taking normal

photos of both sides of the sample as given in Figure 3.5.

3.4.2 six-probe measurement

A resistance table (15 x 15) was obtained by all possible voltage and current
measurement. A total of 225 terms were obtained in a resistance table. Only 90 terms of
4-point measurement would be used in the ﬁtting since the contact resistances were not
concerned for this real material. DC current has been used for the I-V measurement on the
sample. It was important to use low currents to remove the thermoelectric effect of the
high conductivity sample. Each 4—point resistance was measured to check the Ohmic
property of the sample. This was really important to check whether the conductance was
related to the thermoelectric effects. Thermoelectric effects break the symmetry in the 4-
point measurement. So, if these effects are signiﬁcant in the measurement, their
contributions to the resistances should be considered separately. With several different
DC current flows, the I-V curves for all 4-terminal measurements were obtained and ﬁtted
by the linear regression technique to ﬁnd the derivatives that were equivalent to the 4-
terminal resistances. Among 90 terms, some measurements showed bad I-V
characteristics. Some 4—point measurements were observed to be asymmetric. These
terms were removed from the ﬁnal data since I had enough data in the table. The
asymmetry might be normal for 2-point and 3-point resistances because of the non-Ohmic
contacts. However, 4-point resistances should not be asymmetric if the conductance is

only electrical. That asymmetry might be from the thermoelectric conductance. However,

88

that effect has been neglected throughout this work based upon the asymmetry was not
substantial. So, the 4—point resistances were symmetrized by taking the average of two
off-diagonal elements. Finally, this reﬁned resistance table was used to determine the

conductivity tensor elements of the Bi sample.

3.5 Analysis

3.5.1 computer generated sample (2d)

Our method has been tested on a 2d computer model in Figure 3.6. That irregular
plane has six circular probes. Figure 3.6(a) shows the meshed structure and Figure 3.6(b)
presents a voltage distribution on the sample for a measurement given by the FEM
analysis. Table 3.2 shows the prescribed and the ﬁtted conductivity tensors. The ﬁtting
procedures were done by whole resistance table and 4-terminal resistances only. The

second row represents the conductivity elements that have been used to generate the

   

(a) (b)

Figure 3.6 A computer generated sample for the determination of anisotropic
conductivity. This 2d computer model has six circular probes. The contour
map (b) shows the equipotential lines. The meshing was done by ANSYS and
the analysis was done by ABAQUS.

89

4-point terms only in 2d (error-'0.02)

  

 

 

1 _ :
i 7 5 'E -I— 0,,
5.: .. :MIWI-w- ............. :12+6
I 1 . - §~A—o:
4.‘ ..................................................................... : ...........
1 .
,= 3:.
b ; 2... A ~—-—a———a-——A
2L.......,,x ......................... ...............................................
4 l/ -
I". . . . .

 

whole terms in 2d (error=0. 02)

5'.
: f S f S i E-I-o.
5.: . l/._-__. ............ _._ 0,;
‘ =~A~622

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

 

 

. 0
q. . .. . .- ... ..... .. ...._..—........... ..o

 

o—O—O—O
. l i 6 i
step
(b)

Figure 3.7 These ﬁgures show the ﬁtting procedures of the tensor in 2d. The ﬁttings
were performed using only 4-point resistances (a) or whole resistances in the
table (b). Three data points are equivalent to three conductivity tensor terms.

90

 

 

Table 3.2 The comparison of the ﬁtted results of the conductivity tensor elements with
the prescribed values for the computer model in 2d.

 

 

 

 

 

 

 

 

 

0'll 0',2 0'22
Prescribed 5.00 1.40 2.50
4-point terms (a) 5.03 1.42 2.52
Whole terms (b) 5.01 1.43 2.52

 

 

Figure 3.8 Two principal directions were drawn on the surface of 2d computer model.
Three directions (one real and two ﬁtted) represent a remarkable coincidence.

91

simulated resistances. We compared the ﬁtted results in the table with the predeﬁned
tensor. In both cases, our method reproduced the real values remarkably within acceptable
error levels that were associated with the Gaussian noise. For this 2d sample, 2% of noise
had been added to the resistances to simulate the experimental measurement. Figure 3.7
shows the ﬁtting procedures for two data sets (4-point terms (a) and whole terms (b)).
Usually the iterations were observed to be stabilized just after a few steps.

The orthogonal transformations were done to diagonalize the ﬁtted tensors. Two
principal values and two principal directions equivalent to the values were given for the
real tensor and two ﬁtted tensors. The diagonalized tensors and the principal vectors of
each tensor were presented in Eq. (3.26)~(3.28). Eq. (3.26) came from the prescribed
tensor and the rest came from the ﬁtted tensor. These principal directions in each case

were shown on the surface of the model as in Figure 3.8.

5.627 0.000 56, =(0.913 0.409)

00: y _, . (3.26)
0.000 1.873 x,=(—0.409 0.913)
5.67 0.00 “ = 0.912 0.411

,= . f‘ ( ) . (3.27)
0.00 1.88 x,=(—0.411 0.912)
5.66 0.00 ‘ = 0.910 0.414

6,: , f‘ ( ) . (3.28)
0.00 1.87 x,=(-0.414 0.910)

3.5.2 computer generated sample (3d)

I did another test simulation on a 3d computer model with six probes. A meshed
structure of the sample and the voltage distribution on the surface were presented in

Figure 3.9(b). The predeﬁned conductivity tensor and the ﬁtted results were given in

92

Table 3.3. Since we have six terms in a tensor, the number of the ﬁtting parameters is also
six. Two different ﬁttings with different ﬁtting data (4-point terms and whole table) were
done again for this 3d sample. These results by two different data sets (4-point (a) and
whole table (b) in Figure 3.10) showed a good agreement with the predeﬁned tensor
elements. The errors in the results might be acceptable, although they came out higher
than in 2d case. 2% Gaussian noise had been added to generate the simulated resistances.
Usually a few more steps were needed to reach to the Optimum solution compared to the
2d case. It is simply because of the increase of the ﬁtting parameters in 3d. From these
graphs, I realized that the ﬁrst a few steps were critical in the entire process. The rest of
the iterations did not improve the results in most cases. I obtained the principal values and
axes of the conductivity tensor by diagonalizing the ﬁtted tensor. Eq. (3.29) shows the
diagonalized tensor of the prescribed conductivity tensor and the principal vectors. And
Eq. (3.30) and (3.31) were given from the ﬁtted tensors by 4-point terms and whole

resistances respectively.

Table 3.3 The comparison of the ﬁtted results of the conductivity tensor elements with
the prescribed values for the computer model in 3d.

 

O.ll 0'12 022 0'13 0'23 033
Prescribed 8.0 -0.8 1 1.0 1.5 -0.5 9.0

4-point terms (a) 7.640 -0.817 10.19 1.521 -O.613 9.604
Whole terms (b) 7.845 -O.811 10.24 1.619 -0.552 9.366

 

 

 

 

 

 

 

 

 

 

 

 

93

 

(b)

Figure 3.9 The meshing (a) and equipotential lines (b) of a 3d computer generated
sample. Six brick-shaped probes were formed on six facets of the sample.

94

4-point terms only In 3d (error=0. 02)

 

 

 

whole terms in 3d (error=0.02)

 

(b)

Figure 3.10 These ﬁgures show the ﬁtting procedures of the tensor in 3d. The ﬁttings
were performed using only 4-point resistances (a) or whole resistances in
the table (b). Six parameters in a tensor were obtained each step. The
iterations in 3d needed a few more steps to the optimum than in 2d.

95

qﬁ '.--lif'0tl n"'”

6.89 0.00 0.00 2, =(0.822 0.0916 —0.562)
a: 0.00 11.55 0.00 . 32,=(—0.352 0.857 —0.376) (3.29)
0.00 0.00 9.57 2, =(0.447 0.507 0.737)

6.76 0.00 0.00 Jr',=(0.886 0.131 —0.445)
6,: 0.0 11.27 0.00 . 2,=(—0.411 0.666 -0.622). (3.30)
0.00 0.00 9.41 35, =(0.215 0.734 0.644)

6.77 0.00 0.00 x, =(0.854 0.119 —0.507)
6,: 0.00 11.24 0.00 . 224—0441 0.682 —0.583)- (3.31)
0.00 0.00 9.44 2, =(0.277 0.722 0.635)

The sensitivity test of our method has been performed for several simulated
resistance tables that were generated by different noise levels. I monitored how the errors
in the ﬁtted results were related to the noise levels and shown in Figure 3.11. It shows
two curves that correspond to the ﬁtting data of whole table and 4-point terms. They
turned out basically linear as expected, but the derivatives were quite different. The
derivative of the ﬁtting with whole table was somewhat higher than the ﬁtting with 4-
point terms only. This result seems to support the argument that the ﬁtting with 4-point is

more efﬁcient than the ﬁtting with whole table in the determination of the anisotropic

tensor.

3.5.3 real sample (Bi cylinder)

The test for a real material has been done by two different meshings for the same
geometry. Figure 3.12 presents one meshed sample and the voltage distributions of the Bi
cylinder obtained by the FEM analysis. Two ﬁgures obtained from different view angles
show the voltage distribution on the top and side of the sample. Notice that whole

contacts were given as equipotential regions since the six small cylindrical contact pads

96

Fitting Error (%)

Sensitivity Test for 3d Computer Model

40

   

-+- all terms
-- El" 4-point terms;
30 . ‘ .

 

 

25

20

15

+
10 / \
5 _I + . . . ...............
~ ........... 1:1» '3
. ...... Do.
Om'...‘ﬁfﬁ...,......
0 5 10 15

Noise level (%)

Figure 3.11 Sensitivity test of the conductivities ﬁtted by 4-point terms and whole table

has been done to check the relation between the result and the Gaussian
noise.

97

 

Figure 3.12 Equipotential lines of the real sample (Bi) obtained by FEM analysis.

98

were perfect conductors. The well-known conductivity elements of Bi and the ﬁtted
results by our technique are presented in Table 3.4. And the ﬁtting procedures of two
different meshings were shown in Figure 3.13(a) and (b). The real conductivity
components (planar (109.0) and perpendicular (138.0)) of the Bi sample were represented
by the two horizontal solid lines in the ﬁgures. In this test, only 4-terminal resistances
were used in the ﬁtting according to the argument in the previous tests. It has been
realized that the ﬁtting with 4-terminal terms was more effective than the ﬁtting with
whole table. Since the orthogonal transformation has been done in each step, I have only
three data points in the ﬁgures. Three points correspond to the three principal components
of the conductivity tensor. The errors of the ﬁtted conductivities were measured relatively
high compared to the results of the computer models. This result might be acceptable
when the experimental and simulation error levels are considered. However, the
separation between the planar and perpendicular components is not so clear in the results.
IfI take the average of two planar resistivity values (98.75, 126.3), it is given by 112.52
that agrees remarkably with the real value (109.0). However, these two components never
crossed over in the entire ﬁtting procedures. The separation between these two planar

components was quite substantial. Moreover, the ﬂuctuations of the ﬁtted components

Table 3.4 Principal values of the resistivity of Bi. (1045 0 cm)

 

 

 

 

 

 

 

 

 

P11 P22 P33
Real Sample (Bi) 109.0 109.0 138.0
Fitted result 98.75 126.3 152.1
Error (%) 9.4 15.9 10.2

 

 

99

 

 

 

 

 

 

 

 

 

.A-AA
‘ A 3‘ “.I “A"
150- , .
. _A 'Ar-AA AA .5
A ‘ i". , ,0". I... .. ....
E 1 '0' "o "0“ 'o-O' '
"9:: 100 4 E‘IIIII'II I ...[I-I'I'I'I'I'I
Q .. . ....... pp"
50- 5 p
...... 922
.1...» ppaa
0*-T , .......,....T.fj-,
o 5 1o 15 20 25
step
(a)
250-
: A ----- A .
2001 A "A ----- A M‘s-A A
.. - *4
E -1
9
G
to :
2
of :‘
0a I I I l

 

-
q
-
d
a
d
d

(b)

Figure 3.13 Fitting procedures of a resistivity tensor of the real sample. Two meshings
were used to test our method on this sample. Only three principal
components of the tensor were shown in these pictures. Two planar
resistivity components (solid circle and square) are observed stable around

the real value, but the perpendicular component (triangle) shows a
substantial variation, especially in (b).

100

were observed to be signiﬁcant in both ﬁtting procedures as in Figure 3.13. By the
comparison of two ﬁtting procedures, the fluctuation of the perpendicular component was
more serious than the other two planar components.

The experimental errors were believed not to exceed 10%. So, there must be some
other considerable error sources in this method. One of the main error sources might be
the incorrect geometry reconstruction of the sample. Therefore, in next section, I will
describe various stability tests to investigate the possible error sources, the meshings and

the geometry reconstruction.

3.5.4 stability test

First, I selected several different sizes and types of elements for the meshing and
checked if the ﬁtted results depend on the way of meshings. The characteristics of all
meshings used in this test were presented in Table 3.5. Iran the algorithm with these
meshings to obtain the anisotropic conductivity tensors. The ﬁtting procedures were
given in Figure 3.14. The planar components (solid square and circle) fluctuate slightly
around the real value (109.0) that is denoted by the lower solid line for all meshings, but
the variation of the perpendicular component is quite signiﬁcant. Therefore, I realized
that a ﬁner meshing did not improve the ﬁtting result so much. The variation of the
perpendicular component was quite substantial in the simulation.

The other test has been done to ﬁnd the dependence of the result on the positions
of the contact pads. I built up several structures that had contact pads on slightly different
positions by rotating three pads on a plane. The ﬁtted results by these different structures

were presented in Figure 3.15. It was observed that even small rotations (2° ~ 5°) could

101

give substantial differences to the results. Interestingly, the two planar components were

found to cross over at some degree between 2° and 5° . And the perpendicular component

was observed fluctuating signiﬁcantly. It should be emphasized that the geometry

reconstruction is very important especially in anisotropic case.

Meanwhile, I tried to understand why the perpendicular component p3 was not

accurate and showed a serious instability. I thought that the error might be related to the

thickness of the sample. To verify this, another test was done for two identical computer

models except for the thickness. One model (Cylinder 1) was the same shape of the Bi

sample and the other one (Cylinder 2 in Figure 3.16) had different thickness that was

5mm (note that the thickness of Cylinder 1 was 2.0mm). Both samples had six probes at

same positions. I generated the simulated resistances for each model as done before. A

meshing (Bi 10-1) was used to generate the resistances of Cylinder l and a different

meshing (BulklO-l) was used for Cylinder 2. I further tested if the technique could

reproduce the prescribed conductivity components. I used the four meshings described in

Table 3.6 including the meshing that was used to generate the simulated resistances. In

Figure 3.17, the perpendicular component for the thin sample (Cylinder 1) could not be

Table 3.5 Different meshings for the stability test of a cylindrical sample.

 

 

 

 

 

 

Meshing Bi10-1 Bi10-2 Bi4-l Bi4—2
Element type DC3D10 DC3D10 DC3D4 DC3D4
Nodes per elements 10 10 4 4
Total Nodes 4473 1744 2445 728
Number of elements 2492 908 9809 2492

 

 

 

 

 

 

102

 

 

 

 

250- :
1 A a --------p
. ............. j .........
200- A ------ § 5 x p’
: é .....‘ ........ pz
5: 150-:
.92 , ¢ ........................... . .......................... g
E100: - ........................... I ........................... I
4 5
sol
d
0‘ i i i
Bi10-1 Bi4-1 Bi4—2

Figure 3.14 Resistivity ﬁttings were done for the Bi sample using three different
meshings. BilO—l is a ﬁne meshing and Bi4-2 is a rough meshing. It turned
out that the perpendicular resistivity component was so sensitive on the
meshing while the other components were stable.

 

 

 

 

300-
‘ u.- I... px
: .....H p ‘
a V .......
250 . M‘ P, A ...............
I . . IA”
5 zoo-I ‘A'
G : .
«9 .
$3 150- 4Q..- i
‘I; 2.7 """""" . ....... . ....... ..
Q ..l... '- . ....... .4]
100-I .......... l IIIIIII . IIIII ,,,,,
...,
50;
d I: .
I A
0 r I l l I r I 1* I
-3 -2 -1 0 1 2 3 4 5 6
Angle (degree)

Figure 3.15 Resistivity ﬁttings were performed on some modiﬁed structures. The
modiﬁcations were done by rotating the contact pads. The dependence of the
result on the angles was shown in the ﬁgure. Huge variations were observed
for all components. Interestingly, two planar resistivity components ( p1 ,

p v ) crossed over somewhere between 2° and 5° .

103

1'»

4:1

'9

_§_..

~

‘1‘

V

2.4-‘3

 

Figure 3.16 Equipotential lines of thicker cylindrical sample (Cylinder 2). This model

has more layers of elements along the perpendicular direction. That may
enhance the ﬁtting result for that direction.

200 -

 

 

 

 

 

4
q "...-....px
II E I E E "..."pr
150- ‘ g g ""‘m'P,z
E s
9 E
C, ‘ ...... . ................... 1 ................. ...:
.90 100‘ I -. ................... I IIIIIIIIIIIIIII . . 3
1: ~ 2 4., s ""l
of . é ' g .. =
so- ‘
o I I I I
Bil 0-1 Bi10-2 Bi4-1 Bi4-2

Figure 3.17 This graph was obtained from another computer test for the same geometry
of the Bi cylinder (Cylinder 1). The simulated resistances were generated
by the meshing BilO—l. The ﬁtted results from different meshings were

presented. The perpendicular component (triangle) shows still a huge
variation for this computer model.

 

 

 

 

 

200‘ : ; 1 i
- 'Px
. 1‘”
150- ? § 3 up:
A ‘1 ‘ “v -------- i :
. : ..... ‘ ...................
g : a s 7‘
c, : : ............... .
6 G vvvvvvvv W111 ------- a
o1oo- s 5 """""" -
cf
50-
0 r r i i
BulkO Bulk1 Bulk2 Bulk3

Figure 3.18 Another computer test for a thick sample (Cylinder 2) as shown in Figure
3.16. The simulated resistances were generated by the meshing BulklO—l.
The ﬁtting from other meshings could reproduce the all resistivity
components nicely including the perpendicular component.

105

Table 3.6 Different meshings for the stability test of Cylinder2.

 

 

 

 

 

 

Meshing BulklO-l Bulk10-2 Bulk4-1 Bulk4-2
Element type DC3D10 DC3D10 DC3D4 DC3D4
Nodes per elements 10 10 4 4
Total Nodes 9359 4473 2445 728
Number of elements 5727 2492 9809 2492

 

 

 

 

 

reproduced in some different meshings. For a rough meshing (Bi4—2), a noticeable error
was detected for the planar components. However, for the thicker sample (Cylinder 2), all
different meshings could reproduce remarkably the prescribed components as shown in
Figure 3.18. From this test, I realized that the real Bi cylinder was too thin to have enough
information of current path along the perpendicular direction. Most current flows were
formed in the horizontal direction. There was no signiﬁcant current flow into the
perpendicular direction, that is required to determine the conductivity component
corresponding to the direction. So, I suggest that the sample should be a cube-like shape
and the leads should be arranged to detect current paths in all principal directions;

otherwise the method may fail to determine all components of an anisotropic tensor.

3.6 Conclusion

In this chapter, I introduced a technique to determine the anisotropic conductivity
tensor in complex geometry, which has been generalized from the method for the
determination of isotropic conductivity. This method has been tested on computer models

in 2d and 3d with arbitrary shape and real material (Bi). Six probes were required to

106

 

 

identify the six unknown parameters in a symmetric conductivity tensor. The iterative
linearization technique has been adapted to deal with the problem where the resistance in
terms of the tensor components were unknown. This method has worked remarkably in
the determination of the prescribed conductivity tensor for the computer model. However,
for the real sample (Bi cylinder), the technique did not work properly. Some signiﬁcant
errors, especially for the perpendicular component were observed in the resulting
conductivity tensors. After some stability tests performed by various computer models, it
was realized that the resulting errors were inevitable for such a thin sample. For a test on
a computer model that was more cubic, the method turned out to determine the
conductivity much better. This technique has been demonstrated for the conductivity
measurement, however it can be applied more generally to evaluate other properties. I
believe that this algorithm could be used to determine any material property in a complex

geometry.

107

 

Chapter 4 EFFECTIVE CONDUCTIVITY OF A
MEDIUM WITH CIRCULAR INCLUSIONS

4.1 Introduction

In the previous two chapters, I described how to determine a material property
(conductivity) using experimental measurement and theoretical calculation. The samples If“
were actually single materials with isotropic or anisotropic conductivities. Only one
material (diamond or Bi) was considered in the conductivity determination. The

numerical techniques using FEM analysis were introduced and demonstrated to determine

 

the unknown conductivity of samples with a complex geometry. However, in this chapter,
my concern is with the (effective) material properties of a composite material. The study
of material properties of composite materials has been the focus of much research in
condensed matter physics, material science, mechanical engineering, and so on. Due to
the interests in diverse ﬁelds, the study of the effective medium theory on a composite
material has been conducted in various ways. Meanwhile, some efforts have been devoted
to the investigation of the universal relation between various material properties (electric
conductivity, heat conductivity, magnetic permeability, viscosity, even elastic constants)
of the materials. Analytic approaches have been successful in some simple geometries.
However, these approaches are usually limited when some complex structures are
concerned. Therefore, the numerical analysis should be adapted to study the effective
medium theory in more complex geometries. The numerical simulation using FEM could

be introduced for this purpose. To check the performance of this method for composite

108

materials, it should be tested on a composite in simple geometry where analytic solutions
can be obtained.

In next section, I give a brief background of our study on composite materials.
That part describes how the research has been done historically. The motivation of this
work is also mentioned. In Section 4.3, the sample geometry is deﬁned such as a
composite that includes two hyperspherical inclusions. Two resistances from 4-point and
2-point measurements are deﬁned for the geometry. The analytic solutions and
asymptotic forms of the resistances are described, too. The computer simulations using
FEM method are introduced and performed to calculate the resistances numerically in
Section 4.4. The numerical results are presented and compared with the analytic solutions
in Section 4.5.1(2d) and 4.5.2(3d). Furthermore, the relations between the 4-point and 2-
point measurement in the geometry are obtained and discussed in Section 4.5 .3. I will

close this chapter by presenting the conclusion and suggestions in Section 4.6

4.2 A brief history of the study on a composite material

The study of the properties of random mixtures containing a low fraction of
inclusions has been of much interest since Maxwell’s study of a system of spherical
inclusions.38 The dielectric constant of a composite is given by the low volume fraction
expansion as

e=£0(l+[£]c+0(c2)+...) , (4.1)

where c is the volume fraction of the inclusions. Here [8] is the ﬁrst order coefﬁcient, as

109

 

I

80 *

$

Figure 4.1 A composite material consists of two different dielectric constant regions
(80 , 8,.) in dilute limit.

given by Maxwell, is for a d—dimensional sphere

 

s, + (d —1)so

[81:41:54 ] (4.2)

where [3 , the expansion in brackets, is proportional to the polarisability associated with

an isolated spherical inclusion with dielectric constant 8,. in the host medium with
dielectric constant 80. The 0(c) coefﬁcient for the viscosity of a suspension containing a

system of hard spheres has been also found by Einstein,” 40 although these two do not
have universal relations.

Recently, Thorpe obtained the exact results for the 0(c) coefﬁcient of electrical

conductivity for inclusions of other shapes in 2d using conformal mapping.4|
Hetherington et a].42 calculated the conductivity of a sheet containing n-sided polygons
with sharp comers by solving the integral equations numerically. The results were
compared with analytical solutions in the limits when the inclusion is either a hole or a

superconductor. Jeffrey used the Batchelor-Green multipole expansion formalism to

110

 

determine the second order expansion 0(c2 ) in the problem of dielectric inclusions.43

Binns and Lawrenson44 and more recently Djordjevic et al.45 also studied this problem
with the image charge method in 2d, which provided better convergence than the
multipole expansion method. Choy et al. showed in a recent paper that the troublesome
divergent dipole integral term could be removed and the result agreed completely with the
result of Jeffrey and Djordjevic.46

This study of random mixtures has been done also for the mechanical properties
of composite materials that consist of a host and many inclusions. These studies provided
effective mechanical properties as the function of the rigidity of inclusions, the
distribution pattems,“ and the shapes of inclusion. The elastic moduli of a matrix
containing circular holes with three different distributions of honeycomb, triangular, and
random was studied for the percolation theory of Young’s modulus of the composite.48 It

turned out that the relative Young’s modulus El E0 of a two-dimensional sheer

containing circular holes, overlapping or not, is the same for all materials, independent of
the Poisson ratio, for any prescribed geometry as a consequence of the CLM theorem.49
Davis et a1. studied the elastic moduli of composite materials consisting of an isotropic,
elastic matrix with perfectly rigid circular inclusions near random close packing (see
Figure 4.2).50 This study has been done for the practical reason that the stiffness of rigid
ﬁber inclusions could be enhanced remarkably by arranging them inside a host material.5 1
A theory of the bulk modulus was developed by assuming that the elastic energy of the
neck regions was minimized subject to the constraint that average local strain was equal
to the macroscopic strain. The local strain was approximated by a simple function that

depended on the gap distance but not on the surrounding environment, in a concentrated

lll

    

5'4 48'Q4'14~ 4444 411111411” 33%, 4

414841 #:3114435an 141-44414; 4’.

9’ maﬁmﬁﬁmmﬂmﬂmﬂmmﬂ ﬁg

Figure 4.2 Elastic energy distribution in a composite material with many circular
inclusions under hydrostatic pressure. The contour map was given by the

elasticity analysis using ABAQUS. Elastic energy turned out to be distributed
mostly near the necks between two circles.

112

random array of rigid inclusions in an elastic matrix under hydrostatic deformation in 2d.
In that study, a solution of an electric conduction problem for the geometry was applied to
solve the elastic problem according to its universal formalism.

However, most analytic solutions of composite materials could be obtained only
for simple geometry. For more complicated structures, we should do numerical analysis
using computer simulation. The ﬁnite element method (FEM) would be one of the most
efﬁcient techniques to perform this structural analysis. This technique has been
introduced to determine the isotropic and anisotropic conductivity of single materials in
the previous chapters. However, it can be applied also to a composite material. To test
this method on composite materials, the numerical simulation should be done for a simple

composite where some analytic solutions are available.

4.3 Theory

In this section, I provide the analytic solution of the dielectric problem of a
composite that has two hyperspherical inclusions. Two resistances in the geometry are
deﬁned according to the conﬁgurations of voltage and current measurement. The ﬁrst
configuration is that a voltage is applied on the inﬁnite boundaries of the host. Since the
host is very large compared to the size of the inclusions, the applied voltage forms
constant electric ﬁeld near the inclusions. A resistance deﬁned for this conﬁguration is
called 4-point resistance since there is no common lead for the voltage and current
measurements. The second conﬁguration is that the voltage is applied on the centers of
two inclusions. A resistance from this conﬁguration is called Z-point resistance. The

detailed definitions of these resistances will be given in Section 4.3.1 and Section 4.3.2.

113

 

In both cases, the resistances (or effective conductivities) can be calculated as a function
of the distance between the inclusions. In 2d, exact solutions can be obtained by image
charge method 52 and multipole charge method45 when the inclusions are perfect
conductors. In 3d, an empirical approximation and asymptotic solutions are given for the
similar composite structure. An interesting feature is given by the comparison of these
two resistances. The ratio of the two resistances approaches a constant in the close neck

limit.

 

4.3.1 4-point resistance (R4)

Consider two circles with radius a in the inﬁnite plane as shown in Figure 4.3.
They are separated by L and the neck distance is w. The plane is a normal conductor that

has a conductivity 0'0 and the circular inclusions are perfect conductors (0' >> 0'0 ). Two
inclusions are located in a constant electric ﬁeld Eo formed by a voltage applied on the

boundaries of the host. The voltage between the inclusions is calculated for this geometry.
The initial voltage before the formation of the inclusions is given by EOL. However, if
two circular conductors are implanted in this host, they deform the current paths in the

region near the inclusions. Then, a voltage change AV from the initial value is induced

between the two centers of the inclusions.

4.3.1.1 Exact solutions in 2d
Djordjevic et a1. solved this problem for perfect conductor inclusions in 2d using

the multiple image charge method.45 The solution is given in closed form by

114

L

E0
a

Figure 4.3 4-point resistance measurement of two circular inclusions. The voltage
between two inclusions is calculated in terms of the geometrical parameters
when a constant electric field is applied.

 

 

AV = on/LZ --4a2 = EM/(w+2a)2 —4a2 = E0 w(w+4a) , (4.3)

where L = w+ 2a is the distance between two centers, a is the radius of the circular
inclusions. The asymptotic forms of this solution can be checked simply. AV approaches

EOL as L goes to inﬁnity. No voltage difference is measured (AV = 0) when two

inclusions touch (L = 2a).
Besides this induced voltage AV , we should have a current to deﬁne a resistance

for this conﬁguration. Thus, a normalized current I0 has been deﬁned by total current
flow into the cross section of the inclusions in the direction of E0 (see Figure 4.4). The
current density J 0 in the host is given by (JOE, from Ohm’s law. The cross section of the
inclusion is given by 2a in 2d. Therefore, the normalized current I 0 is given by 2a0'0E0.

Thus, from the given AV and the new current 10, a resistance R4 can be deﬁned by

R. = — . (4.4)

115

 

 

 

 

2a

 

 

 

 

 

 

Figure 4.4 Normalized current I 0 is deﬁned by the total current ﬂowing through the

cross section of the inclusions.

This resistance is called a 4-point resistance because it is equivalent to a resistance

measured from the 4-point measurement. The 4-point resistance R4 can be given using

the closed form of AV in Eq. (4.3) and the normalized current 10.

4 10 2116050 -200 Z a

R _AV_E0‘/W(w+4a) _ l w(w+4] (45)

The asymptotic forms in two limits are given by

R4z-l- K , w/a<<1 (4.6)
0'0 lie
1 w

R4 =—(—+2] , w/a>>1. (4.7)
200 (1

4.3.1.2 Empirical approximation and asymptotic solutions in 3d
Before discussing the 3d composite, let us consider the 1d composite that has two

perfect linear conductors as in Figure 4.5. The current density between the inclusions are

116

 

 

 

 

 

 

 

 

 

Figure 4.5 One dimensional composite material with two perfect conductor inclusions.

identical for the two resistance measurements (R4 , R2 ). Since the voltage AV between the
inclusions is simply given by Eow in the constant electric ﬁeld, the 4-point resistance R4

can be represented by

R, =33: EOW =31 . (4.8)
IO OOEO o.0

 

Let us move to the 3d composite structure. In 3d, the inclusions are spheres.

Thorpe suggested an empirical approximation of the voltage AV between the inclusions.53

l

d 274
AV = EOL[l—(ZL£) ] , (4.9)

where d is the dimension. This formula reproduces correctly the exact solutions in 1d and

2d.

AV=EO[L—2a]=Eow , d=l (4.10)

AV=EO\/L2--4a2 , d=2 (4.11)

In 3d (d = 3), the approximation gives

AV = EOL"“(L3 —8a3)”“ . (4.12)

117

 

In addition, the normalized current I0 deﬁned in Section 4.3.1 is given by Irazcr'oE0 in 3d

because the cross section of the spheres is 71212 . Then, the 4-point resistance R4 in 3d is

obtained by

1/4 3_ 31/4 “4 3 1/4
R4=AV=E°L (f 8“) z 1 3+2 1+2 —8 . (4.13)
10 no COED trooa a a

 

 

Meanwhile, the asymptotic solution of the same composite in 3d was given in the
close neck limit (w/ a <<1).46 That was done by Choy et a1. using the multipole image
method. Following this method, the voltage AV in that limit is obtained as a logarithmic

function.

2 -1
AV 4 ” “E0 , w/a «1 (4.14)

 

 

 

a
ln{—)+ const
w

A constant term exists in the formula as the next leading term. Then, the 4-point

resistance is represented by

—1

AV 7‘ , w/a<<l. (4.15)

R = z
4 I0 300a

 

 

 

a
1n — + const
w

Comparing with this logarithmic form in the limit, the approximation in Eq. (4.6) turns

out incorrect; because the asymptotic form of the empirical formula is given by

R

 

_ ﬂ .. 1502a(3w/a)“‘ _ 2.632 [3

1/4
2 , w/a<<1. (4.16)
IO M 00E0 com

a

4.3.2 2-point resistance (R2)
The second conﬁguration for the resistance calculation is much simpler than the

4-point resistance. As presented in Figure 4.6, a voltage difference V0 is applied across

118

 

 

 

W

Figure 4.6 2-point resistance measurement of two circular inclusions. A resistance is
obtained when the current ﬂows between these two inclusions. I:

 

the centers of two inclusions, and then current flows between the inclusions. Most current

...-L

ﬁows through the neck when the inclusions are close. From this current ﬂow, a resistance
is simply given by

V0

R2=T,

(4.17)

where V0 is the applied voltage and I is total current between the inclusions. This

resistance R2 is called 2-point resistance because it is equal to the normal resistance

measurement.

4.3.2.1 Exact solutions in 2d

Exact solutions of this resistance for this geometry has been derived by the image
charge method. Details of this method can be found in the reference.52 Suppose 8 is the
electric permeability of the host. Then, the capacitance C between two circular conductors

in 2d is given by

119

2 —1 2 -1
C = 27:.e[cosh‘"[2l‘7 4]] = 271:13|:cosh"[2W—a2 +303 “H , (4.18)

where L = w + 2a . The resistance R between two circles can be represented as the

function of the geometry parameters (w, a) by

 

2
19:11: 1 cosh"' W—2+—21+1 . (4.19)
C 0'0 271'0'0 2a a

The asymptotic forms of this resistance can be given by

 

 

 

1
R2 == 1/3 +... , w/a «1, (4.20)
7:00 a
R2 = 1 ln[3)+... , w/a >>1. (4.21) ...—
rroo a

One of these asymptotic formulas (Eq. (4.20)) can be obtained more simply using the

approximation that all currents ﬂow through the neck. It is presented in Appendix A.3.

4.3.2.2 Exact solutions in 1d and 3d

In 1d, the 2-point resistance (R2) is simply given by
R2 = —, (4.22)

where w is the neck distance. This resistance (R2) is equal to the 4-point resistance (R4)
in Eq. (4.8).
In 3d, this 2-point resistance can be obtained by the image charge method.52

-1
1 e °° 1

R =———= 27rd sinh -————— , 4.23

2 C00 00 (misinhmﬁ) ( )

n=l

120

where cosh(,B) 2 2i +1 . The capacitance C between two perfect spherical conductors is
a

 

1 l
C: 2 h
Rea sin ””2 :l[slinh( 2n - l)[3 +sinh( Znﬁ {I

(4.24)
1 l

_2 h =2 1
7128618111 (H)I;— =lsinh(nl3) ”8% + sinh(2[3) + sinh(3ﬁ) +

 

Meanwhile, the exact asymptotic expansion in the limit has been done by Choy.46 It is

obtained by

 

—1
R2 = 1 [ln[£]+ const] , w/a <<1. (4.25)

4.3.3 The ratio of R4 and R2

The ratio of the 4-point resistance and 2-point resistance can be given by the exact
form in 2d and asymptotic form in 3d. In 1d, two resistances turned out identical from Eq.
(4.8) and Eq. (4.22). So, the ratio in 1d is 1.

In 2d, two resistances were given by the exact forms. It should be noted that the
two resistances (see Eq. (4.6) and Eq. (4.20)) have the identical forms of square root in
the limit of w/ a <<l . The only difference is the pre-factor that is a constant 71'. So, the
ratio is obtained by

R4

=71? , w/a<<1. (4.26)
R2

In 3d, due to the infinite sums in Eq. (4.20), the asymptotic behavior cannot be

determined in the close neck limit (w/ a << 1) in this analytic solution. However, if we

121

 

 

consider only the leading terms of logarithmic form in Eq. (4.15) and Eq. (4.22), the ratio

of the resistances (R4 and R,) can be given by

R4

7:2
— z —— z 3.29 . (4.27)
R, 3

4.4 Computer Simulation

4.4.1 FEM analysis for 2d composite

To calculate the resistances deﬁned in the previous section for a composite,
numerical simulations were performed using the FEM. Modeling to obtain the meshed
geometry was done by AN SYS, and the FEM analysis of electrical conduction problem
was performed by ABAQUS. Details about these procedures can be found in Chapter 2.
In this problem, the geometry itself is not complicated like the diamond crystal. It consists
of a large square host and two circular inclusions. Therefore, the modeling was somewhat
a simple procedure. However, I had some troubles in the meshing while computer models .
were generated over a wide range of the neck distances. In the close neck limit, the ﬁnite
elements near the neck region should have been extremely small. Usually, the meshing in
this limit has been restricted due to the limitation of the computer memory. To generate a
constant electric ﬁeld, two large pads were attached to the boundaries of the host and
these pads were modeled to be perfect conductors in the analysis.

I made a large square host material (L4 =1.0 , L)‘ 21.0) that included two circular

perfect conductors (a = 0.05 ). In addition, two large pads (P) were attached on the

122

 

 

 

 

 

 

 

 

 

 

 

Figure 4.7 The geometry of 2d composite material with circular inclusions. Huge plates
were attached to the host to form a constant electric ﬁeld in 4-point resistance
measurement. Those plates were removed for 2-point resistance
measurement.

   

(a) (b)

Figure 4.8 Contour graphs of voltage distribution near the neck in 2d. (a) is given from
4—point measurement and (b) is given from 2-point measurement.

123

boundaries to form a constant ﬁeld inside the host for the 4-point resistance calculation.
However, these pads were not used for 2-point resistance calculations. The geometry was
meshed by more than 2500 square elements that have four nodes on the comers.
Normally, there were more than 2500 nodes in the meshed models. Many computer
models over a wide range of the neck distances were generated. In the analysis of the 4-
point resistances, unit conductivity was given to the host material and high conductivity
(106) was given to the inclusions and pads. The resistance was calculated as a function of
neck distance. Two figures in Figure 4.8 represent the voltage distributions inside the

sample. Figure 4.8(a) denotes for a 4-point resistance calculation when a voltage is

 

applied on the host boundary and Figure 4.8(b) shows a 2-point resistance when the
voltage is applied on the centers of the inclusions. In the 2-point resistance calculation, it

should be noted that the current distribution is concentrated near the neck region.

4.4.2 FEM analysis for 3d composite

A numerical simulation using FEM has been done for 3d composite. A large brick
was constructed as the host material (L, =1.0, Ly = L2 = 0.6) and two spheres
(a = 0.05) were used for two spherical inclusions (see Figure 4.9). After meshing this
geometry using tetrahedral elements, the FEM analysis was performed over a wide range
of the neck distances. Figure 4.10(a) and (b) show the generated meshes and the
equipotential distributions for two resistance calculations (R4 and R2 ).

The computer models were meshed by tetrahedral elements that had four nodes

and four facets. Normally, much ﬁner meshings (~1500 nodes, ~7000 elements) than the

124

 

Figure 4.9 Three-dimensional composite material with two spherical inclusions. Similar
large bricks were attached to the host for 4-point measurement.

 

 

Figure 4.10 Contour graph of voltage distribution near the neck in 3d. (a) is given from
4-point measurement and (b) is given from 2-point measurement.

 

Figure 4.11 The meshed structure of the spherical inclusions. For a rough meshing, as
these spheres approach together, the shape of the neck is getting less
realistic because of the sharp ends of the spheres.

meshing in the ﬁgure were used in the real analysis. However, the generation of realistic
shapes of the sphere was seriously limited by the computer memory problem. The
spherical inclusions were modeled as polygons as shown in Figure 4.11. These rough
meshings could generate signiﬁcant numerical errors in the simulation, especially in the

close neck limit (w/ a << 1).

4.5 Analysis

4.5.1 2d composite

From the calculations using FEM code (ABAQUS), 4-point and 2-point
resistances were examined over a wide range of the neck distances (0.01 < w/ a < 10.0 ). I
plotted these calculations and compared them with the analytic solutions as in Figure
4.12. Cross (+) in the graph represents the 4-point resistance calculation obtained from
FEM analysis and the solid line indicates the exact solution given by Djodjevic et a1. as
shown in Eq. (4.5). Square (:1) denotes the FEM result of 2—point resistance and the solid

line means the analytical solution obtained by the image charge method. The results

126

4-point and 2-point Resistances in 2d

   
  
    

+ 4-point (simulation)
exact: eq.(4.5)

 

D 2-point (simulation)
------ exact: eq.(4.19)

R (resistance)

 

 

4-point and 2-point Resistances in 2d

 

 

 

1.2-
+ 4-point (simulation)
exact: eq.(4.5)
1.0-
!] 2-point (simulation)
------ exact: .4.1
A 0.8'l eq( 9)
o
o
c
S
.S’.’ 0.6-
U)
9
a: 0.4-
............. 1:}
02‘ E1 ...... El ........
3.13”
0.0.,.,.,.,ﬁ,.1
0 0 O 2 0 4 0.6 O 8 1 O 1 2
x(=w/a)
(b)

Figure 4.12 Numerical solutions of the 4-point and 2-point resistances were obtained for
2d composite using FEM analysis. They agree well in all regions with the
analytical solutions. (b) shows the detail in the close neck limit.

127

from the simulation show good agreement with the analytical solutions. However, the
simulated 2-point resistances in the far ﬁeld limit were observed to be slightly higher than
the analytical solutions. That is the ﬁnite size effect of the host, i.e. the inclusions were
getting close to the boundaries of the host as the neck distances increase. In Figure

4. 12(b), more details in the close neck limit were shown to verify the asymptotic
behaviors. Even using the rough meshes, the numerical calculation could reproduce the

square root law in 4—point and logarithmic law in 2-point resistances in that limit.

4.5.2 3d composite
Figure 4.13 presents the numerical result and analytical solutions in 3d. The cross

represents the 4-point resistance (R,) from the simulation and the solid line denotes the

empirical formula that is suggested by Thorpe. Choy’s solution in the close neck region is
shown by the dashed line. This asymptotic line is drawn by selecting a constant for the
best ﬁt with the numerical result in the region. The empirical formula agrees well with the
numerical results in most part in the ﬁgure. However, serious disagreements were
detected in the close neck limit as shown in Figure 4.l3(b). In the limit, the numerical
simulation, the empirical formula, and the Choy’s asymptotic solution do not agree very
well. The numerical error should be from the rough meshing in that region.

The square mark shows the simulated 2-point resistance (R2) and the solid line is
the analytic solution given by the inﬁnite series sum. The dashed line denotes the another
asymptotic formula in the other limit (w/ a <<1) given by Choy. The simulated results
were obtained slightly lower than the analytic solution. The rough meshing raised errors

again for the 2-point resistances in the close neck limit. Substantial discrepancy between

128

 

4-point and 2-point Resistances in 3d

   
   
  
    

 

 

 

 

 

 

 

 

 

 

 

8° ' + F1, (simulation)
theory: eq.(4.13)
70‘ .. ”asymptotic: eq.(4.15)
6° 4 D R, (simulation)
3 theory: eq.(4.23)
g 50 7 -------------- asymptotic: eq.(4.25)
m
Iii
'7, 40 ~
2
m 30 ..
20-
10
S E El
0 r f ' I T l '
0 2 4 6 8 10
We
(a)
4-point and 2-point Resistances in 3d
15 -
+ R, (simulation)
‘ —-—— theory: eq.(4.13)
’8‘ 101 asymptotic: eq.(4.15)
c 1
g D R, (simulation)
8 theory: eq.(4.23)
b . - ----------- asymptotic: eq.(4.25)
(I
5 -1
1 ..................... a
. ................. u
Jmﬂw‘ LJ
0 f T f r n r r r
0.0 O 2 0.4 0.6 0 8 1 O
w/a
(b)

Figure 4.13 Numerical solutions of the 4—point and 2-point resistances were obtained
for 3d composite using FEM analysis. They agree nicely in most regions
with the analytical solutions, but show some disagreements in an
asymptotic regime (close limit). (b) shows the detail in that limit.

129

 

 

the numerical result and the exact solution has been detected in that limit

4.5.3 The ratio of 4-point and 2-point resistances

The ratio of two resistances (4-point and 2-point) has been one of the major
concerns in this study. The analytic prediction and the numerical result are presented in
Figure 4.14. The ratio can be determined analytically in 2d since all exact forms of the
resistances are obtained. As expected from the analytical solutions, that ratio approaches E,

It in the close neck limit (w/ a << 1 ). The numerical simulation is observed to approach

this constant in that limit. The prediction of this ratio in 3d was given by 7:2 l 3 from the

 

asymptotic forms done by Choy. However, from the numerical simulation, the ratio turns
out to approach ~8.0. A large error occurred here mainly because the numerical error
caused by the rough meshings. However, the analytic prediction also has error because
the constant terms as the second leading terms in the formula have been neglected in the
derivation of this ratio. It should be tested more carefully with better meshings. This study

could be a subject for the future research.

4.6 Conclusion

In this chapter, I have described two main subjects. One was the demonstration of
the numerical analysis using the FEM for the study of composite materials and the other
was the investigation of the ratio of two resistance calculations (4—point and 2-point).
Both studies were done for a simple composite model that was composed of a large host

and two spherical inclusions of perfect conductors. Two resistance measurements (4-

130

 

 

 

20 - O simulation (2d)
----- theory (2d) A
‘ A simulation (3d)
15 "I'I .............. theory (3d) A
11:” 2 A
r 10 -
“C . A
‘ A A A A A
5 4 O 0 O
. O O
q5'3'3‘3'3’3'3’3'é‘2"—"'—":‘3'3'6'3‘3'3’3'3'02':‘B'EG —————— Q -----------------
0 ... . , . 1
0.01 0.1 1 1 0
w/a

Figure 4.14 The ratio between 4-point resistance and 2-point resistance in 2d and 3d.
Each ratio was observed to approach a constant in the close neck limit
(w/a <<1). The limiting value is given by trin 2d and 71:2 /3 in 3d.

131

 

point, 2-point) were deﬁned in the geometry depending on whether the voltage is applied
to the boundary of the host or the inclusions. For the ﬁrst subject, I presented the
simulated results of these two resistances using the FEM analysis and they were
compared with the analytic solutions. Since the numerical results were observed to be
consistent with the analytic predictions, the FEM could be suggested as a primary step to
investigate the effective medium theory of complex composite materials. For the second
subject, I compared two resistances over a wide range of the neck distances. A simple
relation between these two quantities was found analytically in 2d and numerically in 3d
in the close neck limit ( w/ a << 1). However, the numerical result should be tested with
better meshings. I suggest that this simple relation may provide an useful approach to the

study of composite material properties.

132

 

 

APPENDICES

133

A.1 Reciprocity theorem

This electric reciprocity theorem has been stated in many works since the earliest
work was done by Lord Rayleigh in last century.54 Since then, many publications related
to this topic came out because of theoretical and practical interests.55‘56’S7 This relation is
quite general; it works for all materials with any shapes in any dimensions. It does not
depend on whether the sample is homogeneous, isotropic or not. The sample can be a
composite material that has various shapes of inclusions of different conductivity, even
some holes in it. The only condition of this theorem is that the sample should be Ohmic.
First, this theorem was concerned with electric ﬁeld. However, since the role of magnetic
field in a conductivity tensor was understood by the Onsager relation,” 59 this reciprocity
theorem has been concerned with the interference term in the conductivity tensor by
magnetic ﬁeld, too. This leads to the reverse-magnetic-ﬁeld reciprocity (RMFR) under
the assumption of the linearity of a sample.2| These reciprocity relations have important
implications for electronic-transport properties in electric current and magnetic ﬁeld (Hall
effect, magnetoresistance).

In this appendix, I simply derive this theorem for zero magnetic field under the
assumption that a bulk conductor can be represented by an equivalent resistor circuit.
Suppose a sample in arbitrary shape is electrically equivalent to a circuit of resistors as in
Figure A. 1. All of the resistors are fully connected to each other. In this proof, the sample
and the resistors in the circuit are assumed electrically Ohmic, which means that it does

not show any diode behavior. This is the only requirement of the theorem. Thus, this

134

 

1T

 

 

 

 

 

 

Figure A.1 A sample with four probes and its equivalent circuit. The electric circuit is
composed of many resistors that are fully connected.

theorem does not work for non-Ohmic sample. The laws of Kirchhoff in an electrical

 

circuit can be used to prove it.

Suppose we have a sample with n leads on it for the conduction, then there are r
(="C2 ) resistors in the circuit. First, we need to know the number of independent current

loops in this circuit to ﬁnd out how many linear equations we need. The maximum
number of currents is r, which is same with the number of resistors. And there are n —1
constraints of currents by current conservation on each lead, known as the ﬁrst law of
Kirchhoff. Although there are n vertices for constraints, the two constraints on vertices
that the external current ﬂows in and out are really a single constraint. Then we have d
independent current loops for n leads as

d=r—(n—l)=r—n+l, (A.1)
where r (="C,) is a number of resistors and n is a number of vertices. Because there are r

unknown currents, i], i,, , i, on each resistor, we have n —1 constraints for currents at

 

each vertex by Kirchhoff’ s ﬁrst law of current conservation. Although there are n vertices

135

for constraints, the two constraints on vertices that the external current flows in and out
are really a single constraint. Since r is equal to n(n — l)/2 , the number of independent
loops in the circuit is given by

___ n(n—1)_(n_1)=(n—1)(n—2)
2 2

d (A-Z)

 

which is one for 3 leads, three for 4 leads, six for 5 leads and so on. We also have one

more loop by the external current i0 as well as d loops in the sample. We need to set up

(1 +1 current 100ps for Kirchhoff laws which give (1 +1 current equations. There is one
condition in setting up these loops that each resistor should be included in at least one

current loop. To prove the reciprocity theorem, a resistor R M is selected for external
current ﬂow and Rm” for voltage measurement. A voltage V is applied across RM and the
current io ﬂows in. Suppose Rm is shared by m current loops( i 1+1 , .. , i 1+». ) and RM is

shared by p current loops (i,+l , , i,er ) and the external current i0. Now we have d +1

linear equations from Kirchoff laws

 

 

 

 

/ Rm 0 . 0 0 . _qu —qu . 0 V 10 1 (v)
0 ML] . . . . . . . . 1, 0
O . . Mj+Lj+l . . . . . . ij+1 0 (A 3)
O . . . Mj+ij+m . Mj+m’k+l . o . ij+m = O 0

- qu - - - M k+l, j+m - M k+l,k+l - 4 4 i4+1 0

-qu . . . . . . Mk+pJ¢+P . . ik+p 0

L 0 . . . . . . . . Md,d A id 1 \0/

 

where M00 is R andpterms( M0,k+l’ .. , M0.k+p) in 0-th column are —Rm. In Eq.

I"!

(A.3), M 1. j is a matrix element that represents a sum of resistors where the current i j

136

 

flows in the i-th loop. It is equal to M 1.1 because those two terms represent the same
resistors that are shared by the i-th and j-th current loops. A solution for i j from Eq. (A.3)

is given by

 

(A.4)

where A denotes the determinant of the matrix and A 1. j is a cofactor of M 1. j. . Since the

net current flow 2 i

0:1

in Rm" is given by the sum of all currents of loops pass through it,

j+a

a resistance deﬁned by a voltage across m, n and a current through p, q is given by

m m
Rmn Zij-l-a Rmn 2A0,j+a
an _ =1 =1

i0 V A0.0 A0,0
A

Rm"-

pq _ , (A.5)

 

In this equation, AOIM is equal to _ quiAw’ k +b).(0. I”) because there are p terms of
b=l

— R m in the 0-th column and the other terms are zero except 0-th element R W irrelevant
in calculating A0. 1+4 . [Amman is the determinant of the matrix without the 0-th and k-th

rows and the 0-th and j-th columns] From this relation, we have

m
Rmanq 2 i A(O.k+b),(0,j+a)
Rmn ___ _ a=lb=l , (A.6)
Pq A0,0

 

By switching the leads p, q and m, n [which is equivalent to exchanging k +b , j+ a in

this notation], we have another resistance Rf: deﬁned by a voltage V applied across Rm

137

 

 

 

instead of R124 and a voltage measured across RM . Since A(0.k+b).(0. 1+4) 18 equal to

A“, 1.,“ 110.1412) due to the symmetry of the matrix, we have a symmetric relation by

m
quRmn 2 i A(0,j+a),(0,k+b)

R53, = - “=‘b=‘ = Rm” (A.7)
ADJ) Pq

 

This reciprocity theorem says that the resistance from a voltage across p, q and
current through m, n is equal to another resistance from voltage across m, n and current

through p, q in an equivalent resistor electrical circuit, and therefore in the sample.

A.2 Counting of Independent 4-terminal terms

In this part, I describe how to count the number of 4-terminal independent terms
in a resistance table. I have shown that there are "C, ( = n(n — 1)/ 2) independent terms in

a resistance table for an n-probe measurement in Section 2.2.2. Four probes are enough to
determine the isotropic conductivity because the four probes generate two independent 4-
terminal terms. However, the unknown parameters are increased in an anisotropic
conductivity tensor. Thus, we need more independent terms in a resistance table. To have
more independent terms in a table, we need to have more probes on the sample. How
many probes are needed for an anisotropic conductivity tensor? In Table A. l , the ﬁrst
column shows the number of terminals and the second column is the number of
unknowns when the contacts are Ohmic. The third column is the number of independent
terms in a resistance table as proved in Section 2.2.2. The fourth column is the number of

contact resistances that is equal to the number of leads.

138

1

Table A.1 A table represents the number of independent 4-terminal terms in a resistance
matrix for n-probe measurement.

 

 

 

 

 

 

 

 

 

 

 

Number of Unknowns Independent Contact Independent 4-
Terrninal ( P1,- , contact Elements resistances terminal terms
resistances) (A) (B) ( A - B )
2 8 1 2 0
3 9 3 0
4 10 6 4 2
5 l 1 10 5 5
6 12 15 6 9
7 13 21 7 14
8 14 28 8 20
n 6+n "C, (=n(n_1)/2) n n(n—3)/2

 

 

 

 

 

139

 

 

The number of independent 4-terminal terms in the last column can be given by
following argument. There are A independent terms in a resistance table. Meanwhile, the
contact resistances are independent of the sample resistances. They are independent each
other, too. So, if we subtract the contact resistances (B) from the whole independent terms
(A), the rest (A -B) should be independent of the contact resistances. Therefore, we have
A-B independent 4-terminal terms. Following this argument, we need 6 probes on a II
sample to determine 6 components in an anisotropic tensor, because six probes generate 9

independent 4-terminal terms.

 

A.3 Resistances of Two Perfect Hyperspherical

Conductors

Consider two hyperspheres in the plane with radius a and separated by L so that
the neck distance is w (see Figure A.2). Ifx is the distance up the vertical line from the

center, then we have

x: (a+-’3)sin0 = (1+l)asin0 = (1+1)asin0 , x =l. (A8)
2 2a 2

a
The two spheres have a resistance between them that is mainly given by the neck region.
If we strip the neck into long this differential sections and add up the resistances in

parallel, we ﬁnd the total resistance R as [where w = L - 2a].

 

=00Lxmax Yd-l'x d—de

(A.9)
w + 2a(1— c086)

where 7d factor for the area in d dimensions, that might be 2 in 1d and 271’ in 2d and 471'

140

 

 

 

Figure A.2 Two circular perfect conductors in a host material in the close neck limit.

_d_
in 3d etc. More generally “yd = 27: 2 /I"(-‘-§-). If we look at small 6, we may write

9 z x/ a , so that the integral becomes

d-2
xmax Yd-lx dx

-1
R =0
00

 

2 (A.10)
w+ x lo

if we make a change of variable to y, where x = Jaw y , the integral can be rewritten as

which can be evaluated to give

d__; d-2d
R 1=0'0ayd_1(ara) 2 y y

 

2 (A.11)
1+y

where We have replaced the upper limit in the integral with inﬁnity for narrow necks.
This integral converges at the lower end if d >1 and at the upper end if d < 3 , so that in
general for 1 < d < 3, we have the result

and—3341)
3-d '

R“1 = 0'0 (A.12)

(war) 2

This gives a result in two dimensions as

141

 

R"l =007Na/w (A.13)
Note that as d approaches one from above, R‘1 = 0', / a) as expected, and as d approaches

three form below, we have a divergence. If we use the closed form in 2d,

 

 

R—I=O'0J’xmax yldx :0 Jam,“ cosBdG (A.14)
0 w+2a(1-c089) 0 (l+x)—cosB

Then we have

9m...
(1 1 2 —1
20 +12) arctan ( +1) ° 1.311%
,1 _ 2
(1+1) 1 0 (A.15)

2(l+x) arctan (1+1)

200(1+x) —6max +
(/(l+x)2 —1 2C

In the limit 1 = w/ 2a <<1, we have the same result with Eq. (A. 13).

2 73 6m} ,_

R"1 z 00['9max +—arctan{———~tan = 0'07: 1 (A.16)

71'
727 72 T ““74 w

A more close examination of the integral when d=3, shows that if we use xmax = a for the

R"1=oo(1+x) —0+

 

 

 

I

 

 

 

upper limit, which seems reasonable, we can do the integral in Eq. (A9) in closed form,

to give

 

 

R—1_O. Jxmax szdx _750'0 xmax xdx
— 0 —
0 w+2a(l-cos6) a 0 (l+x)—cos0

=7“, (1+X)zajgmu cosGsin 9d9
0 0 (l+x)—cosO

(A. 17)

 

Let y = cos6 , then

142

0 J1 (1+x)-y

=4000+x)2al—y—(1+x)1n{(1+x)-y}”I“:

=7t'0'0 (1+Z)Za[(ymax _1)+(l+x)ln{(1+X)x_ Ymax }]

== trooa[(ymax — 1) + 1n{1—_—y—m-9-’-‘-}]
Z

(A.18)

 

If we take ymax = 0 (0 :g)’ then

R" z naoa[1n{i}—1]= naoa|:ln{2—a} — 1] = nooa[ln{—a—} + 1n{2}— 1]
7‘ W W (A.19)
= moa[1n{3} — 0.3069]
W

And if ym = =0.7071 (9=%,xm=a), then

1
72
0.2929

x

 

R" =—. trooa[- 0.2929+ ln{

= naoa[1n{3—} —O.8277]
W

The exact asymptotic expansion done by Choy gives

R‘1 z 00a7r[ln(£]— E] z 00aﬂ|:ln[£]— 2.667] . (A.21)
a) 3 60

Thus we get the leading term right in both 2d and 3d, but the next term is wrong.

H = naoa[1n{i}+1n 2 + 1n 0.2929 42929]
W (A20)

143

 

BIBLIOGRAPHY

144

 

BIBLIOGRAPHY

 

' The term resistivity is also used sometimes in this dissertation.
2 This sample was prepared by M. Schnars and B. Golding.

3 S. M. Sze, Physics of Semiconductor Devices (John Wiley & Sons, Inc., New York,
1981), pp.848-851.

4 R. J. Gillespie et al., Chemistry (Allyn and Bacon, Inc., Newton, MA, 1986), p.803.

5 J. R. Olson, R. O. Pohl, J. W. Vandersande, A. Zoltan, T. R. Anthony, and W. F.
Banholzer, Phys. Rev. B 47, 14850 (1993).

6 c. H. Seager and T. G. Castner, J. Appl. Phys. 49, 3879 (1978)

7 J. Werner, “Electronic Properties of Grain Boundaries,” in Polycrystalline
Semiconductors : Physical Properties and Applications, ed. By G. Harbeke, (Springer-
Verlag, Heidelberg, 1985), pp.76-94.

8 L. J. van der Pauw, Phillips Techn. Rdsch. 20, 230 (1958).

9 L. V. Bewley, Two-dimensional ﬁelds in electrical engineering (The MacMillan Co.,
New York, 1948).

'0 Y. Sun, 0. Ehrmann, J. Wolf, and H. Reuchl, Rev. Sci. Instrum. 63, 3752 (1992).

” Any boundary value problem can be classiﬁed into well-posed or ill-posed problem.
Following is a deﬁnition of those problems by Hadamard.ll
Well-posed problem :
(a) There exists a (globally-deﬁned) solution for all (reasonable) data.
(b) The solution is unique.
(c) The solution depends continuously on given data (stability).
Ill-posed problem :
All problems do not satisfy these conditions of well-posed problems.

'2 A. N. Tikhonov and V. Y. Arsenin, Solutions of Ill-Posed Problems (Wiley, New York,
1977).

[3 S. Kubo, in Inverse problems in engineering mechanics, edited by M. Tanaka and H.
D. Bui (Springer-Verlag, Heidelberg, 1993), p.51.

'4 A. R. Day developed an algorithm for this using the conjugate gradient method.

145

 

 

'5 L. Collatz, The numerical Treatment of Diﬁ‘erential Equations (Spring-Vet lag, New
York, 1966).

16 S. H. Crandall, Engineering Analysis (McGraw-Hill, New York, 1956).

'7 K. J. Bathe and E. L. Wilson, Numerical Methods in Finite Element Analysis (Prentice-
Hall, Inc., Englewood Cliffs, NJ, 1976))

'8 G. Strang and G. Fix, An Analysis of the Finite Element Method (Prentice-Hall, Inc.,
Englewood Cliffs, NJ, 1973).

'9 D. S. Burnett, Finite Element Analysis .' From Concepts to Application (Addison-
Wesley, Reading, MA, 1987).

20 B. I. Bleaney and B. Bleaney, Electricity and Magnetism (Oxford University Press,
Oxford, 1957), p.69

21 H. H. Sample, W. J. Bruno, S. B. Sample, and E. K. Sichel, J. Appl. Phys. 61, 1079
(1987).

22 C. Young et al., Ray Tracing Creations (Waite Group Press, Corte Madera, 1994), p.3.

23 W. H. Press, S. A. Teukolsky, W. T. Vetterling, and B. P. Flannery, Numerical Recipes
(Cambridge University Press, New York, 1992), p.653.

2“ ANSYS is a registered trademark of SAS IP, Inc..

25 ABAQUS is the product of Hibbitt, Karlsson and Sorensen, Inc.. For more information,
see ABAQUS User’s Manual HKS, Inc..

26 For the labeling scheme chosen, the contact resistances r,‘ and r4 do not appear in

Table 2.5. A labeling scheme that involves the cyclic permutations of the four indices (
i.e. which replaces 114 with 141 ) would be more useful. The resistance table would
contain all 8 contact resistances that could then be determined from a single set of
measurements.

27 M. D. Jaeger, s. Hyun, A. R. Day, M. F. Thorpe, and B. Golding, Diamond and Related
Materials 6, 325—328 (1997).

28 B. I. Shklovskii and A. L. Efros, Electronic Properties of Doped Semiconductors,
(Spring-Verlag, Berlin, 1984), Chap.4.

2" M. D. Jaeger, Ph. D. Thesis, Michigan State University, 1997.

146

 

 

 

30 C. Kittel, Introduction to Solid State Physics, 6th ed. (John Wiley & SOUS. Inc's New
York, 1986), p.19.

31 N. W. Ashcroft and N. D. Merrnin, Solid State Physics (Pergamon Press, Inc., New
York, 1962).

32 W. P. Mason, Crystal Physics of Interaction Processes (Academic Press, New York,
1966), 13.216.

33 L. D. Landau and E. M. Lifshitz, Electrodynamics of Continuous Media (Pergamon
Press, New York, 1960), p.93.

34 H. R. Schwarz, Numerical Analysis (John Wiley & Sons, New York, 1989), p.320.
35 z. Xiong and A. Kirsch, J. of Comp. and Appl. Math., 42, 109 (1992).

36 G. Arlken, Mathematical Methods for Physicists, 3rd Cd- (Academic P1933, Orlando,
1985), p.510.

37 This sample was prepared by M. Schnars and B. Golding.

38 J. C. Maxwell, Electricity and Magnetism (Clarendon Press, Oxford, 1873).

39 A. Einstein, Ann. Phys. 19 289 (1906).

40 A. Einstein, Ann. Phys. 34 591 (1911).

4‘ M. F. Thorpe, Proc. R. Soc. London Ser. A 437, 215 (1992).

42 J. H. Hetherington and M. F. Thorpe, Proc. R. Soc. London Ser. A 438, 591 (1992).
43 D. J. Jeffrey, Proc. R. Soc. London A 335, 355 (1973).

44 K. J. Binns and P. J. Lawrenson, Electric and Magnetic Field Problems ( Pergamon,
Oxford, 1973), pp.48-53.

45 B. R. Djordjevic, J. H. Hetherington, and M. F. Thorpe, Phys. Rev. B 53, 14862
(1996).

46 T. c. Choy, A. Alexopulos, and M. F. Thorpe, Phil. Trans. R. Soc. Lond. A (1996).

47 J. Chen, M. F. Thorpe, and L. C. Davis, J. Appl. Phys. 77(9) 4349 (1995).

147

 

4“ A. R. Day, K. A. Snyder, E. J. Garboczi, and M. F. Thorpe, J. Mech. Phys. Solids, 40
1031-1051 (1992).

49 A. Cherkaev, K. Lurie, and G. W. Milton, Proc. R. Soc. London A 438, 519-529
(1992).

50 L. C. Davis, K. C. Hass, J. Chen, and M. F. Thorpe, Proceedings of the Conference on
Composite Materials held in Virginia, June 1993.

5' M. F. Thorpe and I. Jasiuk, Proc. R. Soc. London A 438, 531 (1992). r;

52 W. R. Smythe, Static and Dynamic Electricity, 2"d ed. (McGraw-Hill, New York,
1950). 1176-78, p.118-122.

 

53 M. F. Thorpe, unpublished.

:I!-I' -'I Ali-11¢“ 3.“, ’3'-‘ F ...‘r a '

54 J. W. Strutt, Lord Rayleigh, The Theory of Sound, 2“d ed., vol. I (Dover, New York
1945; reprinted from the 2"d ed. Published by Macmillan in 1894), p. 155.

55 H. A. Lorentz, Collected Papers, vol. III (Martinus Nijhoff, The Hague, 1936), pl.
56 J. R. Carson, Bell Syst. Tech. J. 3, 393 (1924).

57 v. H. Rumsey, Phys. Rev. 94, 1483 (1954).

58 L. Onsager, Phys. Rev. 37, 405 (1931).

59 L. Onsager, Phys. Rev. 38, 2265 (1931).

148

   

MICHIGAN STATE UNIV. LIBRARIES
ll11111111111111111111111111111”ll1111111111111111111111
31293016887519