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 long-standing 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 tri-normal probability distribution function, allows 
ellipsoidally symmetric anisotropy in velocity space; the appearance of such anisotropy has been 
an outstanding issue in weakly-coupled 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 weakly-coupled 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 force-balanced equilibrium.  The dust component of such systems has long 
been suspected to be in a force-balanced equilibrium; the appearance of the tri-normal 
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 stereo-PIV 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:  Quasi-neutrality .......................................................................... 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 multi-normal 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 tri-normal 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 sub-systems ................................................................................... 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:  Tri-normal 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 tri-normal velocity space model. ..................................... 86 
Table 5.1:  Velocity space integral results with the tri-normal 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 bi-normal distribution function.  ................................................ 28  
 
Figure 2.3:  Plots showing the surface of constant probability amplitude for the tri-normal 
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 tri-normal 
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 stereo-PIV 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. 471-2)  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 dust-free regime until the observation 
2 
 
of spoke-like 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:  Quasi-neutrality 
 In a plasma environment the system is assumed to be macroscopically quasi-neutral, this 
condition can be expressed as: 
 
 
 
where the subscript s refers to any of the plasma species.  The symbol e (non-subscript) 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 re-written 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 e-n eq n0 +==
?
s
4 
 
minus sign in the equation for q
d
.  With this definition of the dust charge the quasi-neutrality 
becomes: 
 
ddie
ddei
n Z-nn
)n Z-n-(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 Z-n-n
??
? e
c
?
=
?
=??  
1.2 
where 
c
?  is the electrostatic potential near the dust grain.  Next, the quasi-neutrality 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 electron-ion 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 quasi-neutrality condition, Equation 1.2: 
 
?
2
? =
?1
?
0
??
?
?
?
?
 
 
We now assume no zero-order 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 gas-like and the 
portion of the cloud in the red (lower) box is liquid-like. 
 
?
?
?
?
?
?
?
?
?
=?
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 strongly-coupled dusty plasma, the dust grains form a 
Coulomb solid, the so-called dust crystals, as seen in Figure 1.1 (a).  When 
c
?  ? 1 the dust 
component of the plasma is said to be moderately-coupled, 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 weakly-coupled 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, weakly-coupled, 
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 dust-electron electrostatic repulsion, at the same time 
the dust-ion 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 photo-electrons 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 dust-s 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 re-expressed 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 non-zero ?
?
???
 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 quasi-neutrality condition, Equation 1.1: 
 
i
dd
i
e
n
nZ
n
n
,0,0
,0
1?=  
 
The first case is the so-called "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 so-called "dusty plasma" limit.  
In this limit a non-zero 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, non-spherical 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 weakly-coupled 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.  Weakly-coupled 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 weakly-coupled 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 "non-central" 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 tri-normal 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 multi-normal distribution function is described in Section 2.2.2 and a convenient 
coordinate system rotation that simplifies the multi-normal 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 multi-normal 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 d-dimensional Maxwellian probability distribution function has the form: 
 ?
?
(??;?) =
1
(2??
?
2
)
? 2?
exp ?
?(?? ????)
2
2?
?
2
? 2.7 
where ?? is the d-dimensional velocity vector, ??? is the d-dimensional 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 one-dimensional 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 multi-normal distribution function 
 A multi-dimensional generalization of the Maxwellian probability distribution function is 
the multi-normal probability distribution function.  This distribution is also known as the 
Maxwellian distribution function with tensor variance.  The d-dimensional multi-normal 
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 {?
?
,?
?
} sub-space.  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 "bi-normal" distribution 
function and the three dimensional case will be referred to as the "tri-normal" 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 bi-normal and tri-normal distribution functions is an 
extremely convenient way to understand how the inclusion of tensor variance affects the 
distributions.  We begin with the bi-normal distribution.  Figure 2.2 (a) shows the probability 
amplitude for the bi-normal distribution with the parameters ?
?
?
= ?
?
?
= 2 and ?
?
?
?
?
= 0.5,  
which are combined to give the covariance matrix: 
28 
 
 
Figure 2.2:  Plots showing the bi-normal distribution function with the parameters chosen to be: ?
?
?
= ?
?
?
= 2 and 
?
?
?
?
?
= 0.5.  (a) The bi-normal 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 bi-normal 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 tri-normal 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 tri-normal 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 non-zero 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 non-zero correlation 
and anisotropy.  The full expression for the bi-normal 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 bi-normal 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 bi-normal 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 tri-normal 
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 tri-normal 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 tri-normal 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 tri-normal 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 re-parameterizations 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 tri-normal 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 
multi-dimensional 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 tri-normal distribution is the next most simple distribution function.  
The type of analysis presented in this dissertation for the tri-normal distribution could, in 
principal, be repeated for the tensorized generalization of any one dimensional distribution.  The 
37 
 
tri-normal 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 tri-normal distribution function model 
 The use of the Maxwellian velocity distribution function to describe gaseous system 
dynamics dates back to Clerk-Maxwell?s work in the 1860?s.
25
  In his 1867 paper Clerk-Maxwell 
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 tri-normal distribution function, in the principal axis basis 
set, Equation 2.22, without drift terms (which he discusses in the first and last quoted sentences).  
Clerk-Maxwell?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 tri-normal distributions.  In a follow up paper Clerk-Maxwell
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 force-free equilibrium).  He goes on to show that a system in a 
force-balanced 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 force-free equilibrium, the formulation of statistical mechanics by Clerk-Maxwell, 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 "non-scalar" (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 non-scalar 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 Clerk-Maxwell and 
also very briefly by Spitzer
35
).  A literature search has indicated that the more general non-scalar 
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(X|Y,I) = 
prob(hypothesis|data,I), is known as the posterior probability.  On the right hand side of Equation 
2.28 the term prob(Y|X,I) = prob(data|hypothesis, I) is known as the likelihood function of the 
41 
 
