Measurement of the phase space distribution in a complex plasma
by
Ross Kenney Fisher
A dissertation submitted to the Graduate Faculty of
Auburn University
in partial fulfillment of the
requirements for the Degree of
Doctor of Philosophy
Auburn, Alabama
May 6, 2012
Copyright 2012 by Ross Kenney Fisher
Approved by
Edward Thomas, Jr., Chair, Professor of Physics
James Hanson, Professor of Physics
Allen Landers, Associate Professor of Physics
David Maurer, Associate Professor of Physics
ii
Abstract
The development of a diagnostic technique that allows measurement of the six
dimensional phase space distribution for the dust component of a plasma system is presented. In
the course of the measurement and analysis processes a number of longstanding questions
related to the basic properties of dusty plasmas are addressed and explained for the first time.
The work described below demonstrates the importance of using a generalized form of the
standard Maxwellian probability distribution function to model the velocity space portion of
phase space. This generalization, the trinormal probability distribution function, allows
ellipsoidally symmetric anisotropy in velocity space; the appearance of such anisotropy has been
an outstanding issue in weaklycoupled dusty plasmas for several years. The measured velocity
space anisotropy is shown to be both a real effect and an effect that is too large to be accounted
for as a perturbation to the standard spherically symmetric model of the velocity space. The
spatially resolved measurements within the dust cloud and the new model for the velocity space
are then combined to give the spatial distribution of the fluid properties of the dust cloud for the
first time in a weaklycoupled system. The fluid thermodynamic and transport properties
obtained through the process are discussed in detail. The analysis of these newly available
transport and thermodynamic properties clearly shows that the dust component of the system is
in a state of dynamic forcebalanced equilibrium. The dust component of such systems has long
been suspected to be in a forcebalanced equilibrium; the appearance of the trinormal
distribution in velocity space unambiguously confirms the suspicion. The fact that the system is
in a state of dynamic equilibrium is demonstrated by examination of the transport properties of
the dust component and has not been previously demonstrated experimentally.
iii
Acknowledgments
I would like to express my thanks to Dr. Edward Thomas for his support. Without his
guidance the work I have done over the last several years would not have been possible. His
combination of sincere interest in this work and his willingness to let me ?follow my nose? (even
when it must have looked like I was heading towards a dead end) was, and still is, greatly
appreciated.
The many discussions I have had with Dr. Jeremiah Williams provided a great deal of
insight into the diagnostic used in this work. His development of the stereoPIV system during
his time at Auburn laid the foundation for much of this work. His patience with my many
questions is commendable.
I would also like to thank my undergraduate advisor, Dr. Robert Merlino, for giving me
the opportunity to work in his laboratory and for sparking my interest in studying dusty plasmas.
His willingness to continue acting as an unofficial ?outside/objective advisor? has been very
useful and is very much appreciated.
The support and help from the other students working under Dr. Thomas (specifically
Ami Dubois and Ashley Eadon) has been invaluable. I would also like to acknowledge the
faculty in Auburn?s physics department for the assistance and direction that they were willing
and able to provide.
Most importantly, I would like to thank my parents, Bill and Ann Fisher, for their infinite
support and understanding during my time in graduate school. Without their encouragement,
belief, and without the work ethic that they passed on this process would have proven much
more difficult, if not impossible.
iv
Table of Contents
Abstract ......................................................................................................................................... ii
Acknowledgments........................................................................................................................ iii
List of Tables ............................................................................................................................. viii
List of Figures .............................................................................................................................. ix
Chapter 1: Introduction ................................................................................................................ 1
1.1: Historical overview ................................................................................................... 1
1.2: Dusty plasma parameters .......................................................................................... 2
1.2.1: Plasma parameters ..................................................................................... 3
1.2.1.1: Quasineutrality .......................................................................... 3
1.2.1.2: Debye length ............................................................................... 4
1.2.1.3: Plasma frequency ........................................................................ 6
1.2.1.4: Coulomb coupling parameter...................................................... 7
1.2.2: Dust grain charging .................................................................................... 9
1.2.2.1: Charging processes overview ..................................................... 9
1.2.2.2: Dust charging due to electron and ion collection ..................... 11
1.3: Motivation and scope .............................................................................................. 16
Chapter 2: Theory and Background ........................................................................................... 18
2.1: The phase space distribution ................................................................................... 18
2.2: Velocity space distribution models ......................................................................... 22
v
2.2.1: The Maxwellian distribution function ..................................................... 23
2.2.2: The multinormal distribution function ................................................... 26
2.2.3: The principal axis coordinate system....................................................... 30
2.2.4: Other generalized distribution functions .................................................. 34
2.2.5: Applicability of the trinormal distribution function model .................... 37
2.3: Bayesian probability theory .................................................................................... 40
2.3.1: Introduction to Bayesian probability theory ............................................ 40
2.3.2: Bayesian model selection......................................................................... 41
2.3.3: A simple example .................................................................................... 44
2.3.3.1: Penalizing for unknown distribution parameters ...................... 45
2.3.3.2: Finite measurement error .......................................................... 46
2.3.3.3: Calculation of the parameter values.......................................... 48
2.3.3.4: Summary ................................................................................... 50
Chapter 3: Experimental Apparatus ......................................................................................... 52
3.1: 3DPX ...................................................................................................................... 52
3.1.1: 3DPX subsystems ................................................................................... 52
3.1.1.1: Vacuum vessel .......................................................................... 52
3.1.1.2: Gas flow .................................................................................... 54
3.1.1.3: Plasma source ........................................................................... 55
3.1.1.4: The dust..................................................................................... 57
3.1.2: Description of the plasma environment ................................................... 58
3.2: Particle Image Velocimetry .................................................................................... 60
3.2.1: Introduction .............................................................................................. 61
vi
3.2.2: Two dimensional PIV .............................................................................. 62
3.2.3: Stereoscopic PIV...................................................................................... 65
3.2.4: Stereoscopic PIV in 3DPX ...................................................................... 70
Chapter 4: Phase space distribution measurements ................................................................... 74
4.1: Introduction to the experiment ................................................................................ 74
4.2: Single volume element............................................................................................ 75
4.2.1: Configuration space ................................................................................. 76
4.2.2: Velocity space .......................................................................................... 77
4.2.2.1: Calculation of the drift velocity vector ..................................... 77
4.2.2.2: Single particle velocity space ................................................... 78
4.2.2.3: Visualization of the velocity space distribution ........................ 80
4.2.2.4: Parameter value calculation: Maxwellian distribution model . 82
4.2.2.5: Parameter value calculation: Trinormal distribution model ... 84
4.2.3: Comparison of the velocity space distribution models ............................ 86
4.2.3.1: The posterior probability ratio .................................................. 87
4.2.3.2: The effect of large data sets ...................................................... 90
4.2.4: Summary .................................................................................................. 93
4.3: Velocity space model selection............................................................................... 94
4.3.1: Posterior probability ratios....................................................................... 94
4.3.2: Geometrical comparison .......................................................................... 96
4.3.2.1: Velocity space volume .............................................................. 96
4.3.2.2: Variance ratios .......................................................................... 98
4.3.2.3: Geometrical connection ............................................................ 99
vii
4.4: The dust component phase space distribution ...................................................... 101
4.4.1: The number density ............................................................................... 101
4.4.2: The drift velocity ................................................................................... 104
4.4.3: The velocity space covariance tensor, laboratory coordinate system .... 107
4.4.4: The velocity space covariance tensor, principal axis coordinate
system ................................................................................................. 109
Chapter 5: Transport and Thermal Properties of the Dust Component ................................... 117
5.1: The fluid transport equations ................................................................................ 118
5.2: The continuity equation ........................................................................................ 125
5.3: The momentum equation ...................................................................................... 129
5.3.1: Energy densities ..................................................................................... 129
5.3.2: Force densities ....................................................................................... 136
5.4: The energy equation.............................................................................................. 142
5.5: Further discussion of the transport measurements ................................................ 146
5.5.1: A dynamic equilibrium .......................................................................... 147
5.5.2: The heat flux .......................................................................................... 149
5.5.3: Summary ................................................................................................ 152
Chapter 6: Conclusions ............................................................................................................ 154
6.1: Summary ............................................................................................................... 154
6.2: Future work ........................................................................................................... 157
6.3: Concluding remarks .............................................................................................. 158
References ............................................................................................................................... 159
Appendix A: Finite difference derivatives ............................................................................... 162
viii
List of Tables
Table 3.1: Approximate plasma and dust parameters. .............................................................. 59
Table 4.1: PSD parameters with the Maxwellian velocity space model. ................................... 84
Table 4.2: PSD parameters with the trinormal velocity space model. ..................................... 86
Table 5.1: Velocity space integral results with the trinormal velocity distribution function. 119
Table 5.2: Transport equations summary. ............................................................................... 125
Table 5.3: Continuity equation quantities and corresponding SI units. ................................... 126
Table 5.4: Quantities found within the momentum equation and the corresponding SI units. 129
Table 5.5: Average energy density ratios. .............................................................................. 136
Table 5.6: Quantities found within the energy transport equation and the corresponding
SI units .......................................................................................................................... 143
ix
List of Figures
Figure 1.1: Dust clouds with different phases. ............................................................................ 8
Figure 1.2: Schematic of a grazing collision between an electron or ion and a dust grain........ 12
Figure 1.3: Log plot of the dust grain charge number Z
d
, in electrons, vs the ratio of dust
number density to ion number density. ........................................................................... 15
Figure 2.1: Plots of the one, two, and three dimensional Maxwellian distribution function. .... 25
Figure 2.2: Plots showing the binormal distribution function. ................................................ 28
Figure 2.3: Plots showing the surface of constant probability amplitude for the trinormal
distribution. .................................................................................................................... 29
Figure 2.4: An example of a rotation to the principal axis coordinate system in two
dimensions. ..................................................................................................................... 31
Figure 2.5: An example of a rotation to the principal axis coordinate system for the trinormal
distribution function. ...................................................................................................... 33
Figure 2.6: The hierarchy of the Pearson family of distribution functions. . ............................. 35
Figure 3.1: (a) Photograph of the 3DPX device. (b) Schematic view of the connected vacuum
crosses. (c) Detailed schematic of the "experiment end" of 3DPX. ............................. 53
Figure 3.2: (a) and (b) Photographs of the anode inside of the vacuum vessel in the position used
for the experiment. (c) and (d) Photographs of the anode outside of the chamber. (e)
Circuit diagram for the dc glow discharge plasma source. ............................................ 56
Figure 3.3: Schematic drawings of the 2DPIV configuration for a single particle. ................. 63
Figure 3.4: Flow chart for a single 2DPIV calculation for an interrogation cell containing
multiple tracer particles. ................................................................................................. 64
Figure 3.5: (a) Top view of the stereoPIV camera and laser sheet orientation. (b) The
displacement from (a) as seen by camera 1. (c) The displacement in part (a) as seen by
camera 2. ........................................................................................................................ 66
x
Figure 3.6: Images of the calibration plate and the correction with the Scheimpflug adapter. . 67
Figure 3.7: (a) Full camera frame view of a dust cloud cross section. (b) Detailed view of the
camera region containing the dust cloud (the region outlined in yellow in (a)), the yellow
grid indicates the location of the individual interrogation cells. (c) The four camera
images of the cell outlined in red in (b). (d) The instantaneous velocity field measured
for the cloud cross section, projected onto the image plane. ......................................... 71
Figure 3.8: Time series showing the first ten three dimensional velocity vector measurements for
the cell highlighted in Figure 3.7 .................................................................................... 72
Figure 4.1: Number density determination. .............................................................................. 76
Figure 4.2: A plot of the correction factor, CF(N
d
). ................................................................. 79
Figure 4.3: Various histograms of the measured velocity vectors. ........................................... 81
Figure 4.4: Various histograms of the error measurements. ..................................................... 82
Figure 4.5: Plots showing the order of magnitude of the K ratio as a function of the number of
vectors randomly selected from D. ............................................................................... 91
Figure 4.6: Plots showing the parameter and uncertainty values obtained from random samples
of D as a function of the number of vectors within the sample. .................................... 92
Figure 4.7: Histograms of the base 10 logarithm of the K ratio for all volume elements within
the dust cloud. ................................................................................................................ 94
Figure 4.8: Histogram of the order of magnitude of the penalty factor assessed to the
distributions for their unknown parameters and uncertainties. ..................................... 95
Figure 4.9: Scatter plots comparing the velocity space volumes calculated with both models for
each volume element in the dust cloud. ........................................................................ 97
Figure 4.10: Scatter plots showing the variance of the velocity space for the two models. .... 98
Figure 4.11: Plots showing the ratio of the volume of an ellipsoid to the volume of a sphere
whose radius is the average value of the length of the three ellipsoid axes. ............... 100
Figure 4.12: The dust number density distribution. ................................................................ 102
Figure 4.13: Schematic view illustrating the orientation used for figures of scalar field
quantities. ................................................................................................................... 103
Figure 4.14: Spatial distribution of the drift velocity vectors. ................................................ 105
xi
Figure 4.15: Spatial distribution of the drift velocity vector magnitude. ................................ 106
Figure 4.16: Spatial distribution of the velocity space standard deviation values along the ??
?
direction. ................................................................................................................... 107
Figure 4.17: Spatial distribution of the velocity space standard deviation values along the ??
?
direction. ................................................................................................................... 108
Figure 4.18: Spatial distribution of the velocity space standard deviation values along the ??
?
direction. ................................................................................................................... 109
Figure 4.19: Spatial distribution of the velocity space standard deviation values along the ??
1
direction. ................................................................................................................... 110
Figure 4.20: Spatial distribution of the velocity space standard deviation values along the ??
2
direction. ................................................................................................................... 111
Figure 4.21: Spatial distribution of the velocity space standard deviation values along the ??
3
direction. ................................................................................................................... 112
Figure 4.22: Views of the ??
1
vector directions throughout the cloud volume. ....................... 113
Figure 4.23: Views of the ??
2
vector directions throughout the cloud volume. ...................... 114
Figure 4.24: Views of the ??
3
vector directions throughout the cloud volume. ...................... 115
Figure 5.1: The momentum density vector field. ................................................................... 127
Figure 5.2: Particle flux rate through each volume element in the dust cloud structure. ...... 128
Figure 5.3: Drift velocity based energy density. ..................................................................... 130
Figure 5.4: Thermal energy density. ....................................................................................... 132
Figure 5.5: Gravitational potential energy density. ................................................................ 134
Figure 5.6: Comparison of energy density scales. .................................................................. 135
Figure 5.7: Drift force density vector field. ............................................................................ 137
Figure 5.8: Drift force density magnitude. .............................................................................. 138
Figure 5.9: Thermal force density vector field. ...................................................................... 139
Figure 5.10: Thermal force density magnitude. ...................................................................... 140
xii
Figure 5.11: Gravitational force density magnitude. .............................................................. 141
Figure 5.12: Comparison of force density magnitudes. .......................................................... 142
Figure 5.13: Drift velocity based energy flow rate. ................................................................ 144
Figure 5.14: Thermal energy flow rate. .................................................................................. 145
Figure 5.15: Gravitational potential energy flow rate. ............................................................ 146
Figure 5.16: Total measured energy flow rate. ....................................................................... 148
Figure 5.17: Difference in the energy flow rates calculated using the two velocity space
models ......................................................................................................................... 150
Figure 5.18: Comparison of the magnitudes and angles of the heat flux vectors calculated with
the two velocity space models. .................................................................................... 151
1
Chapter 1: Introduction
1.1: Historical overview
The plasma state of matter has been actively investigated in the laboratory for over 130
years. First described by William Crookes in 1879
1
, this fourth ?radiant? state of matter came
about through two mechanisms: The first being a lack of collisions between the constituent
molecules, due to rarefaction of the gas, and second by applying an external force to a gaseous
system and ?coercing them into a methodical rectilinear movement
2
.? (pp. 4712) Crookes?
preferred method of coercion was the application of an electrical bias to various electrodes
within glass tubes from which the air had been evacuated. This is the earliest example of a dc
glow discharge plasma produced in the laboratory; the experiments described in this dissertation
are performed in a device that is quite similar to those described by Crookes.
Interest in the study of plasmas was rekindled in the early part of the twentieth century in
the research labs of the General Electric corporation by Irving Langmuir and Lewi Tonks. The
pair characterized and developed techniques for the mathematical description of plasma systems
for some twenty years, much of which is still used today. Of particular interest from this time
period is Langmuir?s observation of the first laboratory based dusty plasma system, the ?streamer
discharge
3
?, which describes a plasma that contains ?globules? of sputtered tungsten that charge
negatively and interact with the ambient environment. Aside from this early observation, the
thrust of plasma physics research remained mostly in the dustfree regime until the observation
2
of spokelike structures in Saturn?s rings in the early 1980?s. These observations sparked a
renewed and focused interest in dusty plasmas that has continued to this day, with a multitude of
research groups investigating such systems across the globe.
Modern study of dusty plasmas has expanded beyond simple observation of floating
particulate matter; dusty plasma systems have been shown to exhibit, for example, wave
phenomena, order structure formation, and allow detailed study of the thermodynamic properties
of systems in the gaseous, liquid and solid states. The wide range of properties observed within
dusty plasmas makes the study of these systems a robust field of study at the intersection of
chemistry, materials science, fluid mechanics, and fluid/plasma physics. Section 1.2 will give a
brief description of some very basic parameters that can be used to characterize these and other
plasma systems and explain how dust grains acquire charge in a plasma environment. The
motivation for, and scope of, this dissertation is outlined in Section 1.3.
1.2: Dusty plasma parameters
Dusty (?complex?) plasmas are systems composed of electrons, ions, neutral gas, and
charged particulate matter. The inclusion of the fourth plasma species (the ?dust?) is what
separates dusty plasma physics from the standard picture of a plasma system. All plasma
systems can be characterized by a fairly standard set of parameters, including: Densities,
temperatures, screening lengths, and various characteristic frequencies. The inclusion of a dust
component in the derivation of these parameters can have a large effect of the values of these
quantities. Section 1.2.1 contains an overview of how some of these parameters change when
the dust is included and introduces an important dimensionless quantity used to characterize the
dust component. Section 1.2.2 discusses the subject of how dust grains accumulate electrical
3
charge. As will be shown in later chapters, many of the assumptions that are fundamental
components in the derivation of these basic parameters are, at best, only partially valid; but the
discussion of these quantities is quite useful because it facilitates the comparison of dusty
plasmas with other, more standard, plasmas.
1.2.1: Plasma parameters
The introduction of dust grains into a plasma environment changes many of the basic
plasma physics parameters that are used to characterize such systems
4
. In this section several
simplifying assumptions are made about the plasma: The ions are assumed to be singly charged
and positive, the electrons, ions, and dust are assumed to have isotropic Maxwellian velocity
space distributions with no drift, and the dust grains are assumed to have a negative charge. In
what follows, the subscripts d, e, and i denote the dust grain, electron, and ion components of the
plasma, respectively.
1.2.1.1: Quasineutrality
In a plasma environment the system is assumed to be macroscopically quasineutral, this
condition can be expressed as:
where the subscript s refers to any of the plasma species. The symbol e (nonsubscript) is the
elementary (proton) charge, n
s
is the number density of plasma species s, and q
d
is the charge of
a dust grain. The dust grain charge can be rewritten as
dd
Zeq ?= , where Z
d
is the dust grain
charge number. The dust grain charge number is defined to be greater than zero which gives the
ddeiss
n qn en eq n0 +==
?
s
4
minus sign in the equation for q
d
. With this definition of the dust charge the quasineutrality
becomes:
ddie
ddei
n Znn
)n Zn(n e0
=
=
1.1
In some dusty plasma environments a large fraction of the electrons (and to a lesser extent, the
ions) can reside on the surface of the dust grains, an effect known as electron (ion) depletion.
This effect will be discussed further in Section 1.2.2.
1.2.1.2: Debye length
A fundamental property of a plasma is the ability to shield external electric fields. This
Debye shielding is described in terms of a characteristic length scale over which the shielding
occurs. For electrons and ions the Debye length is given by:
?
?,(?,?)
= ?
?
0
?
?
?
(?,?)
?
(?,?)
?
2
where ?
0
is the permittivity of free space, k
B
is Boltzmann's constant, and T
(e,i)
is the electron or
ion kinetic temperature. This characteristic shielding length in a plasma where dust grains are
present can be found as follows: Assume a plasma with an initially uniform background of
electrons and ions. If a single positively charged immobile sphere is placed in the plasma
electrons will be attracted to the sphere and will form a "cloud" around the sphere that acts to
cancel the local electric field introduced by the positive charge of the sphere. The dust Debye
length, ?
D,d
, is a measure of the thickness of this cloud of ions and electrons. The calculation
begins Gauss' law:
5
( )
ddei
00
2
n Znn
??
? e
c
?
=
?
=??
1.2
where
c
? is the electrostatic potential near the dust grain. Next, the quasineutrality relation,
Equation 1.1, is used to eliminate the dust term in Gauss? law, the electron and ion distributions
are assumed to have uniform spatial distributions far away from the grain:
?
2
?
?
=
??
?
0
??
?0
????
???
?
?
?
?
?
???
?0
????
??
?
?
?
?
?
?? (?
?0
??
?0
)?
=
? ?
?0
?
0
?1 ?????
???
?
?
?
?
?
???
? ?
?0
?
0
?1 ?????
??
?
?
?
?
?
??
Next, the standard assumption
?
(?,?)
?
?
?
?
?
(?,?)
? 1 is made and the terms found in the brackets are
expanded in a Taylor series to first order in the small quantity.
c
eDiD
c
eB
e
iB
i
c
Tk
ne
Tk
ne
?
?
?
?
?
?
?
?
?
+=?
?
?
?
?
?
?
?
?
+=??
2
,
2
,
0
,0
2
0
,0
2
2
11
??
??
1.3
If we then assume that the potential from the dust grain falls off as:
?
?
(?) = ?
?0
?????? ?
?,?
? ?
where ?
?0
is the potential at the grain surface (at ? = 0) and ?
?,?
is the dust Debye length, it is
easy to show that ?
?,?
is given by:
1
?
?,?
2
=
1
?
?,?
2
+
1
?
?,?
2
Thus, this characteristic length scale is entirely determined by the ambient electron and ion
environment.
6
1.2.1.3: Plasma frequency
As in the case of an electronion plasma, small displacements of the charged species in a
uniform, low temperature, dusty plasma give rise to stabilizing restoring forces which operate at
rates characterized by a plasma frequency. The calculation of this frequency for a system
containing electrons, ions, and charged dust grains proceeds in the usual manner. Beginning
with the conservation of number density:
0 =
??
?
??
+ ?? (?
?
??
?
)
And the fluid momentum equation (with no pressure gradient or collision terms, discussed in
much greater detail in Chapter 5):
0 =
???
?
??
+ (??
?
??)??
?
+
?
?
?
?
??
Gauss' law with the quasineutrality condition, Equation 1.2:
?
2
? =
?1
?
0
??
?
?
?
?
We now assume no zeroorder drift velocity and small number density and electric potential
fluctuations for all species:
sss
nnn
,1,0
+=
10
??? +=
where
10
?? >> and
ss
nn
,1,0
>> . The above equations can be linearized and combined to give
the following differential equation for ?
2
?
1
:
7
?
2
??
2
?
2
?
1
= ???
?
?
2
?
0,?
?
0
?
?
?
??
2
?
1
which is the differential equation for a simple harmonic oscillator. The plasma frequency of the
system is then:
?
?
2
= ?
?
?
2
?
0,?
?
0
?
?
?
This is the standard expression for the plasma frequency; the only difference, compared to the
usual case of a plasma consisting of electrons and ions, is that the dust component values must be
included when performing the summation.
1.2.1.4: Coulomb coupling parameter
The general thermodynamic state of the dust component can be characterized by the
Coulomb coupling parameter, ?
c
. This parameter is the ratio of a characteristic dust grain?s
electrostatic potential energy to the grain?s thermal kinetic energy. The potential energy used is
the electrostatic potential energy between the characteristic dust grain and a neighboring grain,
this energy is given by:
?=
d
qW
where ? is the Debye?H?ckel potential, which is a solution to Equation 1.3:
?
?
?
?
?
?
?
?
?
=?
dD
r
r
,0
exp
4
1
???
This potential is evaluated at r=a, the average separation distance between dust grains. Using
this definition of the potential energy, the coupling parameter ratio is given by:
8
Figure 1.1: Dust clouds with different phases. (a) A color inverted image of a strongly coupled dusty plasma
5
. (b)
An image of a dust cloud with two phases. The portion of the cloud in the white (upper) box is gaslike and the
portion of the cloud in the red (lower) box is liquidlike.
?
?
?
?
?
?
?
?
?
=?
dDdB
d
c
a
Tak
Ze
,0
22
exp
4 ???
This plasma parameter describes the spatial correlation between neighboring dust grains. In the
limit
c
? >> 1 a grain's electrostatic potential energy is much larger than its thermal kinetic
energy. This type of system is said to be a stronglycoupled dusty plasma, the dust grains form a
Coulomb solid, the socalled dust crystals, as seen in Figure 1.1 (a). When
c
? ? 1 the dust
component of the plasma is said to be moderatelycoupled, in this regime the dust is in a liquid
like state (Figure 1.1 (b), in the lower box outlined in red). Finally, when
c
? << 1 the system is
said to be weaklycoupled and behaves as a Coulomb gas (Figure 1.1 (b), in the upper box
outlined in white). The experiments described in this work fall into the latter, weaklycoupled,
regime.
The derivation of the plasma parameters presented in this section does not vary greatly
from the process used in standard plasma systems, with the exception of the introduction of the
Coulomb coupling parameter.
9
1.2.2: Dust grain charging
There are a number of processes that can affect the charge of a dust grain found within a
plasma environment, these processes fall into two basic categories: The interaction of the dust
grains with the electron and ion populations in the plasma and the interaction of dust grains with
photons. An approximate value for the dust grain charge can be calculated by considering the
current to (and from) the dust grain surface due to the various processes:
?
=
k
k
d
I
dt
dq
1.4
where the sum over k refers to each charging mechanism. In order to calculate the dust grain
charge within a given experiment one must identify the charging processes that are important in
the ambient environment and express each process as a current, I
k
, to (or from) a dust grain. In
this section a brief description will be given for several common charging processes, after which
an equation will be derived for the dust charge based on the charging currents that are important
in the experiment described in chapters three and four.
1.2.2.1: Charging processes overview
When an initially neutral dust grain is introduced into a plasma environment, consisting
of electrons and positively charged ions, it will accumulate a net charge on its surface as a result
of low energy collisions with the ambient electron and ion populations. Because electrons have a
much lower mass than the ions, the electron thermal speed (
ssBsth
mTkv /2
,
= ) is much higher
than the ion thermal speed, this causes a dust grain to collide with the electrons more frequently
and results in a dust grain charge that is net negative. The process is brought to equilibrium
through Coulombic repulsion. After the grain acquires a negative charge the rate of electron
10
collection is decreased due to the resulting dustelectron electrostatic repulsion, at the same time
the dustion collection rate is enhanced due to Coulombic attraction. These electron and ion
collection currents (I
c,e
and I
c,i
) are the dominant charging processes in the experiment described
later.
The second type of charging process that arises due to the interaction of a dust grain with
the background plasma is secondary electron emission. This is the process by which high energy
electrons or ions collide with a dust grain with sufficient energy to cause the ejection of electrons
from the dust grain surface. When a low energy electron or ion collides with a dust grain it
sticks to the surface of the grain (this is the process of collection, described above); if the impact
energy is high enough the incident electron or ion will tunnel into the grain interior. As the
particle tunnels some of its energy is transferred to the dust grain and some of its energy excites
electrons on the surface which allows electrons to escape from the grain. The collision energy
threshold for this process to occur is in the 1 keV range, which is much higher than the electron
and ion energies found in dc glow discharge plasmas.
The third broad class of charging processes is the interaction of the dust component with
photons, via the photoelectric effect. If the energy of a photon incident to a dust grain is larger
than the grain's work function photoelectrons are emitted from the grain. This charging process
can be quite important for dusty plasma systems found in space, where there can be a high flux
of UV radiation.
It is noted that there are additional dust grain charging mechanisms (including field
emission, radioactive emission, neutral impact ionization, and thermionic emission) but these
11
processes require relatively extreme background plasma parameters or environments in order to
become important charging mechanisms.
1.2.2.2: Dust charging due to electron and ion collection
In this section a transcendental equation relating the charge of a dust grain to the
parameters of the ambient plasma environment and the properties of the dust grain is found. The
most straightforward method for calculating the ion and electron collection currents to the
surface of a dust grain is the Orbit Motion Limited (OML) method
4,6
. In this charging process
model several assumptions are made, namely: The dust grains are spherical conductors of radius
r
d
, the dust charge is steady state (and negative), and the velocity distribution functions of the
background electrons and ions are assumed to be isotropic and Maxwellian with no drift
velocity:
( )
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
=
sB
ss
sB
s
sss
Tk
vm
Tk
m
nvf
2
exp
2
2
2/3
,0
?
where v
s
is the velocity of species s far from the dust grain, before the collision.
If a low energy electron or ion collides with a dust grain it will stick to the grain's surface,
for this to occur the incident charged particle must have an impact parameter less than the impact
parameter of a grazing collision, b
s
(a schematic of a grazing collision can be found in Figure
1.2). The value of this critical impact parameter will vary with the dust grain charge, due to
Coulombic attraction (ions) or repulsion (electrons), and the velocity of the incident charged
particle v
s
. The impact parameter can be found by considering the conservation of momentum
and energy:
12
Figure 1.2: Schematic of a grazing collision between an electron or ion and a dust grain.
sfsssss
bvmbvm
,
=
1.5
d
ds
fssss
r
qq
vmvm
0
2
,
2
42
1
2
1
??
+=
1.6
The velocity of the particle after the grazing collision, v
s,f
, can be eliminated from Equations 1.5
and 1.6, which gives an expression for the critical impact parameter:
2
0
2
1
ssd
ds
ds
vmr
eZq
rb
??
+=
Next, the current to the dust grain surface can be found using:
( )
?
=
V
sss
d
ssss
vdvfvqI
3
?
where
2
s
d
s
b?? = is the dusts collision cross section and the region of integration, V, is over the
volume of velocity space from which particles can collide with the dust grain. Due to the
assumption of velocity space isotropy this can be reexpressed as an integral over one velocity
space dimension as:
13
( )
?
?
=
min
232
4
s
v
sssssss
dvvfbvqI ?
where the parameter ?
?
???
is the minimum velocity of a particle s that is required to overcome
the dust grain  s electrostatic repulsion and result in a collision. In this case the dust charge is
assumed to be negative, meaning the minimum ion velocity value, ?
?
???
, is zero because the
electrostatic force is attractive. However, the dust grain  electron electrostatic interaction is
repulsive, resulting in a nonzero ?
?
???
which can be found using the conservation of energy,
Equation 1.6:
ed
d
e
mr
Ze
v
0
2
min
2??
=
with this information the ion and electron collection currents are given by:
i
Tkvm
iid
d
i
iB
i
id
dve
vmr
Ze
v
Tk
m
ner
iBii
2/
0
2
0
2
3
2/3
,0
22
ic,
2
2
1
2
4I
?
?
? ?
?
?
?
?
?
?
?
+
?
?
?
?
?
?
?
?
=
???
?
?
?
?
?
?
?
?
?
+=
iBd
d
i
iB
idic
Tkr
Ze
m
Tk
nerI
0
2
,0
2
,
2
1
2
2
??
?
1.7
And
e
Tkvm
v
eed
d
e
eB
e
ed
dve
vmr
Ze
v
Tk
m
ner
eBee
e
2/
2
0
2
3
2/3
,0
22
ec,
2
min
2
1
2
4I
?
?
? ?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?=
???
?
?
?
?
?
?
?
?
? ?
?=
eBd
d
e
eB
edec
Tkr
Ze
m
Tk
nerI
0
2
,0
2
,
4
exp
2
2
??
?
1.8
For the plasma environment found in in the experiment described in what follows the electron
and ion collection currents are the dominant charging mechanisms.
14
Substitution of Equations 1.7 and 1.8 into Equation 1.4, and by assuming a steady state
dust grain charge, an equation can be obtained for the dust grain charge number, Z
d
:
ecic
d
II
dt
dq
,,
0 +==
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
+=
eBd
d
iBd
d
ie
ei
i
e
Tkr
Ze
Tkr
Ze
mT
mT
n
n
0
2
0
2
,0
,0
4
exp
4
1
????
1.9
This equation can be solved numerically for two scenarios, which arise from algebraic
manipulation of the quasineutrality condition, Equation 1.1:
i
dd
i
e
n
nZ
n
n
,0,0
,0
1?=
The first case is the socalled "dust in plasma" limit, in this limit
idd
nnZ
,0
<< and represents the
case where the dust grains are isolated from one another in the plasma. In this limit there is no
electron depletion and the dust charge, Equation 1.9, is given by:
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
+=
eBd
d
iBd
d
ie
ei
Tkr
Ze
Tkr
Ze
mT
mT
0
2
0
2
4
exp
4
11
????
The second scenario comes when
idd
nnZ
,0
? , this is the socalled "dusty plasma" limit.
In this limit a nonzero fraction of the plasma?s electron population resides on the surface of dust
grains (and the ambient electron population is said to be ?depleted?). The dust grain charge,
Equation 1.9, in this limit is given by the solution to:
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
+=?
eBd
d
iBd
d
ie
ei
i
dd
Tkr
Ze
Tkr
Ze
mT
mT
n
nZ
0
2
0
2
,0
4
exp
4
11
????
15
Figure 1.3: Log plot of the dust grain charge number Z
d
, in electrons, vs the ratio of dust number density to ion
number density. The charge is calculated for r
d
= 1.5x10
6
m, k
B
T
i
= 0.025 eV, and k
B
T
e
= 2.5 eV.
A system with a larger dust number density, n
d
, contains dust grains with lower individual charge
(Figure 1.3) but more of the plasma system's total charge is associated with the dust component
of the plasma.
The OML model for the dust grain charge gives an upper limit for the dust charge
number, Z
d
. All of the charging processes mentioned in Section 1.2.2.1 that were not included in
the calculation of the charge equations decrease the number of electrons on the dust grain surface
(i.e. these processes make q
d
less negative and Z
d
less positive). It is noted that more
sophisticated charging models exist that incorporate the neglected processes and take into
account effects such as electron, ion, and dust grain drift velocities, nonspherical dust grains,
and the effect of the sheath potential structure on the collision process
4,6,7
. These mechanisms
are typically neglected due to the vast uncertainty in the plasma parameters that are required for
even the relatively simple OML description.
0.0 0.2 0.4 0.6 0.8 1.0
1
10
100
1000
Dust number density to ion number density ratio :
n
d
n
i,0
Z
d
:
el
ect
r
o
n
s
Dust Grain Charge vs Dust to Ion Number Density Ratio
16
1.3: Motivation and scope
The motivation for this work can be briefly summarized by noting that, in recent years,
many of the assumptions that went into the derivation of the parameters discussed in Section 1.2
have been found to be overly simplistic. Before any type of rigorous, systematic, theoretical
knowledge of weaklycoupled dusty plasmas can mature to the point where the majority of
observed phenomena can be readily explained, the basic properties of such systems must be
elucidated within the environments in which they are observed. Weaklycoupled dusty plasmas,
particularly those found in dc glow discharge devices, have a number of important basic
properties that require investigation. Many of these properties affect even the most basic
descriptions of such systems; for example, recent work has shown that the distribution of the dust
component of the plasma within velocity space is anisotropic
8?12
and varies by position.
13?15
Additionally, these measurements of the velocity space have resulted in dust temperature
values
16,17
that can reach hundreds, or even thousands, of electron volts (in energy units) which is
orders of magnitude larger than the temperature of any other plasma species found in the
experiments. The experimental results discussed here show that the observed velocity space
anisotropy is both a real and a large effect. The measurements of the velocity space will be
modeled with a distribution that allows the previously unexplained anisotropy to be accounted
for in a physically meaningful way. The results will also show that the variation of various
properties of the dust component can be measured in weaklycoupled dusty plasmas with
unprecedented detail.
The measurements of the dust component properties, in conjunction with the model that
allows the anisotropy, leads to a discussion of the measured transport and distributions of
particles, momentum, energy, forces, heat fluxes, and energy flow rates; all quantities that had
17
been, previously, experimentally unavailable in such a system. The derivation of the transport
equation relations will give more general expressions for several commonly used fluid quantities
than have been previously used to describe plasma fluid systems.
The presentation of this dissertation is as follows. Chapter two contains discussion of the
basic theoretical description of plasmas as fluids in phase space, the new velocity space
distribution function that is used, and a brief introduction to Bayesian probability theory.
Chapter three describes the experimental apparatus in which the experiment was performed and
contains a discussion of the diagnostic used to measure the six dimensional distribution of the
dust component of the plasma within phase space. Chapter four details the process by which the
measurements can be used to construct the spatially resolved phase space distribution and shows
that the anisotropic model of velocity space is required for an adequate characterization of the
dust component. Various thermal and transport properties of the system are discussed in detail
within Chapter five and the various conclusions that are drawn during the course of the
dissertation are summarized in Chapter six.
18
Chapter 2: Theory and Background
In this chapter the basic theoretical constructs required for the discussion of the
experiment, analysis, and results will be reviewed. The chapter begins by describing the phase
space used to describe the observed plasma system. Section 2.2 describes the two probability
distribution functions used to model the velocity space portion of phase space. Bayesian
probability theory is reviewed in Section 2.3; emphasis is given to the processes used to obtain
optimal model parameter values given a data set and the methodology that allows comparison of
theoretical models that will be used extensively in Chapter 4.
2.1: The phase space distribution
The investigation of systems consisting of many small particles has a long history with a
myriad of different specific approaches and methodologies. As the basis for many of these
systems of description lies the idea of a phase space distribution (PSD), which simultaneously
encapsulates information about both the configuration space and velocity (or momentum) space
distributions of the system in question. For the plasma system described below the configuration
space and velocity space both consist of three dimensions; when combined with time, the phase
space of the plasma system consists of seven dimensions. The N individual particles that make
up each of the plasma components (the ions, electrons, neutral gas, and dust) each have a
trajectory through the seven dimensional phase space which, at least in principal, must be known
for a complete deterministic description of the dynamics for each component. The position and
19
velocity of particle j, at time t, is given by the curves ??
?
(?) = {??
?
(?),??
?
(?)}, but due to the large
number of particles, and the extremely complicated nature of the individual particle trajectories,
knowledge of all of the individual ??
?
(?) curves is difficult to measure (except in the case of
carefully prepared simple systems). Even if known, the N six dimensional particle trajectories
would be difficult to use as the basis for an efficient, practical, description of the plasma
component in question. In an effort to provide a tractable method for classifying and extracting
useful information about such a system a simplified, statistics based, view is employed by
treating the plasma components as fluids.
In the phase space fluid picture, knowledge of the trajectories of each particle is replaced
by a statistical description in which a phase space density function, ? (??,??,?), gives the
probability that there are some number of particles within the phase space volume element
?
3
??
3
? at time t. The fluid picture reduces the information required from a six dimensional
vector field for each particle in the plasma component to a single scalar field for all particles
within the system. The fluid description of the dust component of the plasma considered here is
further simplified by the fact that the system is in a steady state (for time scales on the order of
hours or longer) which eliminates the time dimension. The phase space distribution (PSD) is
then:
? (??,??,?) = ? (??)? (??,??) 2.1
where ? (??) is the particle number density within the configuration space volume element ??
3
and ?(??,??) is the function that describes the distribution of velocities within the volume
element. With the fluid picture, several useful quantities are easily obtained from the known
20
PSD by integrating over various components of phase space; for example, the number density
scalar field, ? (??), is found by integrating the PSD over the velocity space:
? (??) = ? ?(??,??,?)?
3
?.
?
??
2.2
Similarly, the total number of particles in the system can be found by integrating over all of
phase space (or in a smaller region, by integrating over the configuration space volume of
interest):
? = ? ? (??) ?
3
?
?
??
= ? ?(??,??,?)?
3
??
3
?
?
??
2.3
The plasma as a fluid model gives direct access to many other physical quantities of interest,
which can be obtained by taking velocity space moments of the PSD. Perhaps the most useful
examples of such quantities are given by the first two velocity space moments of the PSD. The
first velocity space moment, normalized to the number density, gives the drift (average) velocity
vector field of the plasma fluid, ??? (??):
??? (??) =
1
? (??)
? ?? ?(??,??,?)?
3
?.
?
??
2.4
The pressure tensor field is proportional to the second central moment of the PSD:
?
?
(??) =
1
2
?? (?? ????) ? (?? ????)?(??,??,?)?
3
?.
?
??
2.5
where "?" indicates the tensor (outer) product and m is the mass of an individual particle of the
plasma species in question. Similarly, the second "noncentral" moment of the PSD gives the
total energy density tensor, which is the sum of the pressure tensor and the kinetic energy
associated with the drift velocity:
21
?
?
(??) +
1
2
??(??)????(??) ????(??)? =
1
2
?? ?? ????(??,??,?)?
3
?.
?
??
2.6
The process of taking moments of the velocity space proceeds ad infinitum. The fluid picture is
quite powerful and relatively straight forward, in principal, if the number density field and
velocity distribution functions are known.
The basic approach towards applying the fluid model of plasmas is normally quite
different depending if one studies the system in a textbook or experimental manner. As a vast
oversimplification of the processes: The textbook approach starts by specifying the PSD and
finds the resulting fluid parameters while the experimental approach is to measure the fluid
parameters and work backwards towards the PSD. The experimentalist is normally resigned to
such an approach because the measurement tools at his (or her) disposal cannot generally be used
to directly measure the PSD. For example, Langmuir probes of various types can be used to
measure the current drawn to an electrode immersed in a plasma environment as a function of the
voltage applied to the electrode. From the measured relationship between the current and
voltage a large number of assumptions are made regarding the likely properties of the PSD, with
these assumptions quantities such as the average kinetic energy can be extracted for the
electrons. Then, based on the average kinetic energy measurements, one can specify quantities
such as the width of the velocity space distribution (the thermal velocity) or the number density,
which are found in the PSD. The experimental process is filled with assumptions, the validity of
which are very difficult to test.
The dusty plasma system, as described in the following chapters, has several properties
that allow the experimentalist to avoid such problems and approach the system from a
perspective that is closer to that of the textbook. The large radii of the particles that make up the
22
dust component of the plasma allow the number density scalar field to be measured directly by
taking and analyzing photographs of the particles. The large size of the particles also means that
the mass of the individual grains is relatively large, the result is that the dynamics of the dust
system are slow, which is conducive to the direct measurement of the velocities for particles in
the dust component, many avenues exist for the measurement of the dust velocity depending on
the specific characteristics of the system.
8,18?24
The ability to measure the number density and
the velocities of the particles means that the components of the PSD can be measured directly.
In this dissertation a method is discussed that allows the full PSD of the dust component to be
measured within many small volume elements of configuration space within a dust cloud
structure. This type of measurement gives the spatially resolved PSD of the dust fluid and
allows one to proceed through an analysis hierarchy that starts by measuring and modeling the
PSD and working towards using these measured quantities to calculate the various fluid
properties of the system.
2.2: Velocity space distribution models
The phase space of the dust component of a plasma is composed of both configuration
space and velocity space; the configuration space portion gives information regarding the
number density of the particles within a given volume element of configuration space (which
will simply referred to as a "volume element" or, less commonly, as an "element" in what
follows). The number density will be obtained through analysis of digital photographs of the
dust cloud in a process described in detail within Chapter 4. The end result of such analysis is a
scalar value for the number density, which is simply the average number of particles within a
volume element and an associated uncertainty which is used in error analysis. The velocity
space component contains much more information and is modeled by a probability distribution
23
function within each element. Towards this end, both the standard Maxwellian distribution
function
25
and a generalization of the Maxwellian distribution function, the trinormal probability
distribution function, are used to model the measurements of the velocity space. A brief review
of the standard Maxwellian distribution function is given in Section 2.2.1. The less commonly
encountered multinormal distribution function is described in Section 2.2.2 and a convenient
coordinate system rotation that simplifies the multinormal distribution is highlighted in Section
2.2.3. Section 2.2.4 shows the relationship between the distribution functions used to model the
velocity space in this work and many of the other probability distribution functions commonly
used in the sciences within a convenient hierarchy. Finally, in Section 2.2.5 the general
applicability of the multinormal distribution function is discussed and several interesting
historical notes are highlighted.
2.2.1: The Maxwellian distribution function
The standard probability distribution function used to model the velocity space
component of phase space within plasmas is the three dimensional Maxwellian distribution
function. This probability distribution is ubiquitous across all fields in the sciences but will be
discussed in detail to facilitate the discussion of the more mathematically complicated model that
follows. The ddimensional Maxwellian probability distribution function has the form:
?
?
(??;?) =
1
(2??
?
2
)
? 2?
exp ?
?(?? ????)
2
2?
?
2
? 2.7
where ?? is the ddimensional velocity vector, ??? is the ddimensional drift (mean) velocity vector,
and ?
?
2
is the variance of the probability distribution (a scalar quantity). The square root of the
variance (the standard deviation) of the velocity space is the "thermal velocity" which relates the
24
random fluctuations in velocity space to a temperature, T, and the mass, m, of a single particle of
the plasma species in question:
?
?
= ??
?
? ??
2.8
where k
B
is the Boltzmann constant.
The one, two, and three dimensional versions of Equation 2.7 are of particular interest, all
three of these distributions will be referred to as the "Maxwellian" distribution function without
reference to the dimensionality, except in cases where such a reference is convenient in order to
eliminate ambiguity. The functional form of the three distributions is given below, in order of
increasing dimensionality, in the general ?? = ??
?
,?
?
,?
?
? coordinate system:
?
?
(?
?
) = (2??
?
2
)
?1 2?
exp?
?(?
?
??
?
)
2
2?
?
2
? 2.9
?
?
(?
?
,?
?
) = (2??
?
)
?1
exp ?
?1
2
?
?(?
?
??
?
)
2
?
?
2
+
?(?
?
??
?
)
2
?
?
2
?? 2.10
?
?
??
?
,?
?
,?
?
? = (2??
?
2
)
?3 2?
?
exp ?
?1
2
?
?(?
?
??
?
)
2
?
?
2
+
?(?
?
??
?
)
2
?
?
2
+
?(?
?
??
?
)
2
?
?
2
??
2.11
The most important aspect of this distribution function, for this discussion, is that for the two and
three dimensional versions of the Maxwellian distribution the variance is identical in every
vector direction, regardless of the choice of coordinate system orientation. This is the
manifestation of scalar variance of the Maxwellian probability distribution and results in a
circular symmetry for the two dimensional case and a spherical symmetry for the three
dimensional distribution.
25
Figure 2.1: Plots of the one, two, and three dimensional Maxwellian distribution function. In all cases the standard
deviation is two and the offset (drift) is set to zero. (a) The one dimensional Maxwellian distribution. (b) The two
dimensional Maxwellian distribution. The amplitude is colored according to magnitude. The black band around the
peak indicates the contour of constant probability amplitude at the standard deviation. (c) Two dimensional view of
the black band seen in (b). (d) The surface of constant probability amplitude for the three dimensional Maxwellian
distribution is shown as gold. The green planes indicate the coordinate axes. The black circles show the intersection
of the standard deviation surface and the coordinate axis planes.
A geometric picture of the Maxwellian distribution is a useful way to think about the
function and can be built up starting with the onedimensional case, as in Equation 2.9. Figure
2.1 (a) shows the one dimensional Maxwellian distribution function plotted with the parameter
choices u=0 and ?
?
= 2. The circular symmetry of the two dimensional Maxwellian can be
seen in Figures 2.1 (b) and 2.1 (c). Figure 2.1 (b) shows the two dimensional Maxwellian plotted
with the parameters ??? = 0 ??
?
+ 0 ??
?
and ?
?
= 2, the black band around the peak of the
26
distribution is what will be referred to as the "standard deviation contour" (a manifold of
constant probability amplitude). This same contour of constant probability amplitude is shown
in Figure 2.1 (c) as a two dimensional plot, where circular symmetry of the distribution is
immediately apparent. Moving to the case of the three dimensional Maxwellian distribution, as
in Equation 2.11, a plot of the probability amplitude is shown in Figure 2.1 (d), the plot shows
the three dimensional surface of constant probability amplitude for the parameter choices
??? = 0 ??
?
+ 0 ??
?
+ 0 ??
?
and ?
?
= 2. The black circles that lie on the translucent green planes
indicate the intersection of the spherical surface and the three planes formed by the coordinate
axes. The observation that the constant probability surface is spherical is equivalent to the
statement that the distribution in velocity space described by the Maxwellian is spherically
symmetric, or that the velocity space distribution is isotropic.
2.2.2: The multinormal distribution function
A multidimensional generalization of the Maxwellian probability distribution function is
the multinormal probability distribution function. This distribution is also known as the
Maxwellian distribution function with tensor variance. The ddimensional multinormal
distribution function has the form:
? (??;?) =
1
(2?)
? 2?
??
?
?
1 2?
exp?
?1
2
(?? ????)
?
? ?
?
?1
? (?? ????)?
2.12
where ?
?
is the symmetric ? ? ? covariance tensor for the velocity space. The covariance tensors
for the two and three dimensional cases of Equation 2.12, in the general ?? =
??
?
,?
?
,?
?
? coordinate system, are given by:
27
?
?
= ?
?
?
?
2
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
2
? 2.13
?
?
= ?
?
?
?
2
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
2
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
2
?. 2.14
The ?
?
?
terms are the standard deviation values along the ??
?
axis and the ?
?
?
?
?
terms represent the
statistical correlation of the velocity space in the {?
?
,?
?
} subspace. The correlation factors are
dimensionless and have values in the range ?1 < ? < 1. In an attempt to improve clarity, the
two dimensional version of Equation 2.12 will be referred to as the "binormal" distribution
function and the three dimensional case will be referred to as the "trinormal" distribution
function. This naming scheme is somewhat arbitrary due to the fact that Equation 2.12 is also
known as the Maxwellian distribution function with tensor variance. The distributions discussed
in Section 2.2.1 will be referred to as "Maxwellian" and those found in this section as the "d
normal" distributions.
The geometric picture of the binormal and trinormal distribution functions is an
extremely convenient way to understand how the inclusion of tensor variance affects the
distributions. We begin with the binormal distribution. Figure 2.2 (a) shows the probability
amplitude for the binormal distribution with the parameters ?
?
?
= ?
?
?
= 2 and ?
?
?
?
?
= 0.5,
which are combined to give the covariance matrix:
28
Figure 2.2: Plots showing the binormal distribution function with the parameters chosen to be: ?
?
?
= ?
?
?
= 2 and
?
?
?
?
?
= 0.5. (a) The binormal distribution, shown with the probability amplitude colored according to magnitude.
The black band indicates the constant probability amplitude contour at the standard deviation. (b) View of the same
standard deviation contour seen in (a) in two dimensions.
?
?
BN
= ?
4 2
2 4
?
2.15
where the subscript "BN" indicates the binormal distribution. The black band around the peak
of the distribution in Figure 2.2 (a) indicates the standard deviation contour, as in Figure 2.1 (b).
Figure 2.2 (b) shows the same contour projected onto two dimensions. The symmetry for the bi
normal distribution is clearly seen to be elliptical, in contrast with the two dimensional
Maxwellian distribution where the symmetry was circular. Other features of note are that the
standard deviation values along both the ??
?
and ??
?
axes is two, as was the case in Figures 2.1 (b)
and (c), but the inclusion of the correlation term has significantly altered the morphology of the
constant probability amplitude manifold. The effect of the sign of the correlation can also been
seen in the figure; for this example the correlation was chosen to be positive, resulting in the bulk
of the ellipse area to be in the two quadrants where ??
?
and ??
?
have the same sign. A positive
(negative) ?
?
?
?
?
value indicates that if the ??
?
component of a velocity vector is greater (less than)
29
Figure 2.3: Plots showing the surface of constant probability amplitude for the trinormal distribution with
parameters: ?
?
?
= ?
?
?
= ?
?
?
= 2, ?
?
?
?
?
= ?
?
?
?
?
= 0.5, and ?
?
?
?
?
= 0.1. Parts (a) and (b) show the same
distribution, but from different vantage points, which allows the ellipsoidal shape of the constant probability surface
to be seen more clearly.
zero there is an increased probability that the ??
?
component of the same velocity vector is also
greater (less than) zero, compared to the case where the correlation is zero.
The standard deviation surface for a trinormal distribution function is shown from two
vantage points in Figure 2.3. The distribution function parameters chosen for the figure are:
?
?
?
= ?
?
?
= ?
?
?
= 2, ?
?
?
?
?
= ?
?
?
?
?
= 0.5, and ?
?
?
?
?
= 0.1, giving the covariance matrix:
?
?
TN
= ?
4 2 0.4
2 4 2
0.4 2 4
? 2.16
where the subscript "TN" indicates the covariance matrix is for the three dimensional (tri
normal) distribution. The black ellipses in Figure 2.3 show the intersection of the ellipsoidal
surface with the planes formed by the coordinate axes, as in Figure 2.1 (d). The red ellipses are
intended to illustrate the ellipsoidal nature of the standard deviation surface. The standard
30
deviation values along the ??
?
,?
?
,?
?
? axial directions are the same as for the Maxwellian
distribution plotted in Figure 2.1 (d), the inclusion of nonzero correlation, through the off
diagonal tensor elements, clearly has a large effect on the resulting distribution.
2.2.3: The principal axis coordinate system
The geometric picture of the distributions discussed in Sections 2.2.1 and 2.2.2 is
convenient because it allows for a more intuitive interpretation of the tensor variance. Without
the geometric picture it can be quite difficult to clearly see the effects of the nonzero correlation
and anisotropy. The full expression for the binormal distribution is given by inserting Equation
2.13 into Equation 2.12:
?
BN
(v??) =
1
2??1 ??
?
?
?
?
2
?
1 2?
?
?
?
?
?
?
exp?
?1
2(1 ??
?
?
?
?
2
)
?
(?
?
??
?
)
2
?
?
?
2
+
(?
?
??
?
)
2
?
?
?
2
?
2?
?
?
?
?
(?
?
??
?
)(?
?
??
?
)
?
?
?
?
?
?
??
2.17
With the help of the geometric picture this can be simplified with a coordinate system rotation,
which is motivated by the elliptical nature of the contour shown in Figure 2.2 (b). We begin with
the binormal distribution found in Figure 2.2 (b) which is also shown in Figure 2.4. The black
arrows in Figures 2.4 indicate the ??
?
and ??
?
coordinate system axes, the red arrows are along the
major axes of the ellipse. Figure 2.4 (a) shows the standard deviation contour in the same
coordinate system as in Figure 2.2 (b). If the coordinate system is rotated, as indicated by the
blue arrow in 2.4 (a), such that the longest ellipse axis is in the ??
1
direction the distribution is
said to be in the "principal axis" coordinate system.
31
Figure 2.4: An example of a rotation to the principal axis coordinate system in two dimensions. The black ellipse in
(a) and (b) is the same standard deviation contour as in Figure 2.2. (a) Shown in the general coordinate system. The
black arrows show the general coordinate system orientation. The red arrows show the principal axis coordinate
system directions. The blue arrow indicates the rotation direction. (b) The same contour and arrows as in (a),
except plotted in the principal axis coordinate system. The major axes of the ellipse can be seen to lie along the
principal axis ordinates.
The principal axis basis set is the coordinate system where the covariance tensor is
diagonal. The rotation matrix that provides the transformation from the general {?
?
,?
?
}
coordinate system into the principal axis coordinate system {?
1
,?
2
} is composed of columns of
the eigenvectors of the covariance matrix computed in the {?
?
,?
?
} system. In the principal axis
basis set the diagonal elements of the covariance tensor are the eigenvalues of the covariance
matrix in the {?
?
,?
?
} system:
?
?
PA
= ?
?
?
1
2
0
0 ?
?
2
2
? 2.18
The functional form of the binormal distribution in the principal axis coordinate system is
somewhat simplified, when compared to the same distribution in the general coordinate system
(as in Equation 2.17):
32
?
PA
(??) =
1
2??
?
1
?
?
2
exp?
?1
2
?
(?
1
??
1
)
2
?
?
1
2
+
(?
2
??
2
)
2
?
?
2
2
?? 2.19
Where the drift velocity vector must also been rotated into the principal axis coordinate system
with the rotation matrix describe above. In this basis set the correlation is zero, meaning the
velocity space components are statistically independent and the two dimensional distribution
above is simply the product of the two one dimensional Maxwellian distributions along the
??
1
and ??
2
axial directions. The simplification provided by rotation to the principal axis
coordinate system is more dramatic in three dimensions. The functional form of the trinormal
distribution function in the general {?
?
,?
?
,?
?
} coordinate system is given below, by combining
Equations 2.14 and 2.12:
?
TN
??
?
? =
1
(2?)
3 2?
?1 ??
?
?
?
?
2
??
?
?
?
?
2
??
?
?
?
?
2
+ 2?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
1 2?
?
?
?
?
?
?
?
?
?
?
exp ?
?1
2(1 ??
?
?
?
?
2
??
?
?
?
?
2
??
?
?
?
?
2
+ 2?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
)
?
(?
?
??
?
)
2
?
?
?
2
+
(?
?
??
?
)
2
?
?
?
2
+
(?
?
??
?
)
2
?
?
?
2
?
2?
?
?
?
?
(?
?
??
?
)(?
?
??
?
)
?
?
?
?
?
?
?
2?
?
?
?
?
(?
?
??
?
)(?
?
??
?
)
?
?
?
?
?
?
?
2?
?
?
?
?
(?
?
??
?
)(?
?
??
?
)
?
?
?
?
?
?
??
2.20
The rotation matrix to the principal axis basis set is, again, given by columns of the covariance
tensor eigenvectors; the actual covariance tensor elements in the principal axis system are the
eigenvalues, as discussed above in the two dimensional case:
?
?
PA
= ?
?
?
1
2
0 0
0 ?
?
2
2
0
0 0 ?
?
3
2
?. 2.21
33
Figure 2.5: An example of a rotation to the principal axis coordinate system for the trinormal distribution function.
Parts (a) and (b) show the standard deviation surface for the same distribution. (a) The black axes are in the general
coordinate system, the black ellipses show the intersection of the surface with the planes formed by the coordinate
system axes. The red arrows indicate the principal axis directions and the red ellipses show the intersection of the
ellipsoidal surface with the planes formed by the principal axis ordinates. (b) The surface is shown in the principal
axis coordinate system.
The rotation is illustrated in Figures 2.5, which shows the same trinormal distribution surface as
Figure 2.3. Figure 2.5 (a) shows the standard deviation ellipsoid in the {?
?
,?
?
,?
?
} coordinate
system (the black arrows) and Figure 2.5 (b) shows the ellipsoid in the principal axis {?
1
,?
2
,?
3
}
system (red arrows). The convention of aligning the longest ellipsoid axis with ??
1
has been
used, as in the two dimensional example. The simplification of the functional form of the
distribution can be seen by combining Equations 2.21 and 2.12:
?
TN
(??) =
1
(2?)
3 2?
?
?
1
?
?
2
?
?
3
exp?
?1
2
?
(?
1
??
1
)
2
?
?
1
2
+
(?
2
??
2
)
2
?
?
2
2
+
(?
3
??
3
)
2
?
?
3
2
?? 2.22
In the principal axis system the velocity space along each direction is statistically independent
from the other vector directions. The three dimensional distribution is simply the product of
three Maxwellian distributions with different variance values.
34
The rotation to the principal axis coordinate system does not result in any loss of
information because the rotation process is defined in terms of tensor rotation; the distributions
given by Equations 2.20 and 2.22 are exactly the same, one is simply described in a more
convenient coordinate system. In addition to retaining the shape of the distribution, there are two
particularly useful quantities that are invariant when the ?
?
tensor is rotated: The trace and
determinant of ?
?
are each preserved with the rotation:
Tr ??
?
? = ?
?
?
2
+ ?
?
?
2
+ ?
?
?
2
= ?
?
1
2
+ ?
?
2
2
+ ?
?
3
2
2.23
??
?
? = (1 ??
?
?
?
?
2
??
?
?
?
?
2
??
?
?
?
?
2
+ 2?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
)?
?
?
2
?
?
?
2
?
?
?
2
= ?
?
1
2
?
?
2
2
?
?
3
2
2.24
The trace will be seen to be related to the thermal kinetic energy density and the square root of
the determinant is proportional to the volume of velocity space occupied by the distribution.
2.2.4: Other generalized distribution functions
The generalization of the Maxwellian distribution to the trinormal distribution is just one
of many options for an alternate distribution if the spherically symmetric Maxwellian is found to
inadequately describe a system. There are myriad distribution functions used in the sciences to
model systems, many of these distributions fall under the Pearson hierarchy of distribution
functions. Introduced in a series of three papers
26?28
between 1895 and 1916, the Pearson
hierarchy of continuous, one dimensional, parameterized distributions stems from limits of the
solution to the following differential equation for the probability density, p(x):
?? (?)
??
=
?(? + (? ??))
?
2
(???)
2
+ ?
1
(???) + ?
0
? (?) 2.25
35
Figure 2.6: The hierarchy of the Pearson family of distribution functions is shown schematically. The path from the
most general distributions to the Maxwellian distribution is outlined in blue. Several probability distributions that
are commonly encountered in physics are shown.
Solutions to this equation fall into two general classes: When (?
1
2
? 4?
2
?
0
) ? 0 the solutions
result in parameterized distributions with support on limited ranges of the real line. When
(?
1
2
? 4?
2
?
0
) < 0 the solutions for p(x) have support for all real positive and negative values of
x. The hierarchy of these distributions is shown schematically in Figure 2.6. To model the
velocity space we require support for all real values of x, such distributions are of Type IV and
36
the resulting reparameterizations and limits. The path from Equation 2.25 to the Maxwellian
distribution is shown in blue in Figure 2.6, several common distributions (and the limits required
to arrive at the distributions) are also shown.
Without going into great detail, the process required to arrive at the Maxwellian
distribution function contains many simplifying assumptions and the Maxwellian can be
considered the ?least general" (or the most restrictive) of any distribution function in the Type IV
hierarchy. The Maxwellian retains two parameters: ?, the mean of the distribution, and ?, the
standard deviation; all other parameters that describe the location and shape of the distribution
are assumed unimportant when this most simple case is used. The mathematical convenience of
the Maxwellian distribution is a direct result of these simplifications.
The process of starting with a one dimensional distribution and generalizing to multiple
dimensions is, as seen in the transition from the Maxwellian to the trinormal distribution, most
easily done when the variance is equal in all directions (i.e. the transition from the one
dimensional to the three dimensional Maxwellian distribution). If the variance (or some other
shape parameter) is different in one or more vector directions then tensor variance must be
incorporated. The variance can be made tensor with a fairly general method developed primarily
by economists to form what are referred to as "elliptical distributions" which seem to have led to
multidimensional analogs for the one dimensional distributions that fall under the Pearson Type
VII heading
29?31
. The purpose of this discussion is to point out that, after the three dimensional
Maxwellian distribution, the trinormal distribution is the next most simple distribution function.
The type of analysis presented in this dissertation for the trinormal distribution could, in
principal, be repeated for the tensorized generalization of any one dimensional distribution. The
37
trinormal is the next step up in terms of complexity and still retains the benefit of being quite
restrictive in terms of parameterization constraints.
2.2.5: Applicability of the trinormal distribution function model
The use of the Maxwellian velocity distribution function to describe gaseous system
dynamics dates back to ClerkMaxwell?s work in the 1860?s.
25
In his 1867 paper ClerkMaxwell
derives the velocity distribution that would take his name from mostly geometrical
considerations. Of particular interest from this seminal work is Equation 27 and the
accompanying text:
"When the gas moves in mass, the velocities now determined are compounded with the
motion of translation of the gas.
When the differential elements of the gas are changing their figure, being
compressed or extended along certain axes, the values of the mean square of the velocity
will be different in different directions. It is probable that the form of the function will
then be
?
1
(?,?,?) =
?
1
????
3 2?
?
?(
?
2
?
2
+
?
2
?
2
+
?
2
?
2
)
. . . . . (27)
where ?, ?, and ? are slightly different. I have not, however, attempted to investigate the
exact distribution of velocities in this case, as the theory of motion of gases does not
require it.
38
When one gas is diffusing through another, or when heat is being conducted
through a gas, the distribution of velocities will be different in the positive and negative
directions, instead of being symmetrical, as in the case we have considered." (p. 64)
Where the Greek letters ?, ?, and ? are the three velocity space coordinates. One can readily
recognize the above equation as the trinormal distribution function, in the principal axis basis
set, Equation 2.22, without drift terms (which he discusses in the first and last quoted sentences).
ClerkMaxwell?s comment about the differential elements changing figure, due to compression
or extension, is also quite applicable to the preceding discussion of the geometrical picture of the
Maxwellian and trinormal distributions. In a follow up paper ClerkMaxwell
32
used the same
basic geometric representation of the distributions as spheres and ellipsoids and showed that a
system is stable only when spherical symmetry is present (in more modern parlance he is
referring to a system in a state of forcefree equilibrium). He goes on to show that a system in a
forcebalanced equilibrium must exhibit extension or compression along the lines of "principal
tension" or "principal stress" which are the same as the principal axis directions discussed in
Section 2.2.3. There seem to be two main reasons that the shape of the velocity space
distribution function was ignored in the ensuing decades. The first is the fact that the initial
development of statistical physics concerned itself with the very specific case of a system in a
static forcefree equilibrium, the formulation of statistical mechanics by ClerkMaxwell, and
more concretely with that of Boltzmann, was in terms of a phase space that consisted of
configuration space and energy space. The mechanics was largely developed with an energy
distribution and not a velocity distribution. The second reason the anisotropy was omitted in the
development of the statistical fluids picture is that the forces acting on most gases do not
generally result in large deviations from the spherical distribution shape (and such deviations are
39
difficult to measure). These considerations may explain why this more general form of the
distribution seems to have been ignored in the ensuing decades.
The idea of generalizing the spherically symmetric distribution was rejoined some nine
decades later, when strong magnetic fields were added to plasma physics experiments. In such
cases the deviation from the spherically symmetric Maxwellian velocity space distribution was
recognized as an important issue. The anisotropy was accounted for within two basic
hierarchies: That of Braginskii
33
and that of Chew, Goldberger, and Low (CGL)
34
. The
approach outlined in the famous review article by Braginskii was to account for the deviation by
allowing a "Maxwellian plus perturbation" velocity space distribution and through the use of
more complicated transport coefficients, which are allowed to vary based on spatial orientation.
The CGL approach was to treat the velocity space distribution as a combination of a Maxwellian
distribution with "nonscalar" (tensor) variance and additive perturbation distributions. The CGL
approach provides for separate pressures parallel and perpendicular to the magnetic field
direction, due to the different collision rates parallel and perpendicular to the magnetic field that
particles encounter in such environments. The mechanism used to arrive at the anisotropic
pressure tensor required Taylor expansion of the Boltzmann equation in the small quantity that is
the ion mass to charge ratio (which, for the dusty plasma system of interest here, is a quantity
that is much larger than unity). As Chew et al. acknowledge, the general case of a Maxwellian
distribution with nonscalar variance is more complicated than in the model they discuss and had
not been previously treated (except, as we have seen, in passing remarks by ClerkMaxwell and
also very briefly by Spitzer
35
). A literature search has indicated that the more general nonscalar
Maxwellian distribution model remains unexplored to this day. (An excellent account of the
earlier development of the normal distribution was given by Stahl
36
.)
40
2.3: Bayesian probability theory
2.3.1: Introduction to Bayesian probability theory
Bayesian probability theory, as used here, is simply a more coherent formulation of the
standard probability theory. The general idea is that by starting with a more complete
mathematical framework of probability theory much of the ambiguity associated with the zoo of
hypothesis tests and associated test statistics can be avoided in the description of a probabilistic
or statistics based system. Bayesian probability theory is based on four logical rules that govern
the mathematics associated with probabilities: A sum rule (Equation 2.26), a product rule
(Equation 2.27), Bayes' theorem (Equation 2.28), and marginalization (Equation 2.29). Much of
the discussion that follows is adapted from Sivia
37
. Beginning with the sum rule:
prob (??) + prob (?
_
?) = 1
2.26
which means that the probability of X occurring given the known information, I, plus the
probability of X not occurring given the same information is unity. The product rule:
prob (?,??) = prob (??)prob (??,?). 2.27
In words, the probability of X and Y given I is equal to the product of the probability of Y given
I and the probability of X given Y and the information. Bayes' theorem:
prob (??,?) =
prob (??,?)prob (??)
prob (??)
2.28
The terms within Bayes' theorem have special names; for example, suppose X is some
hypothesis and Y is a set of data. The term on the left hand side of Equation 2.28, prob(XY,I) =
prob(hypothesisdata,I), is known as the posterior probability. On the right hand side of Equation
2.28 the term prob(YX,I) = prob(datahypothesis, I) is known as the likelihood function of the
41
data given the hypothesis and information, prob(XI) = prob(hypothesisI) is known as the prior
probability function for the hypothesis given I, and the term prob(YI) = prob(dataI) is known as
the evidence (or marginal likelihood) for the data given I. Finally, the marginalization rule:
prob (??) = ? prob (?,??)??
?
??
2.29
The process of marginalization allows for an unknown quantity ("Y" here) to be accounted for in
the probabilistic analysis of a system without explicit knowledge of the value. These four logical
guides act as the coherent framework for analyzing a system. There are many specific
applications of Bayesian based statistics, but for the purposes of this discussion the framework
allows for a very clear process of comparing two different models for a given data set.
2.3.2: Bayesian model selection
The comparison of models used to describe a given data set is a vast subject with a long
history. Very briefly, the standard (nonBayesian) approach involves several steps: Acquisition
of a data set, hypothesizing that a single model describes the data, selecting a statistical test (from
among a rather large number of possibilities), calculating the appropriate test statistic, and
arbitrarily selecting a cutoff value of the test statistic above or below which the hypothesis will
be either accepted or rejected. The process is rife with ambiguity (for example, which statistical
test to use and the cutoff value to use), is open to false positive and false negative conclusions,
and often has difficulty with large data sets (with ?large? generally meaning data sets with ~500
or more data points). The Bayesian model comparison process avoids some of these issues by
recasting the problem as a comparison of the posterior probabilities of two different models
given data and any related information in hand.
42
Given a data set consisting of N measurements (collectively referred to as "D") the
program for comparing two models ("A" and "B") is to calculate the ratio, K, of the posterior
probability for each of the models:
? ?
prob (??,?)
prob (??,?)
2.30
If the posterior probability for model A, given the data, is greater than the posterior probability
for model B, given the same data, then K>1 and model A is selected as the preferred model. To
add slight complication there is a scale for interpreting values of the ratio
38
. It will suffice for
this discussion to note that values of K > 100 are considered "decisive" (i.e. when the posterior
probability of model A given the data is two orders of magnitude higher than the posterior
probability of model B given the data it can be said that the evidence supporting model A is
decisive compared to the evidence for model B). It will be seen that exact choice for the cutoff
value for K is unimportant in this discussion; the K values that will be encountered in the
ensuing chapters are many orders of magnitude larger than unity. Additionally, the framework
provided by the Bayesian approach, through the marginalization rule, provides an unambiguous
method to account for the finite measurement error in the acquisition of data and a clear method
to penalize a model if it has extra parameters. Accounting for the measurement uncertainty will
prove particularly useful in the comparison of the two models for the velocity space because of
the slight anisotropy in the velocimetry diagnostic that was used in the experiments. The
inclusion of the penalty for extra model parameters also proved to be extremely useful; if two
models are similar but one contains extra parameters, the model with extra parameters will
always provide a better "fit" to the data. By inspecting the Maxwellian velocity distribution
function model (Equation 2.11) and the trinormal velocity distribution function model (Equation
2.20) the need for a method to penalize for extra parameters is immediately apparent: The
43
Maxwellian model contains one independent parameter (?
?
) while the trinormal model contains
six independent parameters (?
?
?
, ?
?
?
, ?
?
?
, ?
?
?
?
?
, ?
?
?
?
?
, and ?
?
?
?
?
).
The process for calculating the ratio of posterior probabilities given two models will be
outlined in the remainder of this section in general terms; the process will then be illustrated in
Section 2.3.3 with a relatively simple example. Returning to the N measured quantities and the
models A and B, the calculation of the K ratio begins by applying Bayes' theorem to the
posterior probability for each model in Equation 2.30:
? = ?
prob (??,?)prob (??)
prob (??)
??
prob (??)
prob (??,?)prob (??)
? 2.31
For all applications in this work the data set used in the comparison is the same for both models,
meaning the evidence terms (the prob(DI)) cancel. The remaining terms in 2.31 can be
rearranged to give the product of the ratio of the likelihood for the data given each of the models
and the ratio of the prior probability for each model:
? = ?
prob (??,?)
prob (??,?)
??
prob (??)
prob (??)
?. 2.32
The ratio of prior probabilities (the rightmost ratio in Equation 2.32) requires knowledge of the
applicability of each model before any of the actual data is considered. This is summarized
nicely by Sivia: "As usual, probability theory warns us immediately that [the value of K]
depends partly on what we thought about the two [models] before the analysis of the data. To be
fair, we might take the ratio of the prior terms [...] to be unity; a harsher assignment could be
based on the track records of the theorists!
37
" In all cases considered here the ratio of prior
probabilities is taken to be one, so as to provide a "fair" comparison. By making the choice that
44
the prior probability ratio is unity, the expression for K reduces to the ratio of the likelihood for
the data set given each of the models:
? =
prob (??,?)
prob (??,?)
2.33
The remaining task is to find the expression for the likelihood for the data given each of the
models, each of these expressions will incorporate finite measurement error and penalization for
the parameters contained in the model. This process is described in detail, for a simple model, in
the next section.
2.3.3: A simple example
To illustrate the process of calculating the likelihood for a data set the following example
will be used: The likelihood function for a one dimensional Maxwellian model will be
calculated for a data set consisting of N data points. The measurement of the data points will
have finite uncertainty and the model will be penalized for an unknown standard deviation value.
The one dimensional Maxwellian distribution is given by Equation 2.9, the various subscripts
and the drift term will be omitted in this discussion:
? (?) =
1
?2? ?
exp?
??
2
2?
2
? 2.34
In this section each of the components of the process will be discussed for this simple case in
some detail, to facilitate the discussion of the much more complex distribution comparison that
will appear in Chapter 4.
45
2.3.3.1: Penalizing for unknown distribution parameters
The first task is to assess a penalty for the unknown parameter ?; the process begins with
marginalization of the likelihood function for the combined data set over ?, followed by use of
the product rule:
prob (??,?) = ? prob (??,?,?)prob (??,?)?? 2.35
The last term in the integrand, prob(? A,I), is the prior probability of ?, given the model. This
prior probability function accounts for any a priori knowledge of ?. Here it will be assumed that
? can take any value 0 < ? ? ?
???
with equal likelihood; that is, the prior probability for the
parameter will be assumed uniform:
prob (??,?) = ?
1
?
max
??
min
,?
min
< ? ? ?
max
0 , elsewhere
?
= ?
1
?
max
, 0 < ? ? ?
max
0 , elsewhere
?
2.36
The form of the prior probability found in Equation 2.36 is to ensure proper normalization.
Next, given the data set and model, there will be a value of the parameter ? that gives the best fit,
?
0
, and an associated uncertainty in the parameter value, ??. If the uncertainty in ? is assumed
to be normally distributed about the optimal value it can be incorporated into the likelihood
function via convolution:
prob (??,?,?) = prob (??
0
,?,?)exp ?
?(???
0
)
2
2(??)
2
? 2.37
Inserting Equation 2.36 and 2.37 into the likelihood function found in Equation 2.35 gives:
46
prob (??,?) =
prob (??
0
,?,?)
?
max
? exp?
?(???
0
)
2
2(??)
2
???
?
max
0
2.38
The integral can be carried out to give:
prob (??,?) = ?
?
2
??
?
max
?erf ?
?
0
?2??
?? erf?
?
0
??
max
?2??
??prob (??
0
,?,?) 2.39
Before finding the values for ?
0
and ?? the finite uncertainty in the measurement process must
be included.
2.3.3.2: Finite measurement error
The combined likelihood for the N measurements given model A and the optimal value
of the standard deviation, ?
0
, is the product of the likelihood values for the N individual
measurements:
prob (??
0
,?,?) = ?prob (?
?
?
0
,?,?)
?
?=1
= ?prob (?
?
?,?
0
,?)
?
?=1
2.40
where ? is the magnitude of the known measurement uncertainty and "A" has been absorbed into
"I", for simplicity. In order to incorporate the measurement error, ?, the assumption is made that
the uncertainty is normally distributed, as is standard. This means that if the k
th
measured value
was found to be v
k
the "true" value of the quantity being measured, ??
?
, was taken from a normal
distribution centered at ??
?
with a standard deviation ?:
prob (?
?
??
?
,?,?
0
,?) =
1
?2??
exp ?
?(?
?
???
?
)
2
2?
2
? 2.41
47
The actual value of ??
?
is, of course, unknown; to incorporate the uncertainty we turn to the
likelihood for an individual measurement given the model (the rightmost term in Equation
2.40), marginalize for ??
?
, and apply the product rule:
prob (?
?
?,?
0
,?) = ? prob (?
?
??
?
,?,?
0
,?)prob (??
?
?,?
0
,?)???
?
?
??
2.42
The first term in Equation 2.42 is given by Equation 2.41, this is the term through which the
measurement uncertainty is incorporated into the analysis. The second term in Equation 2.42 is
the likelihood for the "true" measured value, ??
?
, given the model:
prob (??
?
?,?
0
,?) =
1
?2??
0
exp?
???
?
2
2?
0
2
? 2.43
Combining Equations 2.41 and 2.43 with Equation 2.42 gives the likelihood for the value that
was actually measured, v
k
, given the model and measurement uncertainty:
prob (?
?
?,?
0
,?) =
1
2???
0
? exp?
?(?
?
???
?
)
2
2?
2
?exp?
???
?
2
2?
0
2
????
?
?
??
2.44
Integration then gives:
prob (?
?
?,?
0
,?) =
1
?2?(?
2
+ ?
0
2
)
exp?
??
?
2
2(?
2
+ ?
0
2
)
? 2.45
This expression for the likelihood of the k
th
measurement can then be used with the expression
for the likelihood of the entire data set given the model and ?
0
, Equation 2.40:
prob (??
0
,?,?) = ?(2?(?
2
+ ?
0
2
))
?1 2?
exp?
??
?
2
2(?
2
+ ?
0
2
)
?
?
?=1
2.46
48
To collect the progress to this point, Equations 2.46 and 2.39 are combined to give the likelihood
for the entire data set, given model A:
prob(??,?) = ?
?
2
??
?
max
?erf?
?
0
?2??
?? erf ?
?
0
??
max
?2??
?? ?
?(2?(?
2
+ ?
0
2
))
?1 2?
exp?
??
?
2
2(?
2
+ ?
0
2
)
?
?
?=1
2.47
The terms that appear before the product can be thought of as the penalty that is assessed for the
unknown parameter ?.
2.3.3.3: Calculation of the parameter values
The remaining task is to find the optimal parameter value, ?
0
, and the associated
error,??. To find these quantities the likelihood function for the data given ? and the model is
maximized. This likelihood function is obtained by combining Equations 2.37 and 2.46:
prob(??,?,?) =
exp?
?(???
0
)
2
2(??)
2
??2?(?
2
+ ?
0
2
)?
?? 2?
exp?
?1
2(?
2
+ ?
0
2
)
??
?
2
?
?=1
?
2.48
The natural logarithm of this equation is defined to be L:
? ? ln(prob(??,?,?))
=
?(???
0
)
2
2(??)
2
?? ln ([2?(?
2
+ ?
0
2
)]
1 2?
) ?
1
2(?
2
+ ?
0
2
)
??
?
2
?
?=1
2.49
The optimal value ?
0
is found by minimizing L: The partial derivative of L is taken with respect
to ?
0
, set equal to zero, and evaluated for ? = ?
0
:
49
0 = ?
??
??
0
?
?=?
0
?
0
2
=
1
?
??
?
2
?
?=1
??
2
2.50
The first term on the right is the standard expression for the variance. To find the value for ??,
we return to Equation 2.49, for L, and replace the summation with the previous result:
? =
?(???
0
)
2
2(??)
2
?? ln ??2?(?
2
+ ?
0
2
)??
?
2
2.51
The value for ?? is obtained through the second partial derivative of L with respect to??
0
:
0 = ?
?
2
?
??
0
2
?
?=?
0
??
2
=
(?
0
2
+ ?
2
)
2
?(?
0
2
??
2
)
2.52
This result mirrors what one would expect: The uncertainty in the optimal parameter value
decreases when more data are available (as N becomes large). The fact that ?? becomes large in
the case where the measurement error approaches ?
0
also follows expected behavior (in the
extreme case where the square of the measurement error is larger than ?
0
2
the parameter
uncertainty, ??, becomes imaginary, in such a case one may consider performing a different
experiment). Equation 2.52 contains an additional point worth noting, which is somewhat subtle:
As data sets become large the uncertainty in the optimal parameter values decrease, without
regard to how well the optimal parameter value describes the data. As the number of data points
increases the chance of the optimal parameter value being skewed in an important way by outlier
data points decreases. Regardless of how well a given model ultimately describes the data one
can obtain an optimal value for any parameters, and given enough data points, a low value for
50
the uncertainty in the parameter values. If the model does a poor job in its description of a data
set it will become apparent by a relatively low value of the likelihood of the model, given the
information, not by low parameter uncertainty.
2.3.3.4: Summary
To summarize the process of calculating the likelihood for the example: After the data
set and model were defined, the likelihood function was penalized for the unknown parameter ? .
Finite measurement error was then included, via marginalization, and the process for finding the
optimal value of the unknown parameter and the associated uncertainty in the parameter were
discussed. The likelihood for the data set given the model defined in Equation 2.34 is:
prob(??,?) = ?
?
2
??
?
max
?erf?
?
0
?2??
?? erf ?
?
0
??
max
?2??
?? ?
(2?(?
2
+ ?
0
2
))
?? 2?
exp?
?1
2(?
2
+ ?
0
2
)
??
?
2
?
?=1
?
2.53
The optimal parameter value and the associated uncertainty are obtained by minimizing the first
and second derivatives of the logarithm of the likelihood function for the data given ?, A, and I
(Equations 2.49, 2.50, 2.51, and 2.52), respectively:
? ? ln(prob(??,?,?))
= ?
?(? ??
0
)
2
2(??)
2
? ? N ln ??2?(?
2
+ ?
0
2
)??
1
2(?
2
+ ?
0
2
)
??
?
2
?
?=1
0 = ?
??
??
0
?
?=?
0
0 = ?
?
2
?
??
0
2
?
?=?
0
2.54
51
This process gives the numerator of the expression for the posterior probabilities ratio given in
Equation 2.33. The process would be repeated for the same data set given an alternate model, B.
The models used in the actual data analysis are three dimensional and one contains six
parameters, which results in considerably more complicated results than in Equations 2.53 and
2.54, but the process is essentially the same as that used for the simple model discussed here.
52
Chapter 3: Experimental Apparatus
The experiments described in this dissertation were performed on the Three dimensional
Dusty Plasma eXperiment (3DPX) at Auburn University. 3DPX consists of several different
subsystems that work in conjunction with one another to produce and diagnose the dusty plasma
systems of interest. Section 3.1 describes the 3DPX apparatus and contains a short description of
the plasmas produced within 3DPX. A brief introduction to the particle image velocimetry
diagnostic and its application to dusty plasma systems is given in Section 3.2.
3.1: 3DPX
The 3DPX device is the combination of several different subsystems that work to
provide a repeatable and (relatively) consistent plasma environment in which experiments are
performed. The physical subsystems that make up 3DPX will be discussed in Section 3.1.1;
these systems include the vacuum vessel (Section 3.1.1.1), the gas flow (Section 3.1.1.2), the
plasma source (3.1.1.3), and the dust (Section 3.1.1.4). Section 3.1.2 contains a qualitative
description of the plasma conditions and environment in which the experiments described in
Chapter 4 were performed.
3.1.1: 3DPX subsystems
3.1.1.1: Vacuum vessel
The experimental volume of 3DPX is found within two stainless steel vacuum crosses
connected at International Organization for Standardization (ISO) 100 flange surfaces (which
53
Figure 3.1: (a) Photograph of the 3DPX device. (b) Schematic view of the connected vacuum crosses. (c) Detailed
schematic of the "experiment end" of 3DPX.
feature a nominal tube diameter of 102 mm) compression sealed with orings, a photograph can
be seen in Figure 3.1 (a) and the device is shown schematically in Figures 3.1 (b) and (c). The
section labeled "Experiment End" in Figure 3.1 (b) is shown in more detail in Figure 3.1 (c), and
is easily identifiable in the photograph due to the large rectangular window. The experimental
end of 3DPX is a custom sixway vacuum cross made of stainless steel, the total length of this
section is approximately 17", the inner surface of the cross has a nominally circular cross section
with a diameter of 4". The large window seen in Figures 3.1 was included to allow optical
diagnostic access, the window is approximately 10" long and 3" high, as indicated in Figure 3.1
54
(c). The section labeled "Service End" in Figure 3.1 (b) is a standard 5way vacuum cross with
the same inner cross section as the experimental end.
The large window on the service end of the vessel forms a vacuum tight seal to the main
chamber via a compressed oring. The remaining ports of the two crosses are terminated with
either acrylic windows (for optical access), stainless steel flanges with feedthroughs (to allow
electrical connections to electrodes or gas system access), or with blank stainless steel caps. The
remaining subsystems and diagnostics must all gain access to the interior volume of the
chamber through one of the flange terminations. The fact that the flange terminations are all of
the standard "off the shelf" variety provides for a good deal of flexibility to the experimentalist;
the vacuum vessel basically serves as a physical superstructure on which an individual
experiment can be designed and carried out.
3.1.1.2: Gas flow
The gas flow system consists of two basic systems: Pumping (outflow, or "sink") and
neutral gas inflow ("source"). The vacuum pump used for 3DPX is an ULVAC GLD040 series
oilsealed mechanical roughing pump. The pump is connected to the vacuum vessel structure via
a 1.5" diameter reinforced plastic tube which is connected to the "service" end of the experiment
at one of the flanges. With all vacuum terminations sealed the pump provides a base neutral gas
pressure of approximately 510 mTorr. A neutral gas is introduced into the vacuum vessel
through an additional connection found on the service end of the experiment. A tank of ultra
high purity argon is connected in series with the vacuum chamber and an MKS Type 1179A
mass flow controller. The mass flow controller allows a user selectable flow rate of neutral
argon gas to be introduced into the vacuum vessel. The flow rate of the mass flow controller is
55
adjusted until the source flow rate and sink rate (from the pump) are balanced, which results in a
constant pressure of neutral gas within the experimental volume. Typical neutral gas pressures
for experiments within 3DPX range from 50 mTorr to 400 mTorr (approximately 6.6 to 53.3 Pa).
The mass flow controller gas flow rate is adjusted with a userselectable bias voltage specified
within a custom LabView program developed for 3DPX and is monitored with a thermocouple
based pressure gauge manufactured by the Kurt J. Lesker corporation (model KJC6000).
3.1.1.3: Plasma source
The plasma source used in 3DPX is a dc glow discharge system. This type of plasma
source is, conceptually, quite simple; an electrode is electrically biased with respect to a different
electrode. When the resulting electric field is large enough to overcome the neutral gas
ionization energy the outermost electrons are ripped off of an initially neutral atom resulting in a
free electron and a positively charged ion. The electrode configuration used for this experiment
can be seen in Figures 3.2 (a) and (b), the circular brass disc is two inches in diameter (see
Figures 3.2 (c) and (d) for a closer view). The anode is electrically biased with a Glassman
MK1.5P50L power supply to several hundred volts with respect to the vacuum chamber walls
(see Figure 3.2 (e) for a circuit diagram).
The brass disc acts as the anode in the discharge circuit and draws a constant discharge
current (I
d
) of electrons from the plasma, values of which can be varied between 0 and 50 mA.
The top and lateral surfaces of the anode are covered with a dielectric paste (Torr Seal) to ensure
that the collected current of electrons nominally originates in the region below the anode (this
paste is the white substance seen in Figure 3.2 (d)). The anode is suspended in the experimental
volume with a 1/8" diameter hollow rod which is fed through the flange on the top of the
experimental section of the vacuum vessel and the anode is electrically connected to the power
56
Figure 3.2: (a) and (b) Photographs of the anode inside of the vacuum vessel in the position used for the experiment.
(c) and (d) Photographs of the anode outside of the chamber. (e) Circuit diagram for the dc glow discharge plasma
source.
supply with an insulated wire inside of the support rod. The interior of the chamber also features
a flat stainless steel surface (as seen covered with a layer of particulate matter in Figures 3.2 (a)
and (b)) which is in direct electrical contact with the vacuum chamber walls, the chamber walls
57
and the metallic surface are connected to the laboratory and power supply ground and
collectively act as cathode in the discharge circuit. The flat surface is a feature that has been
added in an attempt to regularize the structure of the electric field within the experimental
volume. The spacing between the bottom surface of the anode and the flat surface on which the
dust sits was 4.7 cm in the experiment described below. The exact separation distance,
orientation of the anode asymmetry, and positions of the metal sheets that form the flat surface
all have large effects on the characteristics of the plasma that is produced. The anode/cathode
orientation used for these experiments was selected for the simple reason that it produced fairly
large dust clouds that were confined to the region near the anode over a fairly wide range of
discharge currents and neutral gas pressures.
3.1.1.4: The dust
The final "subsystem" of the experiment is the particulate matter that acts as the dust
component of the plasma. For these experiments the dust is comprised of 1.5 ? 0.5 ?m radius
silica particles. The grains are nominally spherical and have individual masses of 3.60 ? 10
14
?
0.13 ? 10
14
kg. The particles were acquired from Potters Industries Inc., where they are
manufactured to be used as a dopant in molded plastics fabrication to increase the durability and
strength of the manufactured parts. The dust is introduced into the system by placing an
approximately 5 mm thick flat sheet of the particles on the metallic surface below the anode, as
can be seen in Figures 3.2 (a) and (b).
When a plasma is initiated in 3DPX a small fraction of the dust grains become levitated
when the discharge arcs (transient lightninglike electric field enhancements most likely due to
ionization instabilities). After the initial "kick" up into the plasma environment a grain acquires
charge by collecting a fraction of the ambient ions and electrons, because the electron population
58
has a higher mobility than the ion population the grain collects more electrons than ions which
results in a net negative charge. A simplified qualitative description of the process in which a
grain becomes levitated in the bulk plasma is as follows. If the trajectory of a grain is such that it
falls back into the region near the anode, the force of gravity pulling the grain down towards the
dust pile may, in some cases, be balanced by the upward electrostatic force in the sheath/pre
sheath region above the dust pile. In the vast majority of such events a grain moves with a
velocity that is too large for the electrostatic force to slow the particle to zero velocity before it
passes through the sheath and reenters the dust pile. The dust grains that form the dust cloud
and become longterm participants in the plasma are rare cases, the success rate of injecting
particles into the cloud can be increased by initiating cloud formation at higher neutral pressures
where the dustneutral drag force assists the process by slowing the grains via collisions as they
fall.
3.1.2: Description of the plasma environment
DC glow discharge plasmas are widely used in both industrial and research settings
because they are conceptually simple and fairly easy to produce. One of the main benefits of the
dc glow discharge plasma source is that "islands" of steady state operation can be found within
the parameter space defined by the electrode configuration, neutral pressure, and discharge
current. With the electrode configuration shown in Figures 3.2 (a) and (b) the plasma was found
to be more or less steady state for pressures ranging from ~70  250+ mTorr with discharge
currents ranging from 2  15 mA. The criteria for a "steady state" system consisted of two times
scales: The long time scale criteria were that a dust cloud could form, remain at a given location,
and retain its approximate size for time scales of several hours and the shorter time scale
criterion was that the plasma discharge circuit did not arc, which causes the cloud to either
59
Quantity Value Quantity Value
n
i
10
14
m
3
n
dust
10
10
m
3
q
dust
4400 e

m
dust
3.6 x 10
14
kg
r
dust
10
6
m ?
?,?
10
3
m
?
?,?
10
4
m ?
?,?
10
3 m
k
B
T
i
= k
B
T
n
1/40 eV k
B
T
e
5 eV
?
??,?
3.8? 10
12
eV/m
3
?
??,?
7.5? 10
14
eV/m
3
?
??,?
6.2? 10
15
eV/m
3
Table 3.1: Approximate plasma and dust parameters.
disappear completely or to momentarily "jump" out of position and subsequently "fall" and
violently oscillate about the previous equilibrium position. For the experiment described in
Chapter 4 the electrode configuration was the same as described in Section 3.1.1.3, the argon
neutral pressure was 132 mTorr, and the discharge current was 5.37 mA.
The characteristics of the background plasma species (ions, electrons, and neutral gas) are
difficult to diagnose when dust is present in the system; due to complications related to surface
contamination for probe based diagnostics and because of the low plasma density limitations for
optical diagnostics. Probe based measurements in similar discharge configurations, but without
dust present in the chamber, have given the "typical" approximate plasma parameters listed in
Table 3.1. Acquisition of more accurate background species parameters, particularly for dc glow
discharge plasmas and more generally for any type of discharge in which a significant dust
component is observed, is an ongoing avenue of research in the dusty plasma physics community
and is beyond the scope of this work.
The dust component of the plasma, on the other hand, can readily be investigated within
3DPX. The specific diagnostic technique used for this work, particle image velocimetry, will be
60
discussed in detail in Section 3.2 but several general observations will be noted here in order to
facilitate the discussion that follows. First, the "steady state" criteria described above ensure that
the dust component remains more or less stable for time scales on the order of hours (and likely
longer), this means that any observation or diagnostic can be performed over a period of time
(several hours) on a single dust cloud. This removes the ambiguity associated with the
assumption that repeated formation/destruction/reformation cycles result in identical ensembles
of dust particles. The experimentally selectable parameters used in the experiments described
below (electrode configuration, neutral pressure, and discharge current) were also carefully
chosen to ensure that the dust component formed a single cloud in a single location within the
experimental volume, this was done so that the entire structure could be measured quickly and
efficiently. Finally, the discharge parameters were found where the cloud was as quiescent as
possible, a "stable" cloud, free of waves and other obvious large scale transport, was chosen in an
attempt to simplify the subsequent analysis.
3.2: Particle Image Velocimetry
Dusty plasmas containing large (?
?
? 1 2? ??) dust grains are a novel system in
experimental plasma physics for two main reasons: The large particle size allows direct optical
imaging of the particles and, secondly, the relatively small charge to mass ratio means that the
dust grain dynamics occur on much longer time scales than those of the ion or electron plasma
components. In this section the application of an optical diagnostic known as Particle Image
Velocimetry (PIV) is discussed
39?41
. Due to the fact that PIV is somewhat uncommon outside of
the fluid physics community the measurement technique will be described beginning with its
most basic form, from which the method that is used in these experiments can be explained more
clearly. This section will be organized as follows: Section 3.2.1 gives a brief introduction to
61
PIV, Section 3.2.2 describes two dimensional PIV (2DPIV), Section 3.2.3 contains a description
of stereoscopic PIV (stereoPIV), and Section 3.2.4 discusses the specific stereoPIV system
used in the experiment described in Chapter 4.
3.2.1: Introduction
PIV is a diagnostic technique developed within the fluid physics community to measure
the flow of fluids in a manner that is minimally invasive. Before the advent of PIV, a physical
probe had to be inserted into a fluid flow in order to obtain a quantitative measure of the fluid
velocity at the location of the probe. This practice was not ideal, as the presence of a relatively
large physical probe changes the flow pattern of the system in the region of the measurement.
The solution to this issue was to seed the flow with neutrally buoyant particles ("tracers") which
could be imaged optically. The tracer particles follow the flow pattern of the system and can be
imaged when illuminated with laser light. If appropriate tracer particles are chosen the flow is
not dramatically altered. In order to obtain a measure of the flow pattern in the system of interest
the tracer particles are imaged twice with a known time delay (?t
laser
) between images. By
comparing the two recorded images the displacement vector for groups tracer particles (and thus
the fluid displacement vector) can be calculated, by dividing the measured displacement by the
known time delay between laser pulses the velocity field of the fluid is obtained.
The application of the PIV diagnostic to dusty plasma systems was pioneered at Auburn
University
8,42?44
. The application of PIV has grown beyond Auburn and is now used by many of
the research groups actively studying dusty plasmas. The key advantage of the PIV diagnostic,
when applied to dusty plasmas, is that one of the major difficulties in the standard application of
the technique, the selection of an appropriate tracer particle, is eliminated. In a dusty plasma
62
system the dust itself is the tracer particle. While this makes the experimental process much
more straightforward one must be aware that many of the techniques used by the fluid physics
community to analyze PIV data are not applicable to the measurements of dusty plasmas. This is
due to the fact that in the application within the fluids community the tracer particles are only
present in order to give information about the background fluid, meaning interesting quantities
such as the temperature of the background fluid are calculated with an a priori knowledge of the
properties of the background fluid in question (i.e. the fluid density, diffusivity, etc.). In the
application of PIV to duty plasmas there is no background fluid, the tracer particles themselves
are the system of interest. The application of PIV to dusty plasma systems has the advantage of
adopting a mature and established measurement technique, but has the drawback that much of
the analysis that takes place after velocity vector computation must be carefully and
systematically reconsidered.
3.2.2: Two dimensional PIV
In this section the process of obtaining a two dimensional velocity field using 2DPIV is
discussed, while this specific diagnostic was not used in this work the hope is that by starting
with the most basic form of PIV the more complicated cases discussed in Sections 3.2.3 and
3.2.4 will be more clearly and easily understood. It is noted that there are numerous techniques,
methods, and variations of the 2DPIV measurement scheme, in what follows the socalled
"double frame, double exposure" method is discussed.
As mentioned in Section 3.2.1, the 2DPIV measurement is remarkably simple in
principal: Calculate the displacement field of a seeded flow and divide by ?t
laser
. In practice the
process of finding the displacement field can become somewhat involved, this section describes
63
Figure 3.3: Schematic drawings of the 2DPIV configuration for a single particle. (a) As viewed from above. (b) As
viewed from the plane of the camera, the blue box indicates the camera's field of view.
the process algorithmically. In this section, and in Section 3.2.3, it will be assumed that there is
some general fluid flow that is seeded with appropriately chosen tracer particles.
A two dimensional velocity field on a given spatial plane is found by orienting a laser
sheet so that when the laser fires all of the tracer particles within the plane are illuminated. The
laser is then set to fire two very short pulses separated in time by a delay ?t
laser
. The firing of
each laser pulse is synchronized to a CCD camera oriented so that the CCD array is parallel to
the laser sheet (as in Figure 3.3 (a), for a single particle in the laser sheet, and on the left hand
side of Figure 3.4, for multiple tracer particles). The camera records the intensity of the laser
light scattered by the tracer particles as two separate images (one image for each of the laser
pulses), represented in Figure 3.3 by the positions of the red dot ??
0
and ??
1
. The first image
(frame 0, taken at time t
0
) and second image (frame 1, taken at time t
0
+ ?t
laser
) are each
"snapshots" of the spatial distribution of the tracer particles in the laser plane. The time delay
?t
laser
must be chosen such that all tracer particle motion that occurs in the time between the two
pulses can be approximated as linear, typical values are in the 1  1000 ?sec range, depending on
the flow in question.
Each image is then divided into small regions known as interrogation cells ("cells"), as
seen on the left side of Figure 3.4. The cell size is chosen so that each is large enough to contain
64
Figure 3.4: Flow chart for a single 2DPIV calculation for an interrogation cell containing multiple tracer particles
(the green boxes on the left). (Courtesy of LaVision)
multiple tracer particles but small enough that meaningful spatial resolution is maintained.
Factors such as camera lens, separation distance between the laser plane and camera, and tracer
particle size must be considered when choosing the size of the interrogation cells. Figure 3.4
shows an example of two frames divided into interrogation cells with a seeded flow. The next
step in the process is to calculate the most probable average displacement vector for the tracer
particles within each of the interrogation cells. These quantities are computed by finding the
peak of the spatial correlation function between the frame 0 and frame 1 images within each
discrete cell. The discrete correlation function is given by:
? (dx, dy) = ? ?
0
(?,?)?
1
(? + dx,? + dy) for
??
2
?
(?,?)=0
< (dx, dy) <
?
2
3.1
where n is the spatial step size (the number of pixels on each side of the cell box) and where I
0
and I
1
are the light intensity fields recorded by the CCD array in frames 0 and 1, respectively. In
words, the intensity function in frame 1 is shifted in both the x and y directions (by dx and dy)
and compared to the intensity function in frame 0; the correlation, C(dx,dy), is a measure of the
similarity of the two intensity functions due to a shift of I
1
by (dx,dy). This process is
65
mathematically identical to finding the complex two dimensional Fourier transformation for both
I
0
and I
1
(?
0
~
and ?
1
~
), performing a complex conjugate multiplication (?
0
~
? ?
1
~
*
), and then calculating
the inverse Fourier transformation of the resulting product. The correlation function is a peaked,
as shown in the center of Figure 3.4. The most probable average displacement of the tracer
particles within the interrogation cell is given by the peak location of the correlation function.
This process is then repeated for each interrogation cell on the image plane, giving an
instantaneous most probable average displacement field for the flow at time t
0
. The
displacement field is then converted to a velocity field by dividing each interrogation cell's
instantaneous most probable average displacement by ?t
laser
.
The major limitation of the two dimensional version of the PIV diagnostic is that it can
only measure the velocity field of the flow as projected onto the plane of the camera. This
projection becomes an issue when the flow in question is three dimensional; for example, if the
"real" velocity vector for a fluid element is
??
actual
= ?
?
??
?
+ ?
?
??
?
+ ?
?
??
?
, 2DPIV only gives the
two dimensional vector ?? = ?
?
??
?
+ ?
?
??
?
. To correct for this issue we use stereoscopic PIV.
3.2.3: Stereoscopic PIV
Stereoscopic PIV is a natural extension of the two dimensional version of PIV, this
measurement technique uses two cameras focused on the same spatial region of a flow. The
cameras are oriented at known angles with respect to the laser sheet normal; a schematic of such
a configuration can be seen in Figure 3.5 (a). StereoPIV has three coordinate systems that must
be considered: The coordinate system defined by the laser sheet (the lab frame), the coordinate
system defined by the CCD array of camera 1 (frame 1), and the frame of camera 2 (frame 2). It
is noted that in 2DPIV, as described in Section 3.2.2, the camera is oriented perpendicular to the
66
Figure 3.5: (a) Top view of the stereoPIV camera and laser sheet orientation. (b) The displacement from (a) as
seen by camera 1. (c) The displacement in part (a) as seen by camera 2.
laser sheet, meaning the lab frame and the camera frame are identical. As in the 2D case the
cameras are synchronized to record images of the seeded flow when the flow is illuminated by
the two laser pulses separated by a delay ?t
laser
; but, because the cameras are oriented at different
angles with respect to the laser sheet they will record different two dimensional displacement
fields, as illustrated in Figures 3.5 (b) and (c). It is through the geometric rotations of these
different two dimensional displacement vectors that the three dimensional displacement fields
are obtained.
In order to measure the velocity vectors of tracer particles in units of meters per second
the spatial orientation angles and absolute scaling (camera pixels/meter) of the cameras must be
carefully measured. This is accomplished with the use of a calibration plate, as seen in Figure
3.6 (a). The plate is an array of regularly spaced circular white dots on two planes separated by 1
mm. When the cameras are focused on the plate, the angular orientation and spatial scaling of
each camera can be found by using the known size and separation of the dots on each of the
planes.
67
Figure 3.6: (a) Image of the calibration plate used to determine the absolute orientation and spatial scaling of the
PIV cameras. (b) An image (courtesy of LaVision) illustrating the view of a similar calibration plate without (with)
a Scheimpflug adapter is shown on the left (right). (c) Images showing the view of the calibration plate after de
warping for the calibration used in the experiment. (d) Images showing the superposition of the dewarped images
of the calibration used in the experiment after rotation into the laboratory coordinate system. The image on the left
shows the overlap of the dots in the far (recessed) plane of the plate and the image on the right shows the overlap for
the nearer calibration plate plane. The appearance of multiple dots in the superimposed images for the secondary
plane in each image shows the parallax that is the result of the different camera positions.
Two issues arise due to the fact that the cameras are oriented at angles with respect to the
laser sheet which must be corrected before the displacement fields can be calculated. The first
issue is that the laser plane and the plane of the CCD array (and its focusing lens) are not
parallel, because of this the image is focused only in the region where the focus plane of the
68
camera and the laser plane intersect, the effect can be seen on the left hand side of Figure 3.6 (b).
To correct the focus problem (which is known as the Scheimpflug principal) a Scheimpflug
adapter is placed between the camera lens and the laser plane. When properly adjusted the
Scheimpflug adapter allows the camera to focus on the laser plane across the cameras entire field
of view, as seen on the right hand side of Figure 3.6 (b). The second issue that arises is that the
magnification (zoom) level is not constant across the camera's field of view, as can be seen by
the different dot sizes on the calibration plates in Figure 3.6 (b). During the calibration process a
dewarping function is calculated that corrects for the distortion. Dewarped images of the
calibration plate used in the experiment described below are shown in Figure 3.6 (c). With the
Scheimpflug adapter and dewarping function the images from both cameras can be rotated into
the laboratory coordinate system. Figure 3.6 (d) shows a superposition of the same calibration
plate images as in Figure 3.6 (c) after such a rotation; on the left the dots from the recessed plane
of the calibration plate overlap at the intersection of the red grid lines, the dots on the nearer
plane of the calibration plate do not overlap, showing the effect of the parallax that results from
the two different viewing angles. The image on the right hand side of Figure 3.6 (d) shows the
same but for the other plane of the calibration plate. The calibration is finalized by a process
known as "selfcalibration." This process is based on the idea that if an object in the laser plane
is imaged by both camera 1 and camera 2 at the same time the particle location should be
identical when dewarped and rotated into the lab frame. If the rotated images from camera 1
and camera 2 differ the dewarping functions are corrected and the images are rotated to the lab
frame with the corrected dewarping function. The process is repeated until no further
improvement can be made.
69
With a valid camera calibration, the actual process of calculating the three dimensional
displacement field is very similar to the 2DPIV procedure described in Section 3.2.2. First, each
camera image is dewarped and a two dimensional displacement vector is calculated within an
interrogation cell in the coordinate systems of each camera. The two dimensional displacement
vectors from each camera, ?
?
c1
and ?
?
c2
, are then related to the three dimensional lab frame
displacement vector, ?
?
, in each interrogation cell through a rotation defined by:
?
?
ci
= ?
?
?,?
?
?
3.2
where i refers to either camera 1 or 2 and ?
?
?,?
is the rotation matrix that relates the lab and
camera coordinate systems. This process gives an overdetermined system of four linear
equations with three unknowns (d
x
, d
y
, and d
z
):
?
c1,?
= ?
?
cos (?
1
) + ?
?
sin (?
1
)
?
c2,?
= ?
?
cos (?
2
) + ?
?
sin (?
2
)
?
c1,?
= ??
?
sin (?
1
)cos (?
1
) + ?
?
cos (?
1
)cos (?
1
) + ?
?
sin (?
1
)
?
c2,?
= ??
?
sin (?
2
)cos (?
2
) + ?
?
cos (?
2
)cos (?
2
) + ?
?
sin (?
2
)
3.3
where the angles ?
?
and ?
?
are the angles relating the camera and laser sheet coordinate systems.
The lab frame displacement vector is found by fitting the system of equations using the least
squares method, with the constraint that the uncertainty is distributed evenly across the three
vector components. After this is carried out for all of the interrogation cells the displacement
field is converted to a velocity field by dividing by ?t
laser
.
70
3.2.4: Stereoscopic PIV in 3DPX
The stereoPIV system used on the 3DPX experiment consists of two LaVision Imager
Intense cameras synchronized to a New Wave Research brand Solo laser (Nd:YAG) with a
LaVision PTU 9 programmable timing unit (PTU). The laser illuminates an approximately 2
mm thick plane of the constant ?? plane inside the 3DPX vacuum vessel. The laser and the
cameras are mounted on a translation stage which allows the point of camera focus and laser
plane to be moved to all locations visible through the large window shown in Figure 3.1 (a). The
translation stage, cameras, laser, and PTU are all controlled by computer using the DaVis
software (version 7.2) made by LaVision. The cameras have fields of view that are 1419 pixels
in the ?? (horizontal) direction and 994 pixels in the ?? (vertical) direction. With the configuration
described above, the spatial scale, after dewarping, is approximately 37.60 pixels/mm giving a
field of view that is 3.77 cm in the horizontal direction by 2.64 cm in the vertical direction. A
typical image of a constant z slice of a dust cloud confined below the anode is shown in Figure
3.7 (a). Typical interrogation cell sizes range from 12 x 12 pixels to 128 x 128 pixels, with the
cell size chosen such that >10 dust grains are in each cell within the cloud region. For the
analysis described below the interrogation cells were chosen to be 64 pixels square (shown as the
yellow grid in Figure 3.7 (b)). The four camera images recorded within the cell outlined in red in
Figure 3.7 (b) are shown in Figure 3.7 (c). These four images, along with the fact that ?t
laser
=500
?sec, are then used to calculate the three dimensional velocity vector for the cell, which is shown
projected onto the {??,??} plane in the black box outlined in Figure 3.7 (d). The complete velocity
vector field for the dust cloud cross section pictured in Figure 3.7 (a) is shown in Figure 3.7 (d).
71
Figure 3.7: (a) Full camera frame view of a dust cloud cross section. (b) Detailed view of the camera region
containing the dust cloud (the region outlined in yellow in (a)), the yellow grid indicates the location of the
individual interrogation cells. (c) The four camera images of the cell outlined in red in (b). (d) The instantaneous
velocity field measured for the cloud cross section, projected onto the ???
?
,??
?
? plane.
The process described above gives a "snapshot" of the velocity field for the dust
component of the plasma at the time of the measurement, to investigate the distribution of
velocities the process is repeated many times (1500 for this experiment). Figure 3.8 shows the
first ten of the 1500 three dimensional velocity vectors measured within the volume element
highlighted in Figure 3.7, the velocity vector corresponding to the cell outlined in red in Figure
3.7 (b) is shown as blue in Figure 3.8. The gray arrows in Figure 3.8 show the projections of the
72
Figure 3.8: Time series showing the first ten three dimensional velocity vector measurements for the cell
highlighted in Figure 3.7. The blue vector is that which was calculated from the four camera images shown in
Figure 3.7 (c). The gray arrows on the back (bottom) plane are projections of the three dimensional vectors onto the
{??
?
,??
?
} ({??
?
,??
?
}) plane.
measured three dimensional vectors onto the {??
?
,??
?
} (the vertical, back, plane) and the {??
?
,??
?
}
(the horizontal, bottom, plane) these are intended to show the importance of measuring all three
velocity vector components. After the 1500 measurements of a given cross section are
completed the translation stage holding the cameras and laser are moved 2 mm in the ?? direction
and the process of recording the 1500 velocity field snapshots is repeated for all planes that
contain dust.
The final step involved with the acquisition of the velocity vector data is the
quantification of the uncertainty in the stereoPIV measurement, this is done with the socalled
zero displacement test. The zero displacement test measures a three dimensional velocity vector
in the same way as described above, but with the time between laser pulses set to the minimum
(0.6 ?sec). With the very short time between laser pulses the dust grains do not have a chance to
move between the recordings of the two camera images. The displacement field is calculated as
described above, giving the minimum measureable displacement vector for the diagnostic. The
73
displacement vectors are then converted to m/sec by dividing the displacement vector values by
the laser separation time used in the acquisition of the actual data set (that is, the zero
displacement test displacement vectors are divided by 500 ?sec). When repeated many times
(again, 1500) the zero displacement test gives the minimum measurable velocity vector
distribution in each volume element across the entire dust cloud structure.
74
Chapter 4: Phase space distribution measurements
4.1: Introduction to the experiment
This chapter contains a description of the experiment and analysis which were performed
in order to measure the spatially resolved phase space distribution of a dust cloud within the
3DPX device. The stereoPIV diagnostic (referred to as "PIV" in what follows) was used to
measure both the average distribution of the dust cloud in configuration space and to construct
velocity distributions for small volume elements across the observed structure. To make such
measurements the 3DPX experimental volume was filled with neutral argon gas to a pressure of
132 mTorr (~17.6 Pa). A bias voltage was then applied to the anode, through which a constant
current of 5.37 mA was drawn throughout the experiment. A dust cloud was allowed to form,
grow, and come to an apparent equilibrium during a time period of approximately two hours;
after the growth and stabilization period the cloud was imaged with the stereoPIV diagnostic as
described in Section 3.2. 1500 PIV measurements were taken of each cloud cross section, with a
laser pulse separation time of 500 ?sec, followed by 1500 measurements with a pulse separation
time of 0.6 ?sec (for the zero displacement test). After the two sets of PIV measurements were
recorded for a given spatial cross section the translation stage carrying the PIV laser and cameras
was moved 2 mm in the ?? direction and the process was repeated for all of the planes of constant
?? that were observed to contain dust. In total there were 13 such cross sections resulting in a
total of 39,000 PIV "snapshot" measurements across the dust cloud structure. Each PIV
75
measurement contained 330 velocity vector measurements (for a total of ~12.9 million vectors)
and a total of approximately 156,000 camera images. The PIV velocity vector measurements
were processed with the commercially available PIV analysis software package DaVis (version
7.2) produced by the LaVision corporation.
The remainder of this chapter proceeds as follows: Section 4.2 discusses the PSD
construction process for a single volume element within the dust cloud, Section 4.3 describes the
process that was used to compare the standard Maxwellian distribution and trinormal models of
the velocity space and determine that the trinormal probability distribution function is required
to accurately model the velocity space, and Section 4.4 gives an overview of the resulting
spatially resolved PSD parameter value fields.
4.2: Single volume element
In order to facilitate the discussion of the data analysis process used to construct the
spatially resolved PSD for all 853 volume elements within the dust cloud the analysis procedure
will be demonstrated, in detail, for a single volume element. The volume element that will be
used as an example in this section is the same that was highlighted in Section 3.2.3 and shown in
Figure 3.7.
The PSD contains configuration and velocity space components which are combined to
give the overall PSD, as in Equation 2.1. The process used to determine the configuration space
component, the number density n
d
, is discussed in Section 4.2.1. Section 4.2.2 contains a
description of the method used to find the distribution function parameters for the two models of
velocity space. Section 4.2.3 discusses the test that was used to determine that the trinormal
76
Figure 4.1: (a) On the left, an example of a single camera image of the dust grains observed within the volume
element. On the right, the particles which were identified by the particle counting software are shown as different
colors. (b) A histogram of the number of particles identified within the volume element for all of the 1500 camera
images obtained during the PIV acquisition process.
velocity space distribution function model is required to describe the observations, and Section
4.2.4 contains brief conclusions.
4.2.1: Configuration space
The dust number density was found by examination of the camera images obtained with
the PIV diagnostic. Each PIV measurement results in four camera images, as was seen in Figure
3.7 (c). The image recorded by camera two, with the first laser pulse, was extracted for each of
the 1500 measurements and partitioned into subimages on the same 64 pixel by 64 pixel grid
which was used for the velocity vector reconstruction (as in Figure 3.7 (d)). Figure 4.1 (a) shows
an example of one such subimage for the volume element, the unprocessed image is shown on
the left and the individually identified particles from the same image are shown on the right hand
side of the figure. The particle identification was performed within the computer program
Mathematica. The number of identified particles was counted within the subimage for each of
the 1500 measurements of the volume element. Figure 4.1 (b) shows a histogram of the number
of particles observed within the volume element; the results showed that an average of 48.4?4.8
77
particles were present. The dust number density, n
d
, has units of particles per cubic meter, to
convert the number of particles, N
d
, to the number density, n
d
, we simply divide the average
number of particles by the configuration space volume of the element, 5.8?10
9
m
3
, to arrive at
the resulting number density: n
d
= 8.3?10
9
? 0.8?10
9
m
3
.
4.2.2: Velocity space
The process used for construction of the velocity space component of the PSD is more
complicated than for the number density field. The result of the PIV measurements is two data
sets within each of the volume elements; the first data set is a table of the 1500 three dimensional
velocity vectors that were measured with a laser pulse separation time of 500 ?sec (collectively
referred to as "D"). The second data set is a table of the 1500 velocity vectors that were used to
quantify the error, obtained with a PIV laser pulse separation time of 0.6 ?sec (the data set "Z").
The process used to calculate each of the individual velocity vectors was discussed in detail
within Section 3.2. In contrast to the configuration space distribution, the velocity space
distribution is modeled with a continuous distribution function. The two distribution functions
examined here are the canonical three dimensional Maxwellian distribution function and the tri
normal distribution function, both of which were described in Section 2.2. The remainder of this
section will describe the process used to calculate the parameters for each of the distribution
functions.
4.2.2.1: Calculation of the drift velocity vector
The drift velocity of the dust grains within the volume element is the simplest of the
velocity space distribution parameters that must be calculated. The calculation merely involves
finding the average value for each of the three velocity space components in the data set D, the
78
result is the drift vector, u??. The drift velocity appears in an identical fashion within both the
Maxwellian and trinormal models of the velocity space, because of this the quantity is
calculated in the same way for both distributions. The fact that the calculation process is the
same for both models means that it is unnecessary to incorporate the drift parameters into the
model comparison process discussed later in this Chapter. The values of the drift velocity
components within the example volume element were found to be nearly zero: ??? =
(0.00002 ? 0.00010)??
?
+ (0.00004 ? 0.00013)??
?
+ (0.00001 ? 0.00022)??
?
m/sec.
4.2.2.2: Single particle velocity space
The width characterization of the velocity space distributions obtained from the PIV
measurements must be considered with more care than the offset (drift velocity) characterization.
The individual velocity vectors that the PIV diagnostic system calculates are the result of
correlation integrals; this type of data processing has an inherent bias that favors low
displacement magnitude values (and thus low velocity magnitude values). In effect, this bias
acts to reduce the variance of the actual singleparticle displacement distribution, ?
?
2
, and instead
gives variance values, ?
???
2
, that are smaller in magnitude by a variance ratio factor that is
dependent on the number of particles that are present in the correlation calculation process. An
analytic expression can be found for the correction factor if the exact form of the correlation
function is known for all of the measurements. Estimates of the correlation function from
measured data are, themselves, inherently biased, meaning that obtaining an analytic result for
the correction with the PIV process described in Section 3.2 is impossible.
The issue of correcting for the variance reduction has been addressed by simulation and is
described in detail elsewhere
8
. The simulation study found that the variance values obtained
79
Figure 4.2: A plot of the correction factor, CF(N
d
), that is used to convert the velocities measured by the PIV
diagnostic to the velocity of a single particle within the volume element.
from distributions of the PIV measured velocity vectors were less than the known singleparticle
velocity variance values that were input into the simulation, as expected. It was also shown that
correct variance values could be recovered by mapping the PIV (multipleparticle) velocity space
into a singleparticle velocity space with the use of the correction factor, shown in Figure 4.2,
which depends on the number of dust grains observed within the velocity vector reconstruction
region. If effect, the mapping process rescales all three components of the measured velocity
vector, ??
?,PIV
, into the singleparticle velocity space:
??
?
= (??
?,PIV
????) CF (?
?
(??))? 4.1
where ??
?
is the vector that has been mapped into the singleparticle velocity space and
CF ??
?
(??)? is the correction factor value. For simplicity, from this point forward all references
to "measured velocity vectors" refer to ??
?
, the velocity vector which has already been mapped
into the singleparticle velocity space. It is also noted that the drift velocity calculation is
performed before the data are rescaled. It has been shown, through the analysis of both
simulated dusty plasma systems and with the application of PIV to experiments where other
80
velocimetric methodologies are available
8,43
, that the mean values of the measured velocity space
distributions, as measured by PIV, give the true drift velocity.
4.2.2.3: Visualization of the velocity space distribution
With the corrected values for the velocity vectors we can proceed toward the calculation
of the distribution function parameters. The measured velocity vectors within this volume
element are shown as histograms in Figure 4.3, they are shown to help visualize the distributions
and parameters that will be discussed in what follows. Figure 4.3 (a) shows D as three one
dimensional histograms in the "experiment" (or laboratory) coordinate system (?? = ??
?
,?
?
,?
?
?).
A cursory examination of the histograms in 4.3 (a) would seem to imply that the velocity space is
more or less isotropic, because of the relatively close agreement between the fitted standard
deviation values of the two models (indicated on each of the histogram plots). Figures 4.3 (b),
shows the same data rotated into the principal axis coordinate system (?? = {?
1
,?
2
,?
3
}). The
difference in the two sets of histograms qualitatively illustrates that the apparent isotropy seen in
(a) is not a true feature of the distribution. Figures 4.3 (c) shows two dimensional histograms of
D in the ???
?
,??
?
? subspace, where the elliptical nature of the velocity space distribution can
clearly be seen.
Figure 4.4 shows similar histograms for data set Z. Figure 4.4 (a) shows the error
distribution as one dimensional histograms, in the laboratory frame, which have the same bin
size (resolution) as in Figures 4.3 (a) and (b), from which it can seen that the error distribution
has a much smaller standard deviation than the distributions in D. Figures 4.4 (b), (c), and (d)
81
Figure 4.3: Various histograms of the measured velocity vectors. (a) The three one dimensional velocity space
histograms in the laboratory coordinate system. Note that the three distributions appear to have approximately the
same width. (b) The three one dimensional velocity space histograms after rotation into the principal axis
coordinate system. It can be seen here that the isotropy seen in (a) is not a true feature of the velocity space. (c)
Examples of a histogram in two velocity space dimensions ???
?
, ??
?
?, the two plots show the same data and have the
same color scale. The velocity space anisotropy can be seen clearly in this subspace.
show the error distribution for the three two dimensional velocity subspaces, these histograms
have smaller bins (higher resolution) which allows visualization of the shape of the error
distribution. It is noted that the distributions in Figures 4.3 (c) are not simply different than the
distributions in Figure 4.4 (b) by a scaling factor, but have intrinsically different shapes.
82
Figure 4.4: Various histograms of the error vector measurements. (a) The three velocity space histograms for the
error distribution measurements in the laboratory coordinate system. The histograms have the same resolution (bin
size) as those in Figure 4.3, to show that the error distribution measurements have a much narrower width than those
measured for D. (b), (c), and (d): The three two dimensional histograms of the error distribution measurements.
The histograms here have higher resolution (smaller bin size) than in (a) so as to show the structure. It is noted that
the histogram in (b), for the nullnullnull
null
,nullnull
null
null subspace, is different than the histograms in Figure 4.3 (c) in both width (scale)
and shape.
4.2.2.4: Parameter value calculation: Maxwellian distribution model
The process of calculating the parameter values, assuming the Maxwellian probability
distribution function, is very similar to the process outlined in Section 2.3.3.3. Recall that the
three dimensional Maxwellian distribution is given by Equation 2.11:
83
?
?
??
?
,?
?
,?
?
? =
(2??
?
2
)
?3 2?
exp?
?1
2
?
(?
?
??
?
)
2
?
?
2
+
(?
?
??
?
)
2
?
?
2
+
(?
?
??
?
)
2
?
?
2
??
4.2
This is somewhat simplified because of the fact that we have translated our velocity space into
the local rest (zerodrift) frame, as in Section 4.2.2.2, eliminating ??? from the calculation of the
distribution shape parameters. To find ?
?
, the scalar standard deviation, for the Maxwellian
distribution model we assume that ?
?
can be expressed in terms of an optimal value and a
normally distributed uncertainty: ?
?
? ?
?0
? ??
?
. As with the discussion of the one
dimensional example, outlined in Section 2.3.3.3, the likelihood function of D, given ?
?
, the
model ("M"), the error distribution Z, and the associated information is given by:
prob(??
?
,?,?,?) = exp?
?(?
?
??
m0
)
2
2(??
?
)
2
? ?
?(2?)
3 2?
?
??
?
+ ?
m0
2
?
?
??
??
exp?
?1
2
?(??
?
)
?
? ??
?
+ ?
m0
2
?
?
?
?1
? (??
?
)
?
?=1
?
4.3
Where the components of ??
?
are the measured vectors in D and ?
?
is the velocity space
covariance tensor for Z. The error distribution is modelled with the trinormal distribution
function, the individual elements of ?
?
are given by:
?
ij
= ?(?
ki
?< ?
?
>)(?
kj
?< ?
?
>)??
?
??
?
?
?
?=1
4.4
The summation is over the N
z
velocity vectors that make up the data set Z and is the
average value of the ??
?
component of the vectors in Z.
84
Parameter Value
?
?
(8.3 ? 0.8) ? 10
9
m
?3
??? (0.00002 ? 0.00010)??
?
+ (0.00004 ? 0.00013)??
?
+ (0.00001 ? 0.00022)??
?
m/sec
?
?
0.0219 ? 0.0059 m/sec
Table 4.1: PSD parameters with the Maxwellian velocity space model.
To find the optimal parameter value ?
?0
we maximize the likelihood function given in
Equation 4.3. To compute the ?
?0
value that gives the maximum value of the likelihood
function the computer program Mathematica is used, the program includes the function
"FindDistributionParameters" which numerically calculates the optimal parameter value for the
distribution defined in Equation 4.3 by maximizing the likelihood function. The ?
?0
value
calculated in this way for the volume element discussed in this section was found to be ?
?0
=
0.02191 ?/???. With this value in hand, we can find the uncertainty ??
?
by returning to the
likelihood function, prob(D?
?
,M,Z,I), taking two partial derivatives of its natural logarithm
with respect to ?
?0
, and setting the resulting expression equal to zero, as was described in detail
in Section 2.3.3.3. Substitution of the ?
?0
value into the resulting equation gives the uncertainty
in the parameter estimate: ??
?
= 0.00038 ?/???. So, for the volume element in question the
parameters for the Maxwellian distribution function model were obtained by finding the drift
vector and the scalar variance. The full PSD within this volume element with the Maxwellian
model is specified by the parameters listed in Table 4.1.
4.2.2.5: Parameter value calculation: Trinormal distribution model
The calculation of the distribution parameter values for the trinormal probability
distribution function model proceeds in much the same way as for the Maxwellian model. We
begin with the trinormal distribution function, which is given in the local rest frame by:
85
?
TN
??
?
? = (2?)
?3 2?
??
?
?
?1 2?
exp?
?1
2
(??)
?
? ?
?
?1
? (??)?
?
?
= ?
?
?
?
2
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
2
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
2
?.
4.5
The standard deviation and correlation parameters are all assumed to be distributed as an optimal
value with normal uncertainty: ?
?
?
= ?
?
?
0
? ??
?
?
and ?
?
?
?
?
= ?
?
?
?
?
0
? ??
?
?
?
?
. The resulting
likelihood function is given by:
prob ????
?
, TN,?,??
= ?(2?)
3 2?
???
?
+ ?
?
??
??
exp ?
???
?
?
??
?
?
0
?
2
2???
?
?
?
2
+
???
?
?
??
?
?
0
?
2
2???
?
?
?
2
+
???
?
?
??
?
?
0
?
2
2???
?
?
?
2
+
???
?
?
?
?
??
?
?
?
?
0
?
2
2???
?
?
?
?
?
2
+
???
?
?
?
?
??
?
?
?
?
0
?
2
2???
?
?
?
?
?
2
+
???
?
?
?
?
??
?
?
?
?
0
?
2
2 ???
?
?
?
?
?
2
?exp ?
?1
2
?(?
?
?
)
?
? (?
?
+ ?
?
)
?1
? (?
?
?
)
?
?=1
?
4.6
where the trinormal model is indicated by "TN" and the components of the error tensor, ?
?
, are
calculated as before, with Equation 4.4. The individual optimal parameter and uncertainty values
are found by maximizing the likelihood function as described in Section 4.2.2.3. The results for
the volume element discussed here are summarized in Table 4.2. It is noted that the values for
?
?
?
, ?
?
?
, and ?
?
?
listed in Table 4.2 differ from those indicated in Figure 4.3 (a), this is due to
the fact that the standard deviation values shown in the figure were obtained through curve
fitting, which did not include a consideration of finite measurement error or the full
dimensionality of the velocity space.
86
Parameter Value
?
?
(8.3 ? 0.8) ? 10
9
m
?3
??? (0.00002 ? 0.00010)??
?
+ (0.00004 ? 0.00013)??
?
+ (0.00001 ? 0.00022)??
?
m/sec
?
?
?
0.01883 ? 0.0038 m/sec
?
?
?
0.02299 ? 0.0050 m/sec
?
?
?
0.02418 ? 0.0081 m/sec
?
?
?
?
?
0.627 ? 0.018
?
?
?
?
?
0.016 ? 0.076
?
?
?
?
?
0.192 ? 0.071
Table 4.2: PSD parameters with the trinormal velocity space model.
4.2.3: Comparison of the velocity space distribution models
To this point the PIV measurements have been modeled with both the Maxwellian and
the trinormal velocity space distributions. It can be seen, visually, in Figure 4.3, that the tri
normal model of the velocity space provides a better qualitative fit; but, as discussed in Section
2.3.2, almost any model will provide a better fit to a given data set than an alternate model if it
contains more free parameters. In this section the application of the posterior probability ratio
test outlined in Section 2.3.2 is applied to give a quantitative determination of whether the tri
normal model provides a superior description of the system, even with heavy penalization for the
extra parameters. Section 4.2.3.1 contains a discussion of the process used to calculate the ratio
and its value for the volume element analyzed above and Section 4.2.3.2 discusses the effects of
the large number of vectors that were recorded on the posterior probability ratio and calculated
parameter values.
87
4.2.3.1: The posterior probability ratio
The Bayesian approach to model selection was outlined in Section 2.3.2. Recall that the
ratio of the posterior probabilities for two models given the data set, D, and the accompanying
information, I, is given by:
? ?
prob (TN?,?)
prob (??,?)
4.7
where, in contrast to Equation 2.30, the models are identified as the trinormal model ("TN") and
Maxwellian model ("M") as opposed to A and B. The posterior probabilities are expanded by
use of Bayes' Theorem and a "fair" comparison of the models is made by setting the ratio of prior
probabilities for the models, given the information, to unity. After algebraic manipulation the
expression for K reduces to the ratio of likelihood functions:
? =
prob (?TN,?)
prob (??,?)
. 4.8
Next, the marginalization rule was applied to each of the likelihood functions for all
unknown parameters found in each model (as was discussed in detail within Section 2.3.3.1) and
uniform prior probability functions were chosen. All standard deviation values were assumed to
be within the range 0 < ? ? 0.10 ?/??? and the correlation values were assumed to have
values in the range ?1 < ? < 1. We begin with the likelihood function for the Maxwellian
model, marginalize over ?
?
, and assume normally distributed parameter uncertainty, ?
?
=
?
m0
? ??
?
:
88
prob (??,?) = ? prob (??
?
,?,?)prob (?
?
?,?)??
?
=
prob (??
m0
,?,?)
?
max
? exp?
?(???
m0
)
2
2(??
?
)
2
???
?
max
0
4.9
prob(??,?) =
?
?
2
??
?
0.1
?erf?
?
m0
?2??
?
?? erf ?
?
m0
? 0.1
?2??
?
??prob (??
m0
,?,?).
4.10
The final step needed to find the Maxwellian model's likelihood is to express the posterior
probability, prob(D?
?0
,M,I), in terms of the model and the measurement uncertainty, as was
outlined in Section 2.3.3.2:
prob(??
m0
,?,?) =
?(2?)
3 2?
?
??
?
+ ?
m0
2
?
?
??
??
exp?
?1
2
?(??
?
)
?
? ??
?
+ ?
m0
2
?
?
?
?1
? (??
?
)
?
?=1
?.
4.11
The process is the same for the trinormal velocity space distribution model. The
likelihood function is marginalized for the three standard deviation parameters and the three
correlation parameters, after integration:
prob(?TN,?) = ?
10?
4
?
3
???
?
?
??
?
?
??
?
?
??
?
?
?
?
??
?
?
?
?
??
?
?
?
?
? ?
?erf ?
?
?
?
0
?2??
?
?
?? erf ?
?
?
?
0
? 0.1
?2??
?
?
???erf ?
?
?
?
0
?2??
?
?
?? erf ?
?
?
?
0
? 0.1
?2??
?
?
?? ?
?erf ?
?
?
?
0
?2??
?
?
?? erf ?
?
?
?
0
? 0.1
?2??
?
?
???erf ?
1 ??
?
?
?
?
0
?2??
?
?
?
?
? + erf ?
1 + ?
?
?
?
?
0
?2??
?
?
?
?
?? ?
?erf ?
1 ??
?
?
?
?
0
?2??
?
?
?
?
? + erf ?
1 + ?
?
?
?
?
0
?2??
?
?
?
?
???erf ?
1 ??
?
?
?
?
0
?2??
?
?
?
?
? + erf ?
1 + ?
?
?
?
?
0
?2??
?
?
?
?
?? ?
prob (??
?
?
0
,?
?
?
0
,?
?
?
0
,?
?
?
?
?
0
,?
?
?
?
?
0
,?
?
?
?
?
0
,?,?)
4.12
89
The trinormal distribution likelihood function is given by combining the following:
prob ????
?
?
0
,?
?
?
0
,?
?
?
0
,?
?
?
?
?
0
,?
?
?
?
?
0
,?
?
?
?
?
0
, TN,?? =
?(2?)
3 2?
???
?
+ ?
?
0
??
??
?(2?)
3 2?
???
?
+ ?
?
0
??
??
exp ?
?1
2
?(??
?
)
?
? ??
?
+ ?
?
0
?
?1
? (??
?
)
?
?=1
?
?
?
0
= ?
?
?
?
0
2
?
?
?
?
?
0
?
?
?
0
?
?
?
0
?
?
?
?
?
0
?
?
?
0
?
?
?
0
?
?
?
?
?
0
?
?
?
0
?
?
?
0
?
?
?
0
2
?
?
?
?
?
0
?
?
?
0
?
?
?
0
?
?
?
?
?
0
?
?
?
0
?
?
?
0
?
?
?
?
?
0
?
?
?
0
?
?
?
0
?
?
?
0
2
?.
4.13
With these expressions for the likelihood functions, the ratio K can be expressed as the
product of a factor that penalizes for the distribution parameters, P, and the ratio of the likelihood
functions for the two models given the optimal parameter values:
? = ? ?
prob????
?
?
0
,?
?
?
0
,?
?
?
0
,?
?
?
?
?
0
,?
?
?
?
?
0
,?
?
?
?
?
0
, TN,??
prob(??
m0
,?,?)
? = ?
10?
4
?
3
??
?
?
??
?
?
??
?
?
??
?
?
?
?
??
?
?
?
?
??
?
?
?
?
?50???
?
?erf ?
?
?
?
0
?2??
?
?
?? erf ?
?
?
?
0
? 0.1
?2??
?
?
?? ?
?erf ?
?
?
?
0
?2??
?
?
?? erf ?
?
?
?
0
? 0.1
?2??
?
?
???erf ?
?
?
?
0
?2??
?
?
?? erf ?
?
?
?
0
? 0.1
?2??
?
?
?? ?
?erf ?
1 ??
?
?
?
?
0
?2??
?
?
?
?
? + erf ?
1 + ?
?
?
?
?
0
?2??
?
?
?
?
???erf ?
1 ??
?
?
?
?
0
?2??
?
?
?
?
? + erf ?
1 + ?
?
?
?
?
0
?2??
?
?
?
?
?? ?
?erf ?
1 ??
?
?
?
?
0
?2??
?
?
?
?
? + erf ?
1 + ?
?
?
?
?
0
?2??
?
?
?
?
???erf ?
1 ??
?
?
?
?
0
?2??
?
?
?
?
? + erf ?
1 + ?
?
?
?
?
0
?2??
?
?
?
?
?? ?
?erf (
?
m0
?2??
?
) ? erf (
?
m0
? 0.1
?2??
?
)?
?1
4.14
With Equations 4.11, 4.13, and 4.14 and the parameter values found in Tables 4.1 and 4.2 we can
make a firm quantitative statement about which of the velocity space models best describes the
data. For the volume element discussed in this section the value of the ratio was found to be
90
? = 2.9 ? 10
163
, which includes a penalty factor of ? = 7.3 ? 10
?9
. This value of K is much
higher than the threshold for a definitive conclusion that the trinormal model should be selected
(K>100). The penalty factor is on the order of 10
9
, meaning that the likelihood of the trinormal
model must be more than 10
11
larger than for the Maxwellian model to arrive at such a strong
conclusion (here the trinormal model's likelihood is ~10
172
times greater than that of the
Maxwellian model). The very large value of the K ratio will be shown to be systematic for all of
the volume elements in the dust cloud, but before proceeding, the astonishingly large magnitude
of this K value will be examined.
4.2.3.2: The effect of large data sets
The large value of the posterior probability ratio that was found for the volume element
discussed above seems, at first glance, to be almost too large, even when one considers the fact
that the trinormal model provides a much better qualitative fit to the data (as seen in Figure 4.3).
The large K value quoted in the previous section is a direct result of the large number of vectors
that were measured in the experiment and used throughout the analysis process. The effect of
the large number of measurements is most apparent in Equations 4.3 and 4.6, each likelihood
function has a term similar to:
?(2?)
3 2? ?
?
?
+ ?
m0
2
?
?
?
??
4.15
The standard deviation values are less than one, so these terms grow very quickly as the number
of measurements, N, becomes large. The large number of measurements also eliminates the
effects of small sample size, namely the effects of outlier data points and uncertainty in
91
Figure 4.5: Plots showing the order of magnitude of the K ratio as a function of the number of vectors randomly
selected from D. (a) log
10
K for between 10 and 1000 vectors. (b) A closer view of the boxed region of (a). This
view shows that the two models cannot be easily distinguished when a small number of measurements are made.
The black line is at K=10
0
, where the two models are equally as likely and the blue line is at K=100, where the tri
normal model can be said to be definitively preferred.
parameter estimates. The large data sets mean that strong conclusions about the relative
descriptive power of the distribution models can be drawn without ambiguity.
To test the effect of the number of vector measurements on the value of the posterior
probability ratio and the parameter value calculations, the data set for the volume element
discussed above was broken up into smaller pieces and reanalyzed. Randomly selected subsets
of D containing 10, 20, 30, 50, 75, 100, 200, 300, 500, 750, and 1000 vectors were selected and
analyzed (400 times for each subset size) to give the parameter values and K ratios as a function
of N using the same analysis process described in Section 4.2.2. The base ten logarithm of the
resulting K ratios is shown in Figure 4.5 (a), as a function of the number of vectors in the
randomly selected subsets. The posterior probability ratio's order of magnitude is directly
proportional to the number of vectors measured. Figure 4.5 (b) shows a closer view of the same
data for small values of N. It is noted that the majority of the PIV studies performed in the past
have had N?100 (most commonly with N?50), meaning that, because of the penalization for the
extra parameters, the value of posterior probability ratio in such cases would not indicate that the
92
trinormal model is a better description of the system. The large number of vectors measured
here allows us to definitively discriminate between the models.
Figure 4.6: Plots showing the parameter and uncertainty values obtained from random samples of D as a function of
the number of vectors within the sample. The black dots indicate the average value and the red bars indicate the
standard deviation of the parameter or uncertainty value.
93
The second effect of the large data set is that parameter value estimate uncertainty
decreases as more data is added. In Figure 4.6 the parameter value estimates and the associated
parameter value uncertainties are plotted as a function of N (black dots). The corresponding
standard deviation values of the estimates are also shown as a function of the number of data
points (the red bars). As a general rule, it can be seen that the optimal parameter values reach a
set value with fewer vectors (N~250) than for the uncertainties (N~500). Small parameter
uncertainty is an important factor in the value of the K ratio; if the parameter uncertainty values
are large they can dominate the expression for K through the parameter penalty factor (Equation
4.14). Finally, the more exact parameter values obtained by the inclusion of many vector
measurements in D means that even small differences between the models are magnified, which
makes the model comparison and selection process very clean.
4.2.4: Summary
The process of constructing the PSD within a single volume element of the dust cloud is
straightforward in principal, although somewhat mathematically involved in practice. The
configuration space component of phase space is modeled as a scalar quantity, the number
density, which is obtained through analysis of the camera images acquired with the PIV
diagnostic system. The velocity space was then modeled with both the canonical Maxwellian
distribution function and the trinormal probability distribution function. The optimal parameter
values for the two velocity space distribution models were obtained through the same general
process with the help of Bayesian probability theory. After the parameters for the two velocity
space distribution models were found the Bayesian model comparison/selection process was
employed and it was found that the trinormal model was approximately 163 orders of magnitude
more probable than the Maxwellian model.
94
4.3: Velocity space model selection
The PIV data was analyzed with the process described in Section 4.2 for all 853 volume
elements found within the dust cloud. Section 4.3.1 contains a discussion of the results of the
posterior probability ratio tests that were performed on all volume elements. Section 4.3.2
describes a geometrical method that can be used to visualize the extent of the anisotropy semi
quantitatively.
4.3.1: Posterior probability ratios
Figure 4.7: Histograms of the base 10 logarithm of the K ratio for all volume elements within the dust cloud. (a)
Histogram showing the results for all of the volume elements. (b) Closer view of the same data for small values of
log
10
K. The blue line indicates K=100, no volume elements within the dust cloud were found to have K values
below this cutoff.
The ratio of the posterior probability of trinormal model of the velocity space
measurements compared to that of the Maxwellian model was calculated, as in Section 4.2.3, for
all of the volume elements within the dust cloud. As was seen for the single volume element
described above, the values of K were found to be much greater than 100, indicating that the tri
normal model provides a much better description of the data across the entire dust cloud volume.
The base 10 logarithm of the K values are shown as histograms in Figures 4.7, part (a) shows
that the majority of the volume elements have K ratios in the 10
100
10
300
range. Figure 4.7 (b)
shows the same data but for smaller values of log
10
(K). The smallest value of the posterior
95
Figure 4.8: Histogram of the order of magnitude of the penalty factor assessed to the distributions for their unknown
parameters and uncertainties. The blue dashed line is a guide to the eye.
probability ratio for any volume element within the dust cloud was found to be K=7.3?10
3
,
which is still more than an order of magnitude above the threshold for a definitive conclusion
that the trinormal model provides a superior description of the data (indicated by the blue
dashed line in 4.7 (b)).
The process used to calculate the ratio K included a draconian penalty factor (P, Equation
4.14) for the extra parameters found in the trinormal model of the velocity space. The order of
magnitude of the penalty factor values for all volume elements is shown in Figure 4.8. The
values of the penalty factor were all found to be much less than one, indicating that the extra
parameters in the trinormal model would have to provide a significantly better description of the
observations in order to overcome the penalty and result in large K ratio values. The fact that
these penalty factor values favor the Maxwellian model so strongly makes the K ratio values
seen in Figure 4.7 even more compelling.
96
4.3.2: Geometrical comparison
The posterior probability ratio method of model comparison provides a concrete, but
somewhat abstract, quantitative metric which has allowed us to reach the conclusion that the tri
normal model is required. In this section several semiquantitative methods of comparing the
velocity space models will be given in an attempt to show the discrepancy between the models
more clearly. The first comparison, found in Section 4.3.2.1, looks at the volume of velocity
space that is occupied by each of the two models. The second comparison, the ratio of the
squares of the principal axis ellipsoid lengths to the square of the Maxwellian sphere radius, is
found in Section 4.3.2.2. Section 4.3.2.3 provides a geometrical explanation of the results by
combining the volume and axis length comparisons.
4.3.2.1: Velocity space volume
The volume of velocity space occupied by each of the distributions is the inverse of the
normalization constant used to make the distribution integrate to unity. The Maxwellian
probability distribution function is given by Equation 2.5, from which the volume can be seen to
be:
Vol
?
= (2?)
3 2?
?
?
3
. 4.16
Similarly, the volume occupied by the trinormal distribution is given by Equation 2.16:
Vol
TN
= (2?)
3 2?
??
?
?
1 2?
= (2?)
3 2?
?
?
1
?
?
2
?
?
3
.
4.17
These formulae for the velocity space volume are directly proportional to the volume occupied
by the spheres and ellipsoids discussed in Section 2.2. The volume of velocity space occupied
by the Maxwellian distribution function is equal to the volume of a sphere with radius ?
?
times a
97
Figure 4.9: Scatter plots comparing the velocity space volumes calculated with both models for each volume
element in the dust cloud. (a) All of the volume elements are shown. The solid black line indicates slope one. The
red dashed line indicates a line fitted to the data. (b) A closer view of the boxed region of (a). The blue dot
indicates the volume element discussed in Section 4.2. The green bars are representative error bars.
constant (3?? 2? ). The volume occupied by the trinormal distribution is proportional to the
product of the three ellipsoid axis lengths, which is equal to the volume of an ellipsoid with axes
lengths ?
?
1
, ?
?
2
, and ?
?
3
times the same constant as for the Maxwellian volume. The product of
the ellipsoid axis lengths is equal to the square root of the determinant of the velocity space
covariance matrix. By comparing the volumes of the distributions across the cloud structure we
can see that the two distributions describe very different regions of velocity space.
Figure 4.9 shows a scatter plot of the velocity space volume calculated with each model
for all of the locations in the cloud; the horizontal axis shows the volume calculated from the
Maxwellian model and the vertical axis indicates the volume found assuming the trinormal
model. Figure 4.9 (b) shows a closer view of the region in Figure 4.9 (a) enclosed by the black
rectangle. The Figure shows that the ratio of the volumes lies, more or less, on a line of slope
0.71, meaning that, on average, the Maxwellian model of the velocity space overestimates the
velocity space volume by approximately 29%.
98
Figure 4.10: Scatter plots showing the variance of the velocity space for the two models. The blue dots indicate ?
1
2
,
the red dots show ?
2
2
, and ?
3
2
is displayed as black. The lines indicate linear fits to the three sets of points. (a) Plot
showing the entire parameter space. (b) Plot showing the region indicated by the black box in (a). The error bars
are representative.
4.3.2.2: Variance ratios
The second geometrical comparison looks at the variance values of the trinormal model
and the Maxwellian distribution scalar variance. The scatter plot in Figure 4.10 shows the
Maxwellian model's variance (?
?
2
) on the horizontal axis and the three principal axis variance
values from the trinormal model on the vertical axis (?
1
2
is blue, ?
2
2
is red, and ?
3
2
is black). If
the velocity space were spherically symmetric the three points from each volume element would
have the same value and would all lie on a line of slope 1. The data were fit to a line, the slope
values (ratios) were found to be 1.77, 0.90, and 0.33, for the principal axis 1, 2, and 3 points,
respectively. The large circles in Figure 4.10 (b) indicate the values for the volume element
discussed in Section 4.2. The fact that the ratios of the trinormal distribution variance values are
not the same as the Maxwellian distribution variance is a clear reflection of the fact that the
velocity space of the dust component is highly anisotropic. The error bars shown in Figure 4.10
(b) are representative. It is noted that the uncertainty in the ?
1
2
values is smaller than for ?
2
2
or
?
3
2
because the algorithm that is used to find the eigenvalues calculates ?
1
2
first; the uncertainty
in the calculation of ?
1
2
is propagated into the results for the other variance values.
99
4.3.2.3: Geometrical connection
The velocity space volume ratio and variance anisotropy discussed above can be
combined to give a more complete representation of the difference between the two models. As
noted above, the Maxwellian model overestimates the volume of velocity space occupied by the
dust component. The geometrical picture of these distributions as ellipsoids and spheres
provides an explanation of this discrepancy. The ratio of the velocity space volume occupied by
the trinormal and the Maxwellian models is the ratio of Equations 4.17 and 4.16,:
? =
Vol
TN
Vol
?
=
?
?
1
?
?
2
?
?
3
?
?
3
4.18
When combined with the knowledge that the variance of the Maxwellian distribution is the
average value of the trinormal distribution's variance values the ratio can be written in terms of
two dimensionless axis length values, a and b:
? ? ?
?
1
?
?
?
? ? ?
?
2
?
?
?
? = ??
?
3 ??
2
??
2
4.19
The ratio R is plotted in Figure 4.11 (a), the black curves are at constant values of the volume
ratio (the vertical axis) and are spaced by 0.15 (the largest of which is at 0.90). The same
function and contours are shown in Figures 4.11 (b) and (c). The maximum value of the volume
ratio is one, which occurs when a = b = 1 (i.e. in the limit where the ellipsoidal surface becomes
a sphere). Geometrically, this means that if one has some length L to divide up between the axes
of an ellipsoid, the ellipsoid with the largest volume will be obtained when L is distributed
amongst the three axes evenly. This means that the Maxwellian distribution describes the
velocity space which has the largest possible volume, given a set average velocity variance (in
100
Figure 4.11: Plots showing the ratio of the volume of an ellipsoid to the volume of a sphere whose radius is the
average value of the length of the three ellipsoid axes. (a) Plots of the volume ratio value (vertical, colored) as a
function of the lengths of the two longest ellipsoid axis lengths, normalized to the spherical radius (?a? and ?b?).
The red dot (towards the back of the structure) indicates the location of the velocity space for the volume element
discussed in Section 4.2 in this parameter space. (b) A two dimensional plot of the same function as (a). The
additional black points indicate the location of the velocity spaces for all volume elements from the measurements.
(c) A closer view of (b).
Chapter 5 we will see that the average velocity variance is directly proportional to the energy
density; the two models give identical energy densities within some configuration space volume,
by having the same average axis length, but occupy very different volumes of velocity space).
The black points indicated in Figures 4.11 (b) and (c) show the data from the different volume
elements in the dust cloud in this normalized space (the red dot represents the volume element
101
discussed in Section 4.2). The bulk of the data lie in a group near the second black contour
which is at 75% of the maximum volume ratio, this is in general agreement with the slope of
0.71 from the volume ratio comparison, seen in Figure 4.9. This type of visualization is
convenient because it allows one to see the degree to which the velocity space is anisotropic,
semiquantitatively, in terms of three dimensionless quantities for each volume element (a, b, and
R).
4.4: The dust component phase space distribution
The spatially resolved phase space distribution has now been measured for the dust
component and the process of finding the various model parameters has been described. In
Section 4.3 we arrived at the conclusion that the velocity space component of phase space must
be modeled with the trinormal probability distribution function. In this section the parameters
required to specify the PSD (the number density scalar field, drift velocity vector field, and
velocity space covariance tensor field) are shown as a function of position for the entire dust
cloud structure. These parameter fields have not previously been measured in such a system and
show the phase space of the dust cloud in unprecedented detail. Additionally, a semi
quantitative understanding of the magnitude and spatial variation of the PSD parameter fields is
required before proceeding to the discussion of the fluid and transport quantities in Chapter 5.
4.4.1: The number density
The process used to find the dust number density, as discussed for a single volume
element in Section 4.2.1, was repeated for all of the volume elements across the dust cloud in
which the dust component was observed. The spatial distribution of the number density is shown
in Figure 4.12. The histogram at the top of the figure shows the number of volume elements
102
Figure 4.12: The dust number density distribution. The histogram at the top of the figure shows all of the volume
elements across the cloud, the colors indicate the magnitude of the number density for the plots of constant ?? cross
sections of the dust cloud. The frame at the top left (x=52.91 cm) shows the number density of the cross section as
viewed from beneath the anode looking towards the service end of 3DPX (similar to the vantage point of the
photograph in Figure 3.2 (b)). The other cross sections shown, moving left to right (and top to bottom) in the figure,
are viewed from the same vantage point below the anode moving through the dust cloud away from the anode. The
side of the 3DPX chamber with the large window through which the PIV images are recorded is on the right side of
the cross section plots.
with a given number density and indicates the color scale used in the plots of the scalar fields
shown below the histogram. The plots of the spatially resolved dust number density in Figure
4.12 show planes of constant ??, as viewed from beneath the anode looking away from the PIV
103
Figure 4.13: Schematic view illustrating the orientation used for figures of scalar field quantities. The cuboids
show the volume elements occupied by the dust component of the plasma and are colored by number density (as in
Figure 4.12). (a) View showing the location of the dust cloud in relation to the anode (the brown disk) and the dust
tray (the gray surface at the bottom). This view is equivalent to looking through the large window on the
experimental section of 3DPX (see Figure 3.1 (a)). (b) Closer view of the structure from the same vantage point as
part (a). The translucent gray, vertically oriented, planes show and identify the absolute and relative locations of
several easily identifiable cross sections found in all plots of scalar field quantities within this dissertation.
laser towards the service end of 3DPX. The relative locations of the planes seen in Figure 4.12
are shown in Figure 4.13. Figure 4.13 (a) shows a view of the dust cloud as seen by an observer
looking through the large window in 3DPX (as in Figure 3.1 (a)). Figure 4.13 (b) shows the
volume elements in the cloud (colored by number density as in Figure 4.12 (a)) from the same
vantage point as in 4.13 (a) only closer, also indicated are the locations of several easily
identifiable cross sections used in the plots of scalar field quantities. All scalar field quantities
within this dissertation show the same cross sections of the cloud that are indicated here.
Returning to Figure 4.12, it can be seen that the number density varies fairly continuously
across the dust cloud structure. The cloud features a region of slightly higher number density
104
near it's center (and slightly to the right, at higher ?? values) which decreases towards the outer
boundary of the cloud. Although the densities across the structure are all within, approximately,
an order of magnitude of one another the variation is clearly nonzero and has a regular structure.
It is also noted that although this dust cloud is quite large compared to other clouds observed
within 3DPX, the total volume of the dust cloud (~5 ? 10
?6
?
3
) is still less than one percent of
the experimental volume that is accessible with the PIV diagnostic (~9 ? 10
?4
?
3
). Thus, if we
were to make the standard assumption that the dust component is infinite and homogenous across
the entire volume the resulting number density would be on the order of 10
7
m
3
, which clearly
demonstrates the importance the retaining information about the spatial distribution.
4.4.2: The drift velocity
The drift (mean) velocity vector of the dust component was also calculated within all
volume elements of the cloud structure. Figure 4.14 (a) shows the view of the structure from the
same point that was shown in Figure 4.13 (a) and indicates the views for the vector field plots in
parts (c), (e), and (f), the vantage point for the view of the vector field in part (d) is the same as
in (a), except from a point that is closer to the cloud. The detailed structure of the vector field
can be difficult to see, viewing the vector field from several different perspectives allows the
information to be seen more clearly. All plots of vector field quantities within this dissertation
are from the four vantage points seen in Figure 4.14 (c) through (f) and appear in the same order
in all figures.
105
Figure 4.14: Spatial distribution of the drift velocity vectors. (a) A view from the top of the PIV camera window,
showing the drift vector directions and the cloud orientation with respect to the anode (shown in brown at the top)
and dust tray (the gray surface at the bottom of the plot). The view in part (a) is the same as that in Figure 4.13. The
figure also indicates the vantage points for the plots in parts (c), (e), and (f). (b) Histogram of the magnitude of the
drift vectors. The histogram shows the distribution of drift velocity vector magnitude and indicates the color scale
of the vectors shown below. The plots of the vectors show the magnitude and orientation of the drift velocity
vectors measured within the dust cloud. (c) View from the anode. (d) View from the same vantage point as in (a).
(e) View from below the anode, near the dust tray surface. (f) View from the bottom of the PIV camera window.
106
Figure 4.15: Spatial distribution of the drift velocity vector magnitude. The histogram at the top of the figure shows
the magnitudes of the drift velocity vectors and the color scale used in the plots below. The frames below the
histogram show the magnitude of the drift velocity as a function of position, with the same color scheme as the
vectors plotted in Figure 4.13. The vantage point is the same as in Figure 4.12.
The drift velocity vectors are colored according to magnitude, as indicated by the
histogram in Figure 4.14 (b). The field can be seen to vary in both magnitude and direction in all
three spatial dimensions. The variation of the drift velocity magnitude can be seen more clearly
in Figure 4.15, which shows the magnitude as a scalar field on the same planes of constant ?? as
in Figure 4.12. By taking both of the figures into account it can be seen that the drift velocity
field has a, generally, larger magnitude at lower ?? values (towards the "back" of the chamber, or
107
Figure 4.16: Spatial distribution of the velocity space standard deviation values along the ??
?
direction. The
histogram shown at the top of the figure indicates the color scale used in the plots of constant ?? which are shown
below.
on the left side of the scalar field plots) and at higher ?? values (at points further away from the
anode). It can also be seen that in the regions where the magnitude of the drift velocity is highest
the vector field has a negative ?? (vertical) component, which will be important in Chapter 5.
4.4.3: The velocity space covariance tensor, laboratory coordinate system
The final component required for the specification of the PSD is the velocity covariance
tensor, ?
?
. In the general laboratory coordinate system there are six parameters in the covariance
108
Figure 4.17: Spatial distribution of the velocity space standard deviation values along the ??
?
direction. The
histogram shown at the top of the figure indicates the color scale used in the plots of constant ?? which are shown
below.
tensor (Equation 2.14): Three standard deviation values and three correlation factors. The
process for calculating these parameters within a single volume element was described in Section
4.2.2.5. Figures 4.16, 4.17, and 4.18 show the spatial distribution of the standard deviation
values ?
?
?
, ?
?
?
, and ?
?
?
, respectively, for planes of constant ??. In the figures it can be seen that
?
?
?
and ?
?
?
have approximately the same spatial distribution, while the distribution of ?
?
?
shows
much more spatial variation and a broader distribution of values. In all three cases the
109
Figure 4.18: Spatial distribution of the velocity space standard deviation values along the ??
?
direction. The
histogram shown at the top of the figure indicates the color scale used in the plots of constant ?? which are shown
below.
parameters can be seen to have lower values in the center of the cloud structure and higher
values towards the edges of cloud. In the case of ?
?
?
the higher standard deviation values are
found towards the back of the cloud (at smaller ?? values) and at higher ?? values.
4.4.4: The velocity space covariance tensor, principal axis coordinate system
The correlation values were not shown in the previous section because without the
corresponding standard deviation values they are of limited utility. The information gained with
110
Figure 4.19: Spatial distribution of the velocity space standard deviation values along the ??
1
direction (the largest of
the standard deviation values in the principal axis coordinate system). The histogram shown at the top of the figure
indicates the color scale used in the plots of constant ?? which are shown below.
the inclusion of the correlation is seen most easily by rotating the velocity space coordinate
system in each volume element to the local principal axis basis set. The details of the rotation
process are discussed in Section 2.2.3; briefly, the rotation is to the coordinate system where the
covariance tensor is diagonal. In the principal axis system the shape of the velocity space is
completely specified by the three standard deviation values along the principal axis basis set
directions, ?
?
1
, ?
?
2
, and ?
?
3
. The spatial distributions of the principal axis standard deviation
111
Figure 4.20: Spatial distribution of the velocity space standard deviation values along the ??
2
direction (the second
largest of the standard deviation values in the principal axis coordinate system). The histogram shown at the top of
the figure indicates the color scale used in the plots of constant ?? which are shown below.
values are plotted as scalar fields in Figures 4.19, 4.20, and 4.21. The variation of the standard
deviation values is apparent in all three plots, as opposed to the case of the laboratory coordinate
system, where the variation was almost entirely confined to ?
?
?
. The standard deviation scalar
field plots show the magnitude and the spatial variation of the velocity space width.
To gain a complete view of the information stored in the covariance tensor field the
magnitudes of the standard deviation values must be accompanied by their corresponding
112
Figure 4.21: Spatial distribution of the velocity space standard deviation values along the ??
3
direction (the smallest
of the standard deviation values in the principal axis coordinate system). The histogram shown at the top of the
figure indicates the color scale used in the plots of constant ?? which are shown below.
orientation, in the form of the principal axis basis set directions. Each tensor element comes with
two unit vectors (there is no ?positive? or ?negative? direction, the ellipsoids are
morphologically identical under complete inversion of the Cartesian basis set) which means that,
in principal, each of the lines that are plotted should have an arrow on both ends, these are
omitted here for visual simplicity. These coordinate system axial directions are plotted in
Figures 4.22, 4.23, and 4.24. Figure 4.22 shows the spatial variation of the ??
1
orientation as a
function of position across the dust cloud structure. In the plots of ??
1
the color scale indicates
113
Figure 4.22: Views of the null
1
vector directions throughout the cloud volume. The color scheme is indicated in (a)
and shows the angle the axial direction makes with the vertical. (b) View from the anode. (c) View from the same
vantage point as in (a). (d) View from below the anode, near the dust tray surface. (e) View from the bottom of the
PIV camera window.
the angle between the first principal axis and the vertical, nullnull, direction, as shown in the histogram
in Figure 4.22 (a). The orientation of the first principal axis has a regular structure while
114
Figure 4.23: Views of the null
2
vector directions throughout the cloud volume. The color scheme is indicated in (a)
and shows the angle between the null
2
axis and the null? direction. (b) View from the anode. (c) View from the same
vantage point as in (a). (d) View from below the anode, near the dust tray surface. (e) View from the bottom of the
PIV camera window.
featuring nonzero spatial variation. It is noted that the first principal axis is along the direction
that exhibits the largest velocity space standard deviation; if one were to simply look at the plots
of the velocity space standard deviation values in the laboratory coordinate system (Figures 4.16,
115
Figure 4.24: Views of the null
null
vector directions throughout the cloud volume. The color scheme is indicated in (a)
and shows the angle between the null
3
axis and the null? direction. (b) View from the anode. (c) View from the same
vantage point as in (a). (d) View from below the anode, near the dust tray surface. (e) View from the bottom of the
PIV camera window.
4.17, and 4.18) the expectation would be that this direction would be parallel to null?. However,
examination of Figure 4.22 indicates that the direction of largest velocity variance is generally
116
not in the ?? direction, the inclusion of the offdiagonal terms in the covariance tensor is the only
way in which this direction can be identified.
Figure 4.23 shows the spatial orientation and variation of the second principal axis, ??
2
.
The color shows the angle between the orientation of ??
2
and the ?? axis (for an easier comparison
with ??
3
). The majority of the cloud features a nominally uniform ??
2
direction, except in the
region colored blue. In this region ??
2
is aligned, more or less, with the ?? direction. The very
qualitative explanation for the deviation of the direction of ??
2
orientation in this region is that the
region with blue coloration in Figure 4.23 is the only part of the cloud directly below the anode.
It can also be seen, in Figure 4.20, that this region features lower values of ?
?
2
than the rest of
the cloud. The orientation of the third principal axis is shown in Figure 4.24, the color scheme is
the same as in the plots of the second principal axis direction. Examination of the ??
3
orientation
shows the same interesting structure seen in the plots of ??
2
.
The myriad plots in this section are an attempt to show the variation of a tensor field with
six independent parameters on a three dimensional spatial array; the large amount of information
that must be specified in each volume element precludes a simple, compact, representation. By
combining the standard deviation magnitude scalar fields (Figures 4.19, 4.20 and 4.21) with the
orientation of the velocity space coordinate systems (Figures 4.22  4.24) one can gain a basic
qualitative mental picture of the spatially varying velocity space covariance tensor field.
117
Chapter 5: Transport and Thermal Properties of the Dust Component
In Chapter four the process of measuring and constructing the spatially resolved phase
space distribution for the dust component of a plasma was described in detail. The process led to
three field quantities (the number density scalar field, the drift velocity vector field, and the
velocity space covariance tensor field) that, together, describe the properties of the dust fluid as a
function of position. The parameter fields, shown in Section 4.4, are quite interesting in their
own right, but the true power of the diagnostic, data analysis, and velocity space model
improvements introduced in this dissertation is the ability to take the information for the
individual volume elements and use it to ascertain some of the more general thermodynamic and
transport properties of the dust component. The new measurement capacity gives the ability to
see the internal structure of the dust component with unprecedented detail.
A number of fluid quantities that have previously been inaccessible are now available
with the PSD measurement process described in Chapter 4. Section 5.1 discusses the derivation
of the fluid transport equations starting from the Boltzmann equation and with the assumption
that the velocity space is trinormal. The results of the derivations are three transport equations:
The continuity equation, the momentum equation, and the energy equation. The terms that
appear within each of the three transport equations (and have been measured) will be discussed,
semiquantitatively, in Sections 5.2 ? 5.4. Examination of the spatial distribution and magnitude
of the various densities and flux rates that appear in the transport equations allows one to see the
influence of each individual thermodynamic and transport mechanisms that were measured.
118
Section 5.5 contains further discussion of the results from earlier in the chapter. Specifically,
Section 5.5.1 highlights the fact that the dust component is in a state of dynamic forcebalanced
equilibrium. Section 5.5.2 discusses the difference in the functional form of the thermal heat
flux vector that one obtains when using the trinormal velocity space distribution, as opposed to
the canonical Maxwellian velocity space distribution.
5.1: The fluid transport equations
The fluid equations that govern the thermodynamic and transport properties of the dust
fluid come from the Boltzmann equation
45
. The standard textbook derivations of the transport
equations
33,46?48
are tied at their core to the assumption that the velocity space can be described
with the spherically symmetric Maxwellian distribution function (with perturbations, in the most
general cases). We have seen that the standard velocity space distribution model does a
relatively poor job of describing the system in question, so we must rederive the transport
equations from the beginning. It is noted that in the analysis given here we are interested in the
quantities for which the parameters are experimentally measured, current technological and
diagnostic limitations preclude detailed knowledge of the ion, electron, and neutral population
distributions; because of this experimental inaccessibility all collision terms are omitted.
Additionally, the electric field structure and dust grain charge distribution are both unknown in
the experimental system, leading to the omission of terms related to the electrostatic force that
would normally appear in such derivations. The ignored terms are, presumably, important
mechanisms in the system but because any analysis that includes such terms would be
overwhelmingly dominated by their uncertainties the omission is justifiable. All of the analysis
discussed here is based on the PSD measurements described in the preceding chapters, as a result
of the wealth of information obtained through the PSD measurements we can learn a great deal
119
Integral form Result
? f
d
d
3
v ?
?
? v
?
f
d
d
3
v n
d
u
i
??
?
2
f
d
d
3
v
n
d
(u
i
2
+ ?
v
i
2
)
? v
?
v
?
f
d
d
3
v
n
d
(u
i
u
j
+ ?
v
i
v
j
?
v
i
?
v
j
)
??
?
3
f
d
d
3
v n
d
u
i
(u
i
2
+ 3?
v
i
2
)
??
?
2
?
?
f
d
d
3
v
n
d
(2 u
i
?
v
i
v
j
?
v
i
?
v
j
+ u
j
(?
v
i
2
+ ?
v
j
2
))
? v
?
v
?
?
?
f
d
d
3
v
n
d
(u
i
u
j
u
k
+ u
k
?
v
i
v
j
?
v
i
?
v
j
+ u
j
?
v
i
v
k
?
v
i
?
v
k
+ u
i
?
v
j
v
k
?
v
j
?
v
k
)
Table 5.1: Velocity space integral results with the trinormal velocity distribution function.
about the dust component even while we neglect several of the important mechanisms acting on
the system.
The first three velocity space moments of the Boltzmann equation give rise to equations
that govern density, momentum, and energy transport balances in the plasma fluid. With the
omission of collisions the Boltzmann equation reduces to the Vlasov equation:
0 =
??
?
??
+ v?? ? ?
?
?
?
+
?
?
?
?
? ?
?
?
?
5.1
where f
d
is the dust component PSD, ?
?
is the vector describing any external forces that act on the
dust fluid, and where ?
?
and ?
?
are the gradient operators in configuration and velocity space,
respectively. In what follows the force vector is the sum of the gravitational force and an
omnibus "other" force, ?
?
?????
, which act on the dust component:
?
?
= ??
?
? ?? + ?
?
?????
5.2
where ? ? +9.8 ?/???
2
. To assist with some of the velocity space integration that follows
several particularly useful results are given in Table 5.1.
120
The continuity equation is obtained through the zeroth velocity space moment of the
Vlasov equation: We multiply Equation 5.1 by ??
0
and integrate over the infinite three
dimensional velocity space (all integration limits over velocity space are omitted in this section
but should understood to cover the full infinite range):
0 = ?
??
?
??
?
3
? + ? v?? ? ?
?
?
?
?
3
? + ?
?
?
?
?
? ?
?
?
?
?
3
?
5.3
The derivation procedure for the continuity equation follows the standard recipe, the detailed
steps of which can be found in any plasma physics textbook. The result is the continuity
equation:
??
?
??
= ??
?
? (?
?
???).
5.4
The momentum equation comes from the first velocity space moment of the Vlasov
equation. In order to arrive at an equation that describes the rate of change of momentum density
we multiply Equation 5.1 by ?
?
?? and integrate:
0 = ? (?
?
v??)
??
?
??
?
3
? + ? (?
?
v??)v?? ? ?
?
?
?
?
3
? + ? (?
?
v??)
?
?
?
?
? ?
?
?
?
?
3
?
5.5
The first term on the right hand side ("RHS") of Equation 5.5 gives the standard expression for
the time rate of change for the fluid's momentum density:
? (?
?
v??)
??
?
??
?
3
? =
?
??
(?
?
?
?
???)
5.6
The second term on the RHS of 5.5 is simplified by noting that the spatial gradient of v?? is zero,
this allows the term to be rewritten, with the help of trivial vector identities, as:
? (?
?
v??)v?? ? ?
?
?
?
?
3
? = ?
?
? ??
?
? v??? v???
?
?
3
?? 5.7
121
Where ????? indicates the outer product of ?? with itself, resulting in a rank two tensor quantity.
The integrals can be carried out as in Table 5.1 or, directly, with Equation 2.6:
? (?
?
??)?? ? ?
?
?
?
?
3
? = ?
?
? (2?
?
+ ?
?
?
?
??? ????)
where the pressure tensor is defined as:
?
?
=
1
2
?
?
?
?
?
?
5.9
The remaining term in Equation 5.5 is simplified using the standard application of integration by
parts, the result is:
? (?
?
v??)
F
??
?
?
? ?
?
?
?
?
3
? = ?
?
?
?
? y? ??
?
F
??
other
5.10
Equations 5.6, 5.8, and 5.10 are then reinserted into Equation 5.5, slight algebraic manipulation
gives the momentum equation:
?
??
(?
?
?
?
???) = ?2?
?
? ?
?
??
?
? (?
?
?
?
??? ????) ??
?
?
?
? y? + ?
?
?
?
other
5.11
The final fluid equation of interest is the energy equation. There are two forms of the
energy equation, the scalar and the rank two tensor versions. The scalar energy equation is
obtained by multiplying the Equation 5.1 by the scalar quantity:
1
2
?
?
?? ? ?? =
1
2
?
?
(?
?
2
+ ?
?
2
+ ?
?
2
) =
1
2
?
?
?
2
and integrating over velocity space. The tensor form is found by multiplying by the tensor
quantity:
122
1
2
?
?
?? ??? =
1
2
?
?
?
?
?
2
?
?
?
?
?
?
?
?
?
?
?
?
?
?
2
?
?
?
?
?
?
?
?
?
?
?
?
?
?
2
?
We will restrict our attention to the scalar version. The energy transport equation is then:
0 = ? ?
1
2
?
?
?? ? ???
??
?
??
?
3
? + ? ?
1
2
?
?
?? ? ????? ? ?
?
?
?
?
3
? +
? ?
1
2
?
?
?? ? ???
?
?
?
?
? ?
?
?
?
?
3
?
5.12
The first term on the right hand side is very similar to the integration in Equation 5.8, the result
is:
? ?
1
2
?
?
?? ? ???
??
?
??
?
3
? =
?
??
?Tr ??
?
?? +
?
??
?
1
2
?
?
?
?
?
2
? 5.13
Where Tr(?
?
) is the trace of the pressure tensor and ?
2
= ?
?
2
+ ?
?
2
+ ?
?
2
(which is the trace of the
tensor formed by the outer product of the drift velocity with itself). The second term on the RHS
of 5.12 is where we see the biggest difference with the new velocity distribution model, the terms
that appear here also have the largest notational variation in the literature, so we proceed with
caution:
? ?
1
2
?
?
?? ? ????? ? ?
?
?
?
?
3
? =
1
2
?
?
? ?
(?
?
2
+ ?
?
2
+ ?
?
2
)?
?
(?
?
2
+ ?
?
2
+ ?
?
2
)?
?
(?
?
2
+ ?
?
2
+ ?
?
2
)?
?
? ? ?
?
?
?
?
3
?
The spatial derivative operator can be pulled out of the integral, as in 5.7:
? ?
1
2
?
?
?? ? ????? ? ?
?
?
?
?
3
? =
1
2
?
?
?
?
? ? ?
(?
?
2
+ ?
?
2
+ ?
?
2
)?
?
(?
?
2
+ ?
?
2
+ ?
?
2
)?
?
(?
?
2
+ ?
?
2
+ ?
?
2
)?
?
??
?
?
3
?
123
The nine integrals found in the equation above can all be carried out with the results listed in
Table 5.1, giving:
? ?
1
2
?
?
?? ? ????? ? ?
?
?
?
?
3
? =
1
2
?
?
?
?
?
?
?
?
?
?
?
?3?
?
2
+ ?
?
2
+ ?
?
2
? + ?
?
?
xy
?
?
?
?
+ ?
?
?
xz
?
?
?
?
+ ?
?
?
2
?
?
?
xy
?
?
?
?
+ ?
?
??
?
2
+ 3?
?
2
+ ?
?
2
? + ?
?
?
yz
?
?
?
?
+ ?
?
?
2
?
?
?
xz
?
?
?
?
+ ?
?
?
yz
?
?
?
?
+ ?
?
??
?
2
+ ?
?
2
+ 3?
?
2
? + ?
?
?
2
??
?
?
?
?
=
1
2
?
?
?
?
? ??
?
??
2
??? + ??? ? ?
?
+ ??? Tr ??
?
? + t
?
? ???? ??
?
???
5.14
Where "?" indicates a triple inner product and t
?
is a rank four tensor, defined as:
t
?
? ?
??
?
??
?
??
?
^
?
?
^
?
?
^
?
?
^
?
5.15
Where ? is the Kronecker delta. Equation 5.14 is further simplified by introducing the thermal
heat flux vector, ?
??
:
?
??
?
1
2
?
?
?
?
???? ? ?
?
+ ??? Tr??
?
? + t
?
? ???? ??
?
?? 5.16
It is noted that this somewhat awkward notation is avoided in the tensor version of the energy
equation, where the expression for the rank 3 tensor quantity for the thermal heat flux is:
?
?
?
1
2
?
?
?
?
???? ??
?
+ ?
?
???? + ??
?
?????
?
?
The terms in the parenthesis of the above equation are simply the sum of all unique rank three
combinations of the covariance tensor and the drift velocity vector. With Equation 5.16, the
second term in the energy equation becomes:
124
? ?
1
2
?
?
?? ? ????? ? ?
?
?
?
?
3
? = ?
?
? ?
1
2
?
?
?
?
?
2
???? + ?
?
? ?
??
5.17
Which is the sum of the divergence of the thermal heat flux vector and the divergence of a
quantity that describes a flux of heat which is entirely born from the drifting motion of the dust
fluid. The final term on the right hand side of Equation 5.12 describes the energy flow due to
external forces:
? ?
1
2
?
?
?? ? ???
F
??
?
?
? ?
?
?
?
?
3
? =
1
2
?
?
? ???
?
2
+ ?
?
2
+ ?
?
2
??
?
?
?
?
3
?
5.18
This is simplified by rewriting the integrand with the product rule of partial differentiation:
? ?
1
2
?
?
?? ? ???
F
??
?
?
? ?
?
?
?
?
3
? =
1
2
?
?
? ? ??
?
(?
2
?
?
) ??
?
?
??
2
??
?
?
??
2
??
?
?
??
2
??
?
?
???
3
?
The first term of the integrand on the RHS is integrated by parts to give zero. After the
derivatives in the second term are taken integration yields:
? ?
1
2
?
?
?? ? ???
F
??
?
?
? ?
?
?
?
?
3
? = ??
?
??? ? F
??
= ?
?
?
?
? ?
?
??
?
??? ? F
??
other
5.19
Thus, with the results found in Equations 5.13, 5.17, and 5.19, the energy equation is given by:
?
??
?Tr??
?
?? +
?
??
?
1
2
?
?
?
?
?
2
? =
??
?
? ?
1
2
?
?
?
?
?
2
??????
?
? Q
???
??
?
?
?
??
?
+ ?
?
??? ? F
??
other
5.20
The first three moments of the Boltzmann equation have given three equations that
govern our simplified description of the dust fluid. Each equation has several terms which are
both of interest and whose components have been measured. In the discussion that follows the
125
Continuity Equation
?n
d
?t
= ??
r
? (n
d
???)
Momentum Equation
?
?t
(m
d
n
d
???) = ?2?
r
? P
?
??
r
? (n
d
m
d
??? ????) ? m
d
n
d
g ?? + n
d
?
?
other
Energy Equation
?
?t
Tr ?P
?
? +
?
?t
?
1
2
m
d
n
d
u
2
? = ??
r
? ?
1
2
m
d
n
d
u
2
??????
r
? ?
??
? m
d
n
d
g u
y
+ n
d
??? ? ?
?
other
Pressure Tensor
P
?
=
1
2
n
d
m
d
?
?
???? ???? ?????? ?
??
?
1
2
m
d
n
d
??????
?
+??? Tr??
?
? + t
?
? ???? ? ?
?
??
Table 5.2: Transport equations summary.
spatial distribution and variation of the terms that appear in these equations will be examined.
For convenience, the important results from this section are summarized in Table 5.2.
5.2: The continuity equation
In this section terms that appear in the continuity equation are examined. The quantities
of interest are found on the right hand side of Equation 5.4 and are listed, for convenience, in
Table 5.3.
The term found within the divergence operator of the continuity equation is most easily
examined when multiplied by the dust mass, to give units of momentum density:
??
?
= ?
?
?
?
??? 5.21
The momentum density vector field can be seen in Figures 5.1 (see Section 4.4.2 for an
explanation of the different views of the vector field). The color, arrow head size, and length of
the vectors shown in the figure are on the logarithmic scale indicated in part (a) of the figure, this
scaling is used to allow the structure of the momentum density field to be seen more clearly. It
can be seen in the plots of the vector field (Figures 5.1 (b) through (d)) that the momentum
density has a generally larger magnitude towards the "back" of the chamber (low ?? values) and at
higher values of ??. The vector field is parallel to the drift velocity field, which is shown in
126
?n
d
?t
= ??
r
? (n
d
???)
Quantity Expression Units
Momentum Density ?
?
?
?
??? (kg m/sec) m
 3