data given the hypothesis and information, prob(X|I) = prob(hypothesis|I) is known as the prior 
probability function for the hypothesis given I, and the term prob(Y|I) = prob(data|I) 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 (non-Bayesian) 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 cut-off 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 cut-off 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 tri-normal 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 tri-normal 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(D|I)) 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 right-most 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 
sub-systems 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 sub-systems that work to 
provide a repeatable and (relatively) consistent plasma environment in which experiments are 
performed.  The physical sub-systems 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 sub-systems 
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 o-rings, 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 six-way 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 5-way 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 o-ring.  The remaining ports of the two crosses are terminated with 
either acrylic windows (for optical access), stainless steel flanges with feed-throughs (to allow 
electrical connections to electrodes or gas system access), or with blank stainless steel caps.  The 
remaining sub-systems 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 super-structure 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 (out-flow, or "sink") and 
neutral gas inflow ("source").  The vacuum pump used for 3DPX is an ULVAC GLD-040 series 
oil-sealed 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 5-10 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 user-selectable 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 KJC-6000). 
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 "sub-system" 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 lightning-like 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 re-enters the dust pile.  The dust grains that form the dust cloud 
and become long-term 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 dust-neutral 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/re-formation 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 (stereo-PIV), and Section 3.2.4 discusses the specific stereo-PIV 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 so-called 
"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).  Stereo-PIV 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 stereo-PIV 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 de-warped 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 
de-warping function is calculated that corrects for the distortion.  De-warped images of the 
calibration plate used in the experiment described below are shown in Figure 3.6 (c).  With the 
Scheimpflug adapter and de-warping 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 "self-calibration."  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 de-warped and rotated into the lab frame.  If the rotated images from camera 1 
and camera 2 differ the de-warping functions are corrected and the images are rotated to the lab 
frame with the corrected de-warping 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 de-warped 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 stereo-PIV 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 de-warping, 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 stereo-PIV measurement, this is done with the so-called 
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 stereo-PIV 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 stereo-PIV 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 "snap-shot" 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 tri-normal models of 
the velocity space and determine that the tri-normal 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 tri-normal  
 
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 sub-images 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 sub-image 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 sub-image 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 tri-normal 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 single-particle 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 single-particle 
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 (multiple-particle) velocity space 
into a single-particle 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 re-scales all three components of the measured velocity 
vector, ??
?,PIV
, into the single-particle velocity space: 
 ??
?
= (??
?,PIV
????) CF (?
?
(??))?  4.1 
where ??
?
 is the vector that has been mapped into the single-particle 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 single-particle velocity space.  It is also noted that the drift velocity calculation is 
performed before the data are re-scaled.   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 sub-space. 
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 sub-space, 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 (zero-drift) 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 tri-normal 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 <Z
i
> 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:  Tri-normal distribution model 
 The calculation of the distribution parameter values for the tri-normal probability 
distribution function model proceeds in much the same way as for the Maxwellian model.  We 
begin with the tri-normal 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 tri-normal 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 tri-normal 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 tri-normal 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 tri-normal 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 tri-normal 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 tri-normal 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 tri-normal model should be selected 
(K>100).  The penalty factor is on the order of 10
-9
, meaning that the likelihood of the tri-normal 
model must be more than 10
11
 larger than for the Maxwellian model to arrive at such a strong 
conclusion (here the tri-normal 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 tri-normal 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 re-analyzed.  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 
 
tri-normal 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 tri-normal 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 tri-normal 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 tri-normal 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 tri-normal 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 tri-normal 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 tri-normal 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 semi-quantitative 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 tri-normal 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 tri-normal 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 tri-normal 
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 tri-normal 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 tri-normal 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 tri-normal 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 tri-normal 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 tri-normal 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, 
semi-quantitatively, 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 tri-normal 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 non-zero 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 non-zero 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 off-diagonal 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 tri-normal.  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, 
semi-quantitatively, 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 force-balanced 
equilibrium.  Section 5.5.2 discusses the difference in the functional form of the thermal heat 
flux vector that one obtains when using the tri-normal 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 re-derive 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 tri-normal 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 re-written, 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 re-inserted 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 re-writing 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 semi-quantitative 
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 velocity-based 
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 tri-normal velocity distribution function 
required a slight generalization of the standard definition for the pressure tensor, ?
?
, given in 
Equation 5.9.  The generalization is fairly straight-forward; 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 tri-normal 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 drift-mediated 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 force-balanced 
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 non-zero 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 
non-zero 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 non-zero and non-uniform 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 tri-normal 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 tri-normal 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 tri-normal velocity space model.  The heat flux vector, in 
terms of the tri-normal 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 tri-normal 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 tri-normal 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 tri-normal 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, force-free, 
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 force-balanced equilibrium.  This characterization of the dust fluid as being in a state of force-
153 
 
balanced equilibrium is self-consistent with the observation of the tri-normally distributed 
velocity space.  The ellipsoidal symmetry of the tri-normal velocity space distribution was shown 
to be the direct result of a system being in a force-balanced equilibrium by Clerk-Maxwell 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 weakly-coupled dusty plasma systems that had been unresolved, open, 
questions.  Broadly, this was done with a theoretical component (the application of the tri-normal 
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 tri-normal 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 tri-normal model, which contained a convenient method of geometric 
visualization and discussed the applicability of the tri-normal 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 tri-normal distribution relates to that developed by Chew, 
155 
 
Goldberger, and Low
34
 (?CGL?), which is a widely applied simplification of the tri-normal 
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 tri-normal 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 weakly-coupled 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 tri-normal 
distribution is "allowed" in the context of phase space based Boltzmann equation-type 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 long-known fact that any three dimensional system 
that is not in a force-free equilibrium must be ellipsoidally symmetric in velocity space, as shown 
by Clerk-Maxwell in 1867
32
. 
156 
 
 The second broad advancement was the development of a method that allows the phase 
space distribution of a weakly-coupled dust cloud to be measured with high spatial resolution.  
Previous analyses of weakly-coupled 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 weakly-coupled 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 semi-quantitative 
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 force-balanced equilibrium.  While many of these properties had been 
suspected in weakly-coupled 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 strongly-coupled "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 weakly-coupled regime is not well understood.  The stereo-PIV 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 
weakly-coupled dusty plasma systems.  The application of the tri-normal 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 force-balanced 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, 469-472 (1879). 
3
 I. Langmuir, Science 60, 392-394 (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, 1-10 (2008). 
12
 B. Liu and J. Goree, Physical Review Letters 100, 1-4 (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, 833-837 (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, 1-8 (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, 641-692 (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, 2374-2378 
(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, 814-817 (2010). 
25
 J. Clerk-Maxwell, Philosophical Transactions of the Royal Society of London. 157, 49-88 
(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, 159-188 (2005). 
31
 A. Azzalini, Biometrika 83, 715-726 (1996). 
32
 J. Clerk-Maxwell, The Quarterly Journal of Pure and Applied Mathematics 32, (1867). 
33
 S. I. Braginskii, Reviews of Plasma Physics 1, 205-311 (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, 2551-2567 (2009). 
40
 J. Westerweel, Measurement Science and Technology 8, 1379-1392 (1997). 
41
 R. J. Adrian, Experiments in Fluids 39, 159-169 (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, 892-896 
(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