Particle Flux Rate ?
?
?(?
?
???) sec
?1
?
?3
Table 5.3: Continuity equation quantities and corresponding SI units.
Figure 4.14, but is more uniform in magnitude due to the fact that the regions of the cloud with
higher number density (Figure 4.12) tend to have lower drift velocity.
The divergence of the momentum density vector field is much easier to visualize, as it is
a scalar field. The finite difference method is used to calculate the spatial derivative quantities
found in this chapter, details can be found in Appendix A. This term is particularly interesting if
we take the result of the divergence calculation and multiply by the configuration space volume
of each volume element. The multiplication gives the rate of particle flux out of each element
per unit time, meaning that positive values of this quantity correspond to a loss of particles. The
particle flux scalar field is shown in Figure 5.2 (see Section 4.4.1 for a description of the relative
locations for the various cross sections that are shown), the histogram shows the color scale used
in the scalar field plots found beneath the histogram. The region with the largest particle flux is
localized near the back of the cloud, roughly corresponding to the region with the highest
momentum density magnitude, as seen in Figure 5.1.
127
Figure 5.1: The momentum density vector field. (a) Logarithmic histogram showing the magnitudes of the
momentum density vectors. The color scale used in the histogram indicates the coloration of the vectors in (b)
through (e).
128
Figure 5.2: Particle flux rate through each volume element in the dust cloud structure. (a) Histogram showing the
particle flux rates calculated throughout the dust cloud structure. The color scale in (a) indicates the magnitude of
particle flux rate as a function of position in the scalar field plots below. Positive values shown in this plot indicate
a divergent flux rate.
129
?
?t
(m
d
n
d
???) = ?2 ?
r
? P
?
??
r
? (n
d
m
d
??? ????) ? m
d
n
d
g ?? + n
d
?
?
other
Quantity Expression Units
Drift Energy Density
1
2
?
?
?
?
??? ? ???
J/m
3
Thermal Energy Density
P
?
=
1
2
n
d
m
d
Tr(?
?
)
J/m
3
Gravitational Potential Energy Density m
d
n
d
g y J/m
3
Drift Force Density ?
r
? ?
1
2
n
d
m
d
??? ????? J/m
3
Thermal Force Density
?
?
? ?
?
J/m
3
Gravitational Force Density m
d
n
d
g ?? J/m
3
Table 5.4: Quantities found within the momentum equation and the corresponding SI units.
5.3: The momentum equation
The momentum equation is equivalent to Newton's second law; it relates the time rate of
the change in the dust component's momentum density to any forces acting on the system. There
are two types of quantities found within this equation whose components have been measured:
Energy densities and their spatial derivatives, which give force densities. A semiquantitative
picture of the energy density distribution is quite useful because it gives insight into the general
spatial thermodynamic structure of the cloud and allows one to see the widely disparate
magnitudes of energy density provided by the different terms ("mechanisms") in the momentum
equation. The distributions of the force density vector fields, and their magnitudes, provide
similarly useful information.
5.3.1: Energy densities
The distribution of the energy density associated with the drift velocity is shown in
Figure 5.3. The color scale shown in the histogram indicates the magnitude of the energy density
130
Figure 5.3: Drift velocity based energy density. (a) Histogram showing the energy density values calculated
throughout the cloud attributable to the drift velocity. The energy densities are shown on a logarithmic scale. The
average approximate thermal energy densities of the ion, electron, and neutral plasma populations are indicated on
the histogram for scale. The scalar field plots show the spatial distribution of the energy density using the same
color scheme that is used in the histogram.
in the spatial distribution plots. The approximate thermal energy densities for the neutral (U
n
),
electron (null
null
null), and ion (null
nullnull
null) components are shown in the histogram (using the values listed in
Table 3.1) to give context to the magnitudes of the measured drift velocity energy densities. It is
noted that because the three types of dust component energy density discussed in this section
have such disparate magnitude scales the ion, electron, and neutral thermal energy densities that
appear in the Figures are scaled by a factor which is indicated in the plots (for example, the
thermal energy density of the electrons is null
null
null null7.5null10
nullnull
nullnull/null
null
, this is many orders of
131
magnitude larger than the scale in Figure 5.3, so this value for ?
?
? is multiplied by 10
?8
in
order to have it appear on the histogram). The scale used in the histogram is logarithmic, which
allows easy comparison of the energy density magnitudes provided with the different
mechanisms described in this section. The general spatial structure of the drift velocitybased
energy density is similar to that of the momentum density distribution, the regions on the left of
the plots (low z values, the "back" of the chamber) have higher energy densities than in the
portion of the cloud towards the front of the chamber. The drift energy density distribution can
also be seen to have a larger magnitude and higher uniformity for the larger values of x. The
area with the lowest drift energy density roughly corresponds to the region where the second
principal axis direction is aligned with the ?? direction (the blue region of Figure 4.23).
A quantity of great interest in the dusty plasma physics community is the density of the
thermal energy. The tensor variance included within the trinormal velocity distribution function
required a slight generalization of the standard definition for the pressure tensor, ?
?
, given in
Equation 5.9. The generalization is fairly straightforward; the use of separate ion population
pressure values parallel and perpendicular to an external magnetic field uses the same logic,
although the method presented here is much more general. It is easy to see that the thermal
energy density, defined here as the trace of the pressure tensor, corresponds to the canonical
thermal energy density, ?
th
= 3?
?
?
?
? 2? , when the trinormal velocity covariance tensor is
isotropic (i.e. when the variance is equal to the square of the thermal velocity, ?
?
2
? ?
th
2
=
?
?
? ?
?
? ):
132
Figure 5.4: Thermal energy density. (a) Histogram showing the energy density values calculated throughout the
cloud attributable to the width of the velocity distribution. The energy densities are shown on a logarithmic scale
and cover a higher range of values than in the plots of the drift velocity based energy density. The average
approximate thermal energy densities of the ion, electron, and neutral plasma populations are indicated on the
histogram for scale. The scalar field plots show the spatial distribution of the energy density using the same color
scheme that was used in the histogram.
Tr nullnull
?
nullnullnull
nullnull
1
2
null
null
null
null
Trnullnull
nullnull
null
null
?
null
Tr nullnull
?
nullnullnull
nullnull
3
2
null
null
null
null
null
nullnull
null
Tr nullnull
?
nullnullnull
nullnull
3
2
null
null
null
null
null
Trnullnull
?
nullnullnull
nullnullnull
nullnull
133
The thermal energy density is an especially convenient metric in this context due to the
ambiguity that one encounters when applying the idea of a temperature to systems where the tri
normal model is required to accurately describe the velocity space. It can be seen in the figures
within this section (Figures 5.3 ? 5.5) that the various measured energy densities of the dust
component all have values that are orders of magnitude less than or equal to the approximate
thermal energy densities of the ion, electron, and neutral populations (except in the case of the
gravitational potential energy of the dust fluid, which is has an average energy density
approximately four orders of magnitude larger than the ion thermal energy density). This type of
information is normally discussed within the context of a dust temperature, but the application of
a temperature based energy density metric to a fluid with a velocity space distribution that is not
the spherically symmetric Maxwellian may not be applicable. Coupled with the ambiguity
surrounding the actual applicability of a ?temperature? is the fact that the number density of the
dust is many orders of magnitude less than that of the other plasma components. Discussing
values of thermal energy density, instead of temperatures, seems to be a more appropriate
measure of energy content.
The thermal energy density distribution measured for the cloud is shown in Figure 5.4.
The color scale and approximate values for the other plasma components are, again, indicated in
the histogram to give an idea of the scale of the dust component thermal energy density. The
numerical range of the energy density used in the figure is narrower, but several orders of
magnitude higher, than those that were measured for the drift energy density.
The third type of energy density which is accessible through the PSD measurements is
the gravitational potential energy density (with zero potential energy defined to be on the surface
134
Figure 5.5: Gravitational potential energy density. (a) Histogram showing the gravitational potential energy
density values calculated throughout the cloud. The energy densities are shown on a logarithmic scale and cover a
higher range of values than in the plots of the drift velocity or velocity distribution width based energy densities.
The average approximate thermal energy densities of the ion, electron, and neutral plasma populations are indicated
on the histogram for scale. The scalar field plots show the spatial distribution of the energy density using the same
color scheme that was used in the histogram.
of the dust tray), the spatial distribution is shown in Figure 5.5. The scale for these energy
density values is the highest of the three types of energy density. As one would expect, the
gravitational potential energy is generally highest near the top of the cloud, except near the edges
where the number density is low.
The general structure of the spatial distributions of the various energy densities is quite
interesting by itself, but perhaps the more important application of this type of analysis is the
135
Figure 5.6: Comparison of energy density scales. The horizontal axis shows the gravitational potential energy
density calculated for all volume elements in the dust cloud. The vertical axis indicates the drift velocity based
energy density as black points and the thermal energy density as red points. The data are plotted logarithmically
along both axes. The dashed lines indicate where the data would lie if the thermal/drift energy densities were equal
to the gravitational potential energy multiplied by the indicated power of ten. The blue and black points found
above the black dashed line show where the approximate locations of the ion and neutral populations in this
parameter space.
ability to compare the magnitudes of the three measured energy density mechanisms. Figure 5.6
shows the thermal (red dots) and drift (black dots) energy density values for each volume
element in the dust cloud plotted on the vertical axis and shows the gravitational potential energy
density plotted on the horizontal. In the Figure it is clear that, of the three energy sources
considered, the gravitational potential energy mechanism dominates, due to the high mass of the
dust grains. The values of the thermal energy density are approximately two orders of magnitude
smaller than the gravitational potential energy, the energy density associated with the drift
appears to be almost negligible. The approximate values for the argon neutrals (the black dot
labeled Ar
n
) and the argon ions (the blue dot labeled Ar
+
) are also shown, the electrons are
omitted because their gravitational potential energy density is approximately four orders of
magnitude lower than minimum value that appears on the plot. The average ratios of the thermal
and drift energies for the dust component and the ratio of the thermal energy to the gravitational
136
Plasma Species Ratio Value
Dust Drift/Gravity 9.3 ? 10
?7
Dust Thermal/Gravity 1.6 ? 10
?2
Neutral Thermal/Gravity 8.2 ? 10
2
Ion Thermal/Gravity 1.2 ? 10
3
Electron Thermal/Gravity 1.3 ? 10
11
Table 5.5: Average energy density ratios.
potential energy of the ion, electron, and neutral components of the plasma are listed in Table
5.5.
5.3.2: Force densities
The terms that actually appear in the momentum transport equation are the spatial
derivatives of the energy density fields discussed in Section 5.3.1. The force densities are plotted
as both three dimensional vector fields and as scalar fields to show the magnitude. The
magnitudes of the force densities will be shown to have a hierarchy that is similar to that which
was seen for the energy densities.
The force that is entirely associated with the drift velocity is shown as a vector field in
Figure 5.7. The magnitude of the vectors and the projection of the force density vector direction
on the plane of the figure is shown in Figure 5.8, the arrows on the scalar field plots are all of the
same length, regardless of magnitude. The spatial distribution of the drift force density is similar
to that of the corresponding energy, with larger values towards the back of the chamber. The
magnitudes of the force density range from approximately 10
10
to 10
7
N/m
3
.
137
Figure 5.7: Drift force density vector field. (a) Histogram showing the magnitude of the force density vectors
calculated throughout the cloud attributable to the drift velocity divergence. The forces densities are shown on a
logarithmic scale. The vector field plots show the spatial distribution of the drift force density using the same color
scheme that was used in the histogram.
138
Figure 5.8: Drift force density magnitude. (a) Histogram showing the magnitude of the force density vectors
calculated throughout the cloud attributable to the drift velocity divergence. The magnitudes are shown on a
logarithmic scale. The scalar field plots show the spatial distribution of the drift force density magnitude using the
same color scheme that was used in the histogram. The arrows superimposed on the scalar field plots are all of the
same length and are shown to indicate the direction of the vector field in the plane of the cross section. The arrows
colored red have a positive horizontal component and the black arrows have a negative horizontal component.
The divergence of the pressure tensor (the "thermal force density") is shown in Figures
5.9 and 5.10. The magnitude of the thermal force densities ranges from approximately 10
5
to
10
4
N/m
3
, approximately three orders of magnitude higher than for the drift force densities. A
particularly interesting feature of the spatial distribution can be seen by inspecting the arrows on
the plots of the scalar field. At low x values the null?null
?
vector field has nullnull (vertical) components
that are all pointing upward, as one moves towards larger values of x the magnitude of
139
Figure 5.9: Thermal force density vector field. (a) Histogram showing the magnitude of the force density vectors
calculated throughout the cloud attributable to the divergence of the velocity space distribution width and
orientation. The force densities are shown on a logarithmic scale. The vector field plots show the spatial
distribution of the thermal force density using the same color scheme that was used in the histogram.
140
Figure 5.10: Thermal force density magnitude. (a) Histogram showing the magnitude of the force density vectors
calculated throughout the cloud. The force density magnitude is shown on a logarithmic scale. The scalar field
plots show the spatial distribution of the thermal force density magnitude using the same color scheme that was used
in the histogram. The arrows superimposed on the scalar field plots are all of the same length and are shown to
indicate the direction of the force vector field in the plane of the cross section. The arrows colored red have a
positive vertical component (they point upward) and the black arrows have a negative vertical component.
thermal force density increases and the orientation changes, becoming predominantly downward.
This behavior can also be seen in the plots of the vector field, particularly in Figures 5.9 (c) and
(e).
The final force density field which can be directly computed from the PSD measurements
is the gravitational force density. The vector plots are omitted due to the fact that the field is
141
Figure 5.11: Gravitational force density magnitude. (a) Histogram showing the magnitude of the gravitational force
density vectors calculated throughout the cloud. The magnitude is shown on a logarithmic scale. The scalar field
plots show the spatial distribution of the gravitational force density magnitude using the same color scheme that was
used in the histogram. Arrows are omitted in this figure, as the vectors are uniformly oriented downward.
uniform in direction (down). The magnitude of the gravitational force density has a fairly
narrow distribution and large values compared to the drift and thermal force densities.
A comparison of the magnitudes of the force densities acting on each volume element is
shown in Figure 5.12. The thermal and drift force density magnitudes are plotted vertically and
the magnitude of the gravitational force density is shown on the horizontal. In contrast to the
comparison of the energy density magnitudes, the magnitudes of the three types of force are
closer to one another. In several cases the thermal force density is larger than the gravitational
142
Figure 5.12: Comparison of force density magnitudes. The horizontal axis shows the magnitude of the gravitational
force density calculated for all volume elements in the dust cloud. The vertical axis indicates the magnitude of the
drift velocity based force densities as black points and the thermal force densities as red points. The data are plotted
logarithmically along both axes. The dashed lines indicate where the data would lie if the thermal/drift force density
magnitudes were equal to the gravitational force density multiplied by the indicated power of ten.
force density and for almost all of the volume elements the magnitudes of the thermal and
gravitational force densities are within three orders of magnitude of one another. As was seen
for the energies, the contribution from the drift is extremely low and seems to play a rather small
role in the overall force balance.
5.4: The energy equation
The scalar energy equation gives information about the rate at which energy flows
through each volume element in the dust cloud. In this section the derivative terms will be
examined, plots of the actual heat flux vector fields are omitted. The quantities shown in this
section are the actual values of the energy flowing into each volume element and not densities
(as with the particle flux rate in Section 5.2).
143
?
?t
Tr?P
?
?+
?
? t
?
1
2
m
d
n
d
u
2
?=??
r
??
1
2
m
d
n
d
u
2
??????
r
??
??
?m
d
n
d
g u
y
+n
d
??? ? ?
?
other
Quantity Expression Units
Drift Energy Flow Rate Density ?
r
??
1
2
m
d
n
d
u
2
???? (J/sec)m
3
Thermal Energy Flow Rate Density
?
r
??
??
(J/sec)m
3
Gravitational Potential Energy Flow Rate Density m
d
n
d
g u
y
(J/sec)m
3
Table 5.6: Quantities found within the energy transport equation and the corresponding SI units.
Figure 5.13 shows the divergence of the flux of heat that is a result of the drifting motion
of the dust fluid. Positive values of the energy flow rate indicate that the "drift kinetic energy" of
the fluid decreases as it flows through the volume element (a positive divergence value). The
magnitude of the drift energy flow rates are relatively low and are highly peaked about zero. The
spatial distribution indicates that the region with the most negative driftmediated energy flow is
near the back edge of the cloud, where the fluid gains energy as a result of the change in drift.
The region with the largest positive values for the drift energy flow rate nominally
correspond to the interface of the blue region of the cloud in Figure 5.3, where the drift energy
density is lowest, and the remainder of the cloud where the drift energy density is higher and
more uniform.
The divergence of the thermal heat flux vector, ?
??
, is shown in Figure 5.14. It is noted
that the scale used in the figure is approximately four orders of magnitude wider than was used
for the drift energy flow rate plots. Positive values of the thermal energy flow rate are generally
found in the interior of the cloud. The region of the cloud at low z values (on the left of the
scalar field cross section plots) can be seen to have a concentration of negative values, this
indicates that the energy transport provided by pressure differences is flowing into these volume
144
Figure 5.13: Drift velocity based kinetic energy flow rate. (a) Histogram showing the rate of energy flow due to the
drift velocity based mechanism. Positive values indicate a divergence (loss) of energy for a given volume element.
The scalar field plots use the same color scheme as in (a) and show that the majority of the energy transport due to
this mechanism is confined to the rear of the cloud (i.e. at low z value, on the left hand side of the cross sections).
elements (i.e. a negative divergence, or a convergence, of heat flux) as a result of the thermal
heat flux mechanism.
The rate of gravitational potential energy flow is shown in Figure 5.15. The distribution
of this quantity is highly asymmetric about zero and is skewed towards negative values, which
indicates that, for this cloud, gravitational potential energy flows out of the blue colored regions
in Figure 5.15 due to the presence of the drift velocity vector field, which points down in the blue
145
Figure 5.14: Thermal energy flow rate. (a) Histogram showing the rate of energy flow due to the divergence of the
heat flux vector. Positive values indicate a divergence (loss) of energy for a given volume element. The scalar field
plots use the same color scheme as in (a) and show that the energy transport due to this mechanism is distributed
throughout the cloud volume.
region. The reason for this can be seen by examination of the drift velocity vector field (Figure
4.14) and its magnitude (Figure 4.15): The magnitude of the y component of the drift velocity
field becomes larger, but remains negative, as one moves from the low ?? values to higher ??
values. The regions with larger vertical drift velocity magnitude can be seen to correspond to the
regions with a higher gravitational potential energy redistribution rate. As a dust grain drifts to a
region that is situated lower in the cloud it loses gravitational potential energy, which must be re
146
Figure 5.15: Gravitational potential energy flow rate. (a) Histogram showing the rate of gravitational potential
energy flow. The negative (blue) values indicate a convergence (gain) of energy for the particles within the volume
element due to the falling motion caused by the drift velocity vector field. The scalar field plots use the same color
scheme as in (a) and show that the majority of the energy transport due to this mechanism is confined to the side of
the cloud furthest from the anode (i.e. at high x values).
deposited into some other form of energy (most likely in the form of electrostatic potential
energy in this experimental configuration).
5.5: Further discussion of the transport measurements
The process of examining each of the individual measurable terms that appear in the
transport equations is now complete. The process has given a great deal of insight into the inner
workings of the system. Section 5.5.1 describes how the evidence we have considered to this
147
point allows one to reach the conclusion that the system is in a state of dynamic forcebalanced
equilibrium. Section 5.5.2 discusses the thermal heat flux quantity and describes how the two
models of velocity space result in different functional forms of the heat flux vector. Section
5.5.3 contains a summary.
5.5.1: A dynamic equilibrium
The motivation for the first conclusion that we can draw from the transport analysis
comes from consideration of the terms within the energy equation. To see the combined effect
of the measured energy flow mechanisms in the system we return to the energy equation, 5.20:
?
??
Tr??
?
? +
?
??
?
1
2
?
?
?
?
?
2
? =
??
?
? (
1
2
?
?
?
?
?
2
???) ??
?
? ?
??
??
?
?
?
gu
?
+ ?
?
??? ? ?
?
other
5.22
The total measured time rate of change of energy content in each volume element is the sum of
the terms on the right hand side of the equation (we exclude the energy flow associated with the
"other" forces, as they are, by definition, unknown):
??
?
? (
1
2
?
?
?
?
?
2
???) ??
?
? ?
??
??
?
?
?
g u
?
This quantity is plotted in Figure 5.16. It is noted that the scale in Figure 5.16 is the reverse of
that used in the preceding plots of the energy flow rate. If it were incorrectly assumed that there
were no other sources of potential energy storage in the system, a positive value would indicate
that the trace of the pressure tensor and the drift velocity based energy density, both found on the
left hand side of Equation 5.22, would increase as a function of time. An alternate way to think
148
Figure 5.16: Total measured energy flow rate. (a) Histogram showing the rate of gravitational potential energy
flow. The scale used in figure is the opposite of the previous plots; positive values indicate a gain of energy by the
particles within the dust fluid as they pass through a given region of the cloud. The scalar field plots use the same
color scheme as in (a) and show that the majority of the measurable energy transport is confined to the rear of the
cloud (i.e. at low z value, to the left hand side of the cross sections).
about this is that if the energy density content within a given volume element is assumed to be
constant in time (if the left hand side of Equation 5.22 is zero) the values of the energy flow
indicated in Figure 5.16 are those that are provided by the forces in the system which are
experimentally inaccessible.
Comparison of the scalar fields shown in Figure 5.16 to those in Figure 5.15 indicates
that gravitational potential energy flow is the dominant measurable energy flow mechanism. It is
149
also important to note that the gravitational mechanism is nonzero only if the drift velocity is
included. The preponderance of evidence to this point has, perhaps, suggested that the drifting
motion of the dust fluid is an inconsequential component of the overall system dynamics. The
drift component of the energy density field for the system was an average of nearly 10
8
times
smaller than the gravitational potential energy density. The force arising from the divergence of
the drift energy density tensor field was some six orders of magnitude smaller than the
gravitational force. The energy flow rates that arose purely from the drifting motion (Figure
5.13) were found to be extremely low and peaked about zero. However, if the drift velocity was
naively assumed to be zero, the gravitationally mediated energy redistribution mechanism would
also be zero, as would all other terms on the right hand side of the energy equation. While the
nonzero drift velocity field plays only a small role in the first two transport equations it is the
mechanism that drives the flow of energy in the system.
The observation of an energy flow that is the result of a nonzero and nonuniform drift
velocity field is the characteristic sign of a system that is in a state of dynamic equilibrium. The
standard practice of ignoring the drift velocity in the analysis of such systems results in the
omission of important physics. The drift velocity field cannot be observed without the new
method that allows the phase space distribution to be measured with spatial resolution.
5.5.2: The heat flux
The tensor variance in the trinormal model leads to a difference in the functional form of
thermal heat flux vector, ?
??
(this is also true if the rank two tensor version of the energy equation
is examined, in which case the heat flux quantities are rank three tensors). In Section 5.1 the
150
Figure 5.17: Difference in the energy flow rates calculated using the two velocity space models. (a) Histogram
showing the magnitude of the differences. (b) Histogram showing the energy flow rate magnitude as calculated
with the trinormal velocity space model (red) and the associated uncertainty in the values (gray) to give an idea of
scale to the other portions of the figure. The scalar field plots show the spatial distribution of the difference in
energy flow rates and show that there is minimal coherent spatial structure in the calculated differences.
vector form of ?
??
was derived for the trinormal velocity space model. The heat flux vector, in
terms of the trinormal velocity distribution function parameters, is given by:
?
??
TN
=
1
2
?
?
?
?
?
?
?
(3?
?
2
+ ?
?
2
+ ?
?
2
) + ?
?
?
xy
?
?
?
?
+ ?
?
?
xz
?
?
?
?
?
?
?
xy
?
?
?
?
+ ?
?
(?
?
2
+ 3?
?
2
+ ?
?
2
) + ?
?
?
yz
?
?
?
?
?
?
?
xz
?
?
?
?
+ ?
?
?
yz
?
?
?
?
+ ?
?
(?
?
2
+ ?
?
2
+ 3?
?
2
)
? 5.23
151
Figure 5.18: Comparison of the magnitudes and angles of the heat flux vectors calculated with the two velocity
space models. (a) Percentage difference (normalized to the trinormal model's heat flux vector magnitude) in the
calculated vector magnitudes. (b) Distribution of the angular orientation differences in the calculated vectors.
The result one obtains with use of the Maxwellian model, with the same energy density and drift
velocity, is given by:
?
??
?
=
5
2
?
?
?
?
?
?
2
??? =
5
6
?
?
?
?
(?
?
2
+ ?
?
2
+ ?
?
2
)???
5.24
The two expressions for the heat flux vector result in different vector fields, as can be seen in
Figure 5.17. Figure 5.17 shows the difference in the energy flow rates that are found when the
different velocity space models are used:
? (? ? ?
??
) ? ? ? ?
??
TN
?? ? ?
??
m
5.25
Figure 5.17 (b) shows the magnitude of the energy flow rate assuming the trinormal model, and
its uncertainty, for scale (see Figure 5.14 for the full spatial distribution of ? ? ?
??
TN
). The median
uncertainty of the energy flow rate values was found to be approximately 800 eV/sec, the
resolution for the histogram shown in 5.17 (a) is one third of this uncertainty. The histogram
resolution was chosen so that in the scalar field plots of ? (? ? ?
??
) a volume element that shows
an approximately white color has a calculated energy flow rate that is nominally the same for
both models, within the limits of the propagated uncertainties. The difference in the resulting
calculated energy flow rates can be seen to be pervasive throughout the cloud structure.
152
The actual thermal heat flux vectors obtained with the different velocity space models
also show differences. Two convenient metrics are used to show the difference in the heat flux
vector calculation results: The percentage difference of the magnitudes of the vectors and the
difference in the angle between the calculated vectors. Figure 5.18 (a) shows a histogram of the
percent difference in the magnitudes of the heat flux vectors obtained with each model and the
angle between the calculated vectors for each of the volume elements is shown in the histogram
in (b).
5.5.3: Summary
The fact that the heat flux vector field exhibits measurable differences when the more
general model of velocity space is applied underscores the importance of modeling the velocity
space data with the trinormal velocity space distribution. This is especially true when one
considers that the dust cloud discussed in this work was specifically chosen because of the fact
that it appeared to exhibit very little transport. The a priori expectation was that the data
analysis process would reveal little particle, momentum, or energy transport. The semi
quantitative observations discussed in this Chapter have shown, a posteriori, that these quantities
are not zero. Additionally, even in a case such as that described here, where the magnitude of
the drift velocity was very small, the standard assumption that the cloud is in a static, forcefree,
equilibrium was unambiguously shown to inadequately describe the state of the cloud.
Examination of the energy flow distribution measured within the system clearly showed that the
equilibrium is dynamic. The measured force density vector field showed that the gravitational
force, a force that arises externally, dominates over the contribution of the thermal and drift
forces that come from within the dust component itself; this indicates that the system is in a state
of forcebalanced equilibrium. This characterization of the dust fluid as being in a state of force
153
balanced equilibrium is selfconsistent with the observation of the trinormally distributed
velocity space. The ellipsoidal symmetry of the trinormal velocity space distribution was shown
to be the direct result of a system being in a forcebalanced equilibrium by ClerkMaxwell in
1867. The ability to resolve, measure, and quantify the internal transport properties, even in
cases such as this where the magnitude of the transport is relatively small, gives weight to the
power of the diagnostic and analytic processes that have been developed, applied, and discussed
herein.
154
Chapter 6: Conclusions
The work presented in this dissertation has addressed many of the issues related to the
characterization of weaklycoupled dusty plasma systems that had been unresolved, open,
questions. Broadly, this was done with a theoretical component (the application of the trinormal
model of the velocity space) and an experimental component (the development of spatially
resolved phase space distribution measurements). The combination of these developments has
led to the ability to directly measure and deduce the transport and thermal properties for the dust
component of such plasma systems. Section 6.1 briefly summarizes the work that has been
performed. The limitations of the components of this work will be discussed in Section 6.2,
which motivates a discussion of possible future work. Brief concluding remarks appear in
Section 6.3.
6.1: Summary
Perhaps the most broadly applicable of the topics that have been discussed is the
development of the trinormal probability distribution function as the model for the velocity
space portion of phase space. Section 2.2 gave an overview of both the canonical Maxwellian
distribution and the trinormal model, which contained a convenient method of geometric
visualization and discussed the applicability of the trinormal model in the context of both the
larger hierarchy of probability distribution functions and in a historical context. What was not
discussed, in detail, is how the trinormal distribution relates to that developed by Chew,
155
Goldberger, and Low
34
(?CGL?), which is a widely applied simplification of the trinormal
model. Briefly, in the method of CGL, and in subsequent slight generalizations, the appearance
of an anisotropic pressure tensor comes about by expanding the Boltzmann equation (with the
assumption that the velocity space is described by the isotropic Maxwellian distribution function)
in powers of some small parameter, generally the ion mass to charge ratio m
i
/q
i,
along a single
coordinate axis. The dust mass to charge ratio is much larger than one, meaning that such a
formulation is not applicable to dusty plasmas. The fact that the trinormal model can be used to
obtain an equivalent, and indeed more general, anisotropic pressure tensor without reliance on
such an expansion means that part of this work could certainly be of use beyond the realm of
dusty plasmas.
In the narrower context of weaklycoupled dusty plasmas, the application of the tri
normal model for the velocity space has addressed and answered several questions that were
previously open: Must the velocity distribution be spherically symmetric? Is the measured
velocity space anisotropy a real effect, or the product of uncertainty in the measurement? What
mechanism could produce such anisotropy? It was shown, by explaining that the trinormal
distribution is "allowed" in the context of phase space based Boltzmann equationtype analysis,
that the velocity space is not required to be spherically symmetric. The application of Bayesian
probability theory based statistical tests quantitatively showed that the anisotropy is real by
taking the finite measurement error of the diagnostic into account and also showed that the tri
normal distribution does a vastly better job of describing the observations. The mechanism that
brings about the anisotropy comes from the longknown fact that any three dimensional system
that is not in a forcefree equilibrium must be ellipsoidally symmetric in velocity space, as shown
by ClerkMaxwell in 1867
32
.
156
The second broad advancement was the development of a method that allows the phase
space distribution of a weaklycoupled dust cloud to be measured with high spatial resolution.
Previous analyses of weaklycoupled dusty plasma systems gave the following: A single
average number density, a single drift velocity vector (if this was not assumed zero), and a single
velocity space distribution width characterization. As seen in Section 4.4, these parameters vary
as a function of position throughout dust cloud structures. The spatial variation was evident
before the type of analysis presented here was given, by visual inspection of dust cloud
structures, but had not been previously explored rigorously. This work gives a method that
allows both the configuration and velocity space components of the PSD to be measured with
high spatial resolution.
The measurement of the transport and thermodynamic quantities discussed in this work
are new developments in the analysis of weaklycoupled dust plasma systems. Previous work
lacked both an adequate description of the velocity space and the spatial resolution that are
required to ascertain this type of information. Additionally, the systematic removal of the
assumption that the system is in a static, force free, equilibrium had not been attempted. The
analysis and derivation of the equation that describes the transport of energy also revealed that
the standard form of the heat flux vector is not applicable to systems such as the one described
here. Specifically, the following thermodynamic and transport quantities were previously
unavailable: The particle flux, the momentum density vector field, the energy density scalar (and
tensor) fields, the force vector fields, the heat flux vector (and tensor) fields, and the scalar (and
tensor) fields that describe the rate of energy flow. These quantities built a semiquantitative
picture of the internal structure of the cloud which allowed the following conclusion to be drawn:
The drift velocity plays an integral part in the system dynamics, because of this the cloud is in a
157
state of dynamic forcebalanced equilibrium. While many of these properties had been
suspected in weaklycoupled dusty plasmas, they had never been explored experimentally.
6.2: Future work
The work presented here opens the door to the precise analysis of any number of
phenomena observed within dusty plasmas. However, there are several aspects of the process
that could be investigated and improved. In terms of the information measured and discussed
here, the two areas ripest for improvement are: A better understanding of the configuration
space distribution of the dust grains and volumetric PIV measurements. In this work the
configuration space was simply modeled as a scalar field; one of the most commonly studied
dusty plasma systems is the stronglycoupled "crystalline" system. Such systems contain a high
level of dust grain spatial correlation, it may be beneficial to improve the model of the
configuration space to include such effects, especially since the transition between the strongly
and weaklycoupled regime is not well understood. The stereoPIV diagnostic has the limitation
that it can only measure a velocity field within a thin planar laser sheet, this means that data must
be acquired over long periods of time. After each cross section of the cloud is measured many
times, the PIV laser and cameras are moved and the process repeats itself. Dust clouds must be
present, and nominally stable, for long periods of time in order to allow the diagnostic process to
take place (it took approximately two and a half hours to record the data for the cloud described
above). By acquiring volumetric measurements of the cloud, with, perhaps, tomographic,
holographic, or plenoptic based velocimetry techniques, the velocity field of the entire cloud
volume could be measured simultaneously, allowing the investigation of more transient
phenomena.
158
Knowledge of spatially resolved background plasma species parameters and the electric
field structure within the experimental device together represent the most pressing area for
improvement. Average values for the background plasma parameters are known, at best, to
within an order of magnitude and without any awareness of spatial variation. If such information
were available the transport equations for the combined dust, ion, electron, and neutral system
could be used to ascertain very detailed properties relating to, for example, energy exchange
between the plasma species. The electric field structure is completely unknown in the device,
this is problematic because, after gravity, the electrostatic force is presumably the most important
mechanism that acts on the dust component.
6.3: Concluding remarks
This work has answered a number of the questions relating to the internal structure of
weaklycoupled dusty plasma systems. The application of the trinormal probability distribution
function as the model for the velocity space of the dust component allowed the observations to
be described adequately. The new velocity space model, coupled with the development of
spatially resolved measurements of the PSD, unambiguously showed, for the first time, that such
systems are in a state of dynamic forcebalanced equilibrium. The ability to draw such a
conclusion is especially promising for future work because the dust cloud structure discussed
here was chosen because it appeared to exhibit very little transport. The fact that the
improvements to the model of the velocity space and the diagnostic system allowed this
information to be ascertained in such a case speaks volumes to the sensitivity and precision that
is now possible in the study of these systems.
159
References
1
W. Crookes, in The British Association for the Advancement of Science (Elsevier, 1879).
2
W. Crookes, Proceedings of the Royal Society of London 30, 469472 (1879).
3
I. Langmuir, Science 60, 392394 (1924).
4
P. K. Shukla and A. A. Mamun, Introduction to Dusty Plasmas (IOP Publishing, LTD, Bristol,
2002).
5
V. Nosenko, J. Goree, and a. Piel, Physics of Plasmas 13, 032106 (2006).
6
F. F. Chen, Plasma Diagnostic Techniques (Academic Press, 1965).
7
J. G. Laframboise and L. W. Parker, Physics of Fluids 16, (1973).
8
J. D. Williams, Ph.D. Dissertation, Measurement of the Thermal Properties of a Weakly
Coupled Complex (Dusty) Plasma, Auburn University, 2006.
9
J. D. Williams and E. Thomas, Physics of Plasmas 13, 063509 (2006).
10
J. D. Williams and E. Thomas, Physics of Plasmas 14, 063702 (2007).
11
B. Liu, J. Goree, and Y. Feng, Physical Review E 78, 110 (2008).
12
B. Liu and J. Goree, Physical Review Letters 100, 14 (2008).
13
R. Fisher and E. Thomas, Physics of Plasmas 18, 113701 (2011).
14
R. Fisher and E. Thomas, IEEE Transactions on Plasma Science 38, 833837 (2010).
15
V. E. Fortov, O. Vaulina, O. F. Petrov, M. Vasiliev, A. Gavrikov, I. Shakova, N. Vorona, Y.
Khrustalyov, A. Manohin, and A. Chernyshev, Physical Review E 75, 18 (2007).
16
E. Thomas, R. Fisher, and R. L. Merlino, Physics of Plasmas 14, 123701 (2007).
17
R. A. Quinn and J. Goree, Physics of Plasmas 7, 3904 (2000).
18
I. Aranson and L. Tsimring, Reviews of Modern Physics 78, 641692 (2006).
160
19
O. Arp, D. Block, M. Klindworth, and a. Piel, Physics of Plasmas 12, 122102 (2005).
20
C. Boesse, M. Henry, T. Hyde, and L. Matthews, Advances in Space Research 34, 23742378
(2004).
21
a. V. Ivlev, H. M. Thomas, G. E. Morfill, V. I. Molotkov, A. M. Lipaev, V. E. Fortov, and T.
Hagl, New Journal of Physics 8, 25 (2006).
22
T. Miksch and A. Melzer, Physical Review E 75, 5 (2007).
23
I. Pilch, T. Reichstein, and A. Piel, Physics of Plasmas 15, 103706 (2008).
24
T. Reichstein, I. Pilch, and a. Piel, IEEE Transactions on Plasma Science 38, 814817 (2010).
25
J. ClerkMaxwell, Philosophical Transactions of the Royal Society of London. 157, 4988
(1867).
26
K. Pearson, Philosophical Transactions of the Royal Society of London. Series A. 186, 343
(1895).
27
K. Pearson, Philosophical Transactions of the Royal Society of London. Series A. 197, 443
(1901).
28
K. Pearson, Philosophical Transactions of the Royal Society of London. Series A. 216, 429
(1916).
29
J. Owen and R. Rabinovitch, The Journal of Finance 38, 745 (1983).
30
A. Azzalini, Scandinavian Journal of Statistics 32, 159188 (2005).
31
A. Azzalini, Biometrika 83, 715726 (1996).
32
J. ClerkMaxwell, The Quarterly Journal of Pure and Applied Mathematics 32, (1867).
33
S. I. Braginskii, Reviews of Plasma Physics 1, 205311 (1965).
34
G. F. Chew, M. L. Goldberger, and F. E. Low, Proceedings of the Royal Society of London.
Series A, Mathematical and Physical Sciences 234, 112 (1956).
35
L. Spitzer Jr, The Physics of Fully Ionized Gases (Interscience Publishers, New York, 1955).
36
Saul Stahl, Mathematics Magazine 79, 96 (2006).
37
D. S. Sivia and J. Skilling, Data Analysis: A Bayesian Tutorial, Second Ed. (Oxford
University Press, USA, New York, 2006).
161
38
H. Jeffreys, The Theory of Probability (3e) (Oxford, 1961), p. 432.
39
R. Lindken, M. Rossi, S. Grosse, and J. Westerweel, Lab on a Chip 9, 25512567 (2009).
40
J. Westerweel, Measurement Science and Technology 8, 13791392 (1997).
41
R. J. Adrian, Experiments in Fluids 39, 159169 (2005).
42
E. Thomas, J. D. Williams, and J. Silver, Physics of Plasmas 11, L37 (2004).
43
E. Thomas, J. D. Williams, and C. Rath, IEEE Transactions on Plasma Science 38, 892896
(2010).
44
E. Thomas and J. D. Williams, Physics of Plasmas 13, 055702 (2006).
45
L. Boltzmann, Sitzungsberichte Der Akademie Der Wissenschaften Wien II 66, 275 (1872).
46
F. F. Chen, Introduction to Plasma Physics and Controlled Fusion Volume 1: Plasma Physics,
2nd ed. (Springer Science, New York, 1974).
47
D. A. Gurnett and A. Bhattacharjee, Introduction to Plasma Physics With Space and
Laboratory Applications (Cambridge University Press, Cambridge, 2005).
48
H. Goedbloed and S. Poedts, Principles of Magnetohydrodynamics (Cambridge University
Press, New York, 2004).
162
Appendix A: Finite difference derivatives
The distribution functions, ? (??,??), discussed in this work is measured on a discrete
spatial array. The form of the gradient operator that is used, primarily in Chapter 5, is described
here. The spatial array has regular spacing in the three configuration space dimensions, these
spacing are ?x, ?y, and ?z (?? = ?? = 0.171 ?? and ?? = 0.200 ??). If ? (??
0
,??) =
? (?
0
,?
0
,?
0
??) is the measured phase space distribution in the volume element located at ?? = ??
0
then the distribution function one unit to the right (in the ?? coordinate) is given by:
?(?
1
,?
0
,?
0
,??) = ?(?
0
+ ??,?
0
,?
0
,??)
The discrete partial derivative used here is the standard finite difference partial derivative:
?? (?)
??
= Lim
?x?0
? (? + ?x) ?? (? ??x)
2 ?x
The spacing of the data considered here is small, so the limit is removed.
?? (?
0
)
??
?
? (?
0
+ ?x) ?? (?
0
??x)
2 ?x
=
? (?
1
) ?? (?
?1
)
2 ?x
This definition of the partial derivative operator is used to for all calculations involving
derivatives in this dissertation. The two forms of this operator used in this text (and the gradient,
as a simple example) are listed below.
The gradient of the scalar quantity ?(??
0
) is calculated as:
163
?
?
? (?
0
,?
0
,?
0
) =
?
?
?
?
?
1
2 ??
(?(?
1
,?
0
,?
0
) ?? (?
?1
,?
0
,?
0
))
1
2 ??
(?(?
0
,?
1
,?
0
) ?? (?
0
,?
?1
,?
0
))
1
2 ??
(?(?
0
,?
0
,?
1
) ?? (?
0
,?
0
,?
?1
))
?
?
?
?
?
B.1
The divergence of the vector field ??(??
0
) = ?
?
(??
0
)?? + ?
?
(??
0
)?? + ?
?
(??
0
)?? is calculated as:
????(??
0
) =
1
2 ??
??
?
(?
1
,?
0
,?
0
) ??
?
(?
?1
,?
0
,?
0
)?
+
1
2 ??
??
?
(?
0
,?
1
,?
0
) ??
?
(?
0
,?
?1
,?
0
)?
+
1
2 ??
(?
?
(?
0
,?
0
,?
1
) ??
?
(?
0
,?
0
,?
?1
))
B.2
The divergence of a symmetric rank two tensor, ?
?
(??
0
), is given by:
164
????
?
(??
0
)? =
???
1
2 ??
??
??
(?
1
,?
0
,?
0
) ??
??
(?
?1
,?
0
,?
0
)?
+
1
2 ??
??
??
(?
0
,?
1
,?
0
) ??
??
(?
0
,?
?1
,?
0
)?
+
1
2 ??
(?
??
(?
0
,?
0
,?
1
) ??
??
(?
0
,?
0
,?
?1
))?
+???
1
2 ??
??
??
(?
1
,?
0
,?
0
) ??
??
(?
?1
,?
0
,?
0
)?
+
1
2 ??
??
??
(?
0
,?
1
,?
0
) ??
??
(?
0
,?
?1
,?
0
)?
+
1
2 ??
??
??
(?
0
,?
0
,?
1
) ??
??
(?
0
,?
0
,?
?1
)??
+???
1
2 ??
??
??
(?
1
,?
0
,?
0
) ??
??
(?
?1
,?
0
,?
0
)?
+
1
2 ??
??
??
(?
0
,?
1
,?
0
) ??
??
(?
0
,?
?1
,?
0
)?
+
1
2 ??
(?
??
(?
0
,?
0
,?
1
) ??
??
(?
0
,?
0
,?
?1
))?
B.3