UBC Theses and Dissertations

UBC Theses Logo

UBC Theses and Dissertations

Simulation of drifting gas-phase ions of arbitrary shape Chen, Xiu 2003

Your browser doesn't seem to have a PDF viewer, please download the PDF to view this item.

Item Metadata

Download

Media
831-ubc_2003-0204.pdf [ 2.65MB ]
Metadata
JSON: 831-1.0061288.json
JSON-LD: 831-1.0061288-ld.json
RDF/XML (Pretty): 831-1.0061288-rdf.xml
RDF/JSON: 831-1.0061288-rdf.json
Turtle: 831-1.0061288-turtle.txt
N-Triples: 831-1.0061288-rdf-ntriples.txt
Original Record: 831-1.0061288-source.json
Full Text
831-1.0061288-fulltext.txt
Citation
831-1.0061288.ris

Full Text

Simulation of Drifting Gas-Phase Ions of Arbitrary Shape 2. by Xin Chen B.Sc, Nanjing University, China, 1999 A THESIS SUBMITTED IN PARTIAL F U L F I L M E N T OF THE REQUIREMENTS FOR T H E D E G R E E OF M A S T E R OF SCIENCE in The Faculty of Graduate Studies (Department of Chemistry) We accept this thesis as conforming to the required standard  T H E UNIVERSITY OF BRITISH COLUMBIA April 19, 2003 © Xin Chen, 2003  In presenting this thesis in partial fulfilment of the requirements for an advanced degree at the University of British Columbia, I agree that the Library shall make it freely available for reference and study. I further agree that permission for extensive copying of this thesis for scholarly purposes may be granted by the head of my department or by his or her representatives. It is understood that copying or publication of this thesis for financial gain shall not be allowed without my written permission.  Department of Chemistry The University Of British Columbia Vancouver, Canada Date  Abstract  ii  Abstract This thesis is mainly devoted to the description of the first attempt to extend the simulation of drifting ions from the treatment of linear ions to that of arbitraryshaped ones. The simulation is based on molecular dynamics using classical equations of motion. The simulation results for N O drifting in argon are presented. The accuracy of an approximate distribution function is examined in detail. It is found to give results that are generally qualitatively correct, and for many properties, semiquantitatively correct as well. The agreement, however, is not as good as that found for systems with ion-bath gas ratios greatly exceeding unity. The ion-bath gas mass ratio for N O - A r is close to but less than unity, thus increasing the importance of inelastic collisions. The perpendicular velocity coupling term ^ ( ^ l l ) minimum which displays a general behaviour of drifting ions missed in the study of NO -He. Strong velocity-angular momentum coupling is found, and in particular the quadrupolar alignment parameter as a function of the velocity parallel with the field takes on, with decreasing velocities, values that start as negative, become positive, and subsequently decay toward negative. To the best of our knowledge, this is the first report of the decay of this alignment towards negative values at the low end of the velocity distribution. In addition to N O - A r results, simulations of H 2 0 - H e are also presented. This system was used to test modifications of the code to simulate ions of arbitrary shape. An ab initio potential energy surface for H 2 0 - H e was calculated. The mobilities at different field strengths is reported. +  +  n  a  s  a  +  +  +  +  Contents  iii  Contents Abstract  ii  Contents  iii  List of Tables  iv  List of Figures  v  I  1  1  2  Thesis Experimental and Theoretical Backgrounds  2  1.1 1.2 1.3  2 7 8  Simulation Methodology  18  2.1 2.2  18 20 20 23 29 34  2.3 2.4 3  The Drift Tube Technique Applications of Ion Mobility Spectrometry (IMS) in Chemistry . . The Kinetic Theory of Ion Mobility  Frames of Reference Dynamic Simulation 2.2.1 Forces and Torques 2.2.2 Ion Trajectory Evolution Simulation Boxes Approximate Distribution Functions  Results NO+ -  41  3.1  3.2  4  3.1.1 3.1.2 H 0+ 3.2.1 3.2.2 2  Ar  General Description of the Calculation Testing Approximate'Distributions - He PES of the H 0 + - H e Complex Other Results  Conclusions  Bibliography  2  41 41 42 54 55 58 60 63  List of Tables  iv  List of Tables 2.1 Classification of ions according to principal moments of inertia. . 20 3.1 Values of calculated parameters for N O - A r at 300 K for a range of reduced field strengths[37] 3.2 Values of calculated parameters for NO -He at 300 K for a range of reduced field strengths [20] 3.3 Values of calculated parameters for N O - A r at 300 K for a range of reduced field strengths[37] 3.4 Values of fitted parameters for N O - A r at 300 K for a range of reduced field strengths. Most of the coefficients have estimated errors of ± 3 in the last reported digit.[37] +  43  +  44  +  44  +  46  List of Figures  v  List of Figures 1.1 1.2 1.3  Schematic diagram of a drift tube. [3] The arrival-time spectra of K in nitrogen[4, 5] Arrival-time spectra of N ignoring the conversion reaction if the drift tube is so short that the N doesn't have enough time to react with N 2 during the drifting. [6] Arrival-time spectra of N with conversion reaction to form the stable ion, Ng".[6] +  3 5  +  +  1.4 2.1 2.2 2.3  2.4 2.5  6  +  Schematic view of spherical coordinates and definitions of r, 9,4>. Definitions for the Euler angles <j>, ip, and 6 Schematic view of periodic boundary conditions. The parameter rcut is the cutoff radius that is normally applied when calculating the force between two atoms. As can be seen, an atom may interact with one in the neighbouring cell (which is an image of one of the atoms in the simulation cell) because it is within the cutoff radius Drift velocities, as determined by the MD method, as a function of E/N for NO+ drifting in helium at 300 K . [22] (a) Effective total temperatures as a function of drift velocity for N O drifting in helium at 300 K . (b) Partitioning of the rotational and translations! energy between perpendicular and parallel directions. [22]  6 22 25  32 35  +  3.1  3.2  37  Plots of f (v\\) for N O drifting in argon at 300 K with a reduced field strength of 112 Td. The results from molecular dynamics simulations are given by open circles, while those from Eqs. (3.2)(3.4) are given by the dashed, dotted, and solid lines, respectively. 45 Plots of t^i>||)/(u x) for NO+ drifting in argon at 300K for three reduced field strengths. The results from molecular dynamics simulations at 75, 112, and 250 Td are given by the open triangles, circles, and squares, respectively. Also shown is the fitted curve (dashed line) at 112 Td predicted by Eq. (3.5) with the values of /3 listed in Table 3.4. The solid line denotes simulation results for N O drifting in helium at 300 K for a reduced field strength of 100 Td 48 +  v  2  m  +  List of Figures 3.3  3.4  vi  Plots of 2itv±g (vx) for N 0 drifting in argon at 300 K with a reduced field strength of 112 Td plotted at the five values of v\\ indicated. The results from molecular dynamics simulations are given by open circles, while those from Eqs. (2.35), (2.39), and (2.40) are given by the dashed, dotted, and solid lines, respectively. 49 Plots of J]_(vx)/(J ) (lowerpanel) and (v_i_)/(J )(upper panel) as a function of v± for N O drifting in argon at 300 K from molecular dynamics simulations. Results for reduced field strengths of 75, 112, and 250 T d are drawn as open triangles, circles, and squares, respectively. The lines drawn through the symbols are meant to aid the eye only 50 Plots of J]_(v\\)/(J ) (lower panel) and J {v\\)/{J ) (upper panel) as a function of v for N O drifting in argon at 300 K from molecular dynamics simulations. Results for reduced field strengths of 75, 112, and 250 Td are drawn as open triangles, circles, and squares, respectively. The lines drawn through the symbols are meant to aid the eye only 51 Plot of T (v .,v \)/T as a function of v and v\\ for N O drifting in argon at 300 K from molecular dynamics simulations at a reduced field strength of 112 Td 52 Plots of quadrupole alignment parameters as a function of angular momentum, a^\j), (left most panel), parallel velocity, % ( \\)' (middle panel), and perpendicular velocity, ar (v±), (right most panel) for N O drifting in argon at 300 K for three reduced field strengths. The results from molecular dynamics simulations at 75, 112, and 250 Td are given by the open triangles, circles, and squares, respectively. Predictions of Eq. (2.48) using Method II[37] at 75,112, and 250 Td are given by the solid, dotted, and dashed lines, respectively 53 Two-dimensional slice of the PES for the H2 0 -He system as function of r and 6 taken with fixed a = 10° 56 Two-dimensional slice of the PES for the H20 -He system as a function of r and 6 taken with fixed a = 0° 57 Two-dimensional slice of the PES for the H20 -He system as a function of r and a taken with fixed a — 90° 57 Two-dimensional slice of the PES for the H2 0 -He system as a function of r and a taken with fixed 8 = 60° 58 The drift velocities of H 0 drifting in helium as a function of reduced strength both from experiment [7] and MD calculation. 59 +  Vv  2  2  +  3.5  2  2  2  +  3.6 3.7  r  r  1  +  ]  ±  v  0  +  3.8 3.9 3.10 3.11 3.12  +  +  +  +  +  2  Part I Thesis  Chapter 1. Experimental and Theoretical Backgrounds  2  Chapter 1  Experimental and Theoretical Backgrounds In this section, I will introduce some background to the mobility of ions from both experimental and theoretical perspectives.  1.1  The Drift Tube Technique  Drift tubes have been used for many decades to study gas-phase ion collisions and reactivity. The first drift-tube mass spectrometers were developed in the early 1960s at Georgia Tech[l] and Bell Telephone Laboratory[2]. A conventional ion drift tube, see Fig. 1.1, usually consists of an enclosure containing gas, an ion source positioned on the axis of the enclosure, a set of electrodes that establish a uniform axial electrostatic field along which the ions drift, and a current-measuring collector in the gas at the end of the ionic drift path. The drift field E causes the ions of any given molecular composition to swarm through the gas with a drift velocity and diffusion rate determined by the nature of the ions and gas molecules, the field strength, and the gas pressure and temperature. As an analytic tool, we can extract useful information from the output signals. The drift velocity, and longitudinal and traverse diffusion coefficients are the key characteristics of the ion swarm. For the determination of these properties, the ion source is operated in a repetitive, pulsed mode, and the spectrum of  Chapter 1. Experimental and Theoretical Backgrounds Ceramic Insulator Elecfrospray Source  RGA  Ouadrupole Mass Spectrometer  \  G10 Spacer  3  Collision Oynode  Ion Gate  -Hi ii II ) I | M I I I I I I I I I t I I I I I I M I I i I I i I I I M I I I I 1 M I I I i I I j  t  % l I 1 II M I I M I ! ! I I < I I  I U I I I I I I I 1 U  1  M I I I I  i  > < I I I I I I  • I  j  UJ Diffusion Pump ^ Guard Rings  Diffusion Pumps  i  =  i  MicroChannel Plates  Ion Leas  Figure 1.1: Schematic diagram of a drift tube. [3]  arrival times of the ions at the collector is measured electronically. I will now clarify these key physical properties for the ion swarms. After ions arrive at the end plate and pass through the orifice, out of the drift tube, the ions will move into an evacuated region that contains a mass spectrometer and an ion detector. The mass spectrometer can be used to filter out specific ions for study. Mass filters can also be placed before the drift tube to select particular ions. The arrival-time spectrum can be mapped separately for each kind of ion arriving at the end of the drift tube. T h e pressures at which ionic drift-tube experiments have been performed have ranged from about 2.5 x 10~ torr to above 1 atm. O n the assumption that a 2  cross section for collisions of the ions with gas molecules is around 5 0 x l 0 ~ c m , 1 6  the ionic mean free path corresponding to a pressure of 2 . 5 x l 0 cm.  The mean free path at 1 atm is 7 . 4 x l 0 ~  6  - 2  2  torr is 0.23  cm. A t pressures lower than  Chapter 1. Experimental and Theoretical Backgrounds  4  the minimum quoted here the ions would make too few collisions to reach a macroscopic steady state velocity unless an uncommonly long apparatus were used. The ratio of the drift field intensity to the gas number density, E/N, in the drift-tube measurement varies from 0.3 to 5000 Td. One Townsend(Td) is equal to 1 0  - 2 1  V m . The ratio E/N is a very important parameter that determines 2  the average ionic energy acquired from the electric field, that is, the field energy. There are two situations that may happen in the drift tube when making a measurement: 1. The primary ions in the swarm don't react with each other or with the bath gas. In this situation, it is easy to accurately determine ionic drift velocities. Data can be collected for each ion one by one. The pressure in the drift tube should be fairly high to avoid excessive broadening of the arrival-time spectrum due to diffusion, and also to make an accurate measurement of the pressure. As an example, the spectrum of the arrival times for K  +  in nitrogen is shown in Fig. 1.2. The figure in the upper  panel displays seven spectra corresponding to seven different ion source positions. Their arrival times are recorded on a single drawing. 2. The primary ions in the swarm react with the bath gas molecules. A typical example is the experiment on nitrogen ions drifting in nitrogen gas [6]. This is a simple example in which a known primary ion can be converted to a known secondary species that doesn't react further, and retains itself within the drifting time. The conversion reaction is  N+ + 2 N —>N+ + N . 2  2  The effect can be seen by examining the spectrum of the arrival times of N  +  without(Fig. 1.3) and with(Fig. 1.4) the conversion reactions.  Chapter 1. Experimental and Theoretical Backgrounds  MtMVU. TUB OGCftOHMaoa  Figure 2-0-2. (») Arrival-time spectra of "nonreacting" K * ions in nitrogen. Recorded at seven ion source positions (dots) compared with the spectra calculated by the analysis of Section 2-6 (solid curve) (Moscley et al., 1969a). T  i—i—i—i—r—i—i—i—i—i—i—i—i—i—i—i—IT  AB RIVAL TIME (MICROSECONDS  Figure 2-2-1 (b) Arrival-time spectra of "nonreacting" K * ions in nitrogen. Recorded for a drift distance of 4.177 cm (histogram) compared with the spectrum calculated by the analysis of Section 2-6 (smooth curve). Reactions were insignificant at the gas pressure used here (Thomson et al„ 1973).  Figure 1.2: The arrival-time spectra of K+ in nitrogen[4, 5]  5  Chapter 1. Experimental and Theoretical Backgrounds  MM  ' T - ' " ' I T  m  f  1 li \ 1 1 S  /  |  i  ]" V |  I  |'  K* IN NITBOOEN N 121 Td PRESSURE ftOTSTOSR DRIFT DISTANCE 626 CM HISTOGRAM: EXPERIMENTAL DATA SOLID CURVE: ANALYTICAL PROFILE —  1  / J  :  i  \  /  mt m  »  -  > 2  -  m  a  -  1  m  6  /  198  • ,„ m  E —fc-—1 , l 1 j  i  mm  !  mm  1 i no  i  ARRIVAL 10U  1 i T - i. 1 i 1 _L. ten noo itm  1 i  m  m  WICROtf CONM}  Figure 2-2-3. Experimental arrival time spectrum for N ' ion* in nitrogen when conversion to N,, * is negligible compared with the spectrum calculated by the analysis of Section 2-6. In the calculation of the analytical profile a rate coefficient of 2.0 x I0" cm*-s" was used for the reaction converting N " to SM * (2-2-1) (Moseley et. al.. 1969b). Js  1  S  Figure 1.3: Arrival-time spectra of N ignoring the conversion reaction if the drift tube is so short that the N doesn't have enough time to react with N during the drifting. [6] +  +  2  . 1 1 !  itoo  \  i  N,'  HO  i X  j.;  i  « * AND  l  BOB  M  rat-  1  m  i  i  IN ( 1 , tm  TORR  an  -  0HIFT DISTANCE  1  3T.H CM _  1  "™  I  & ~ m 8  «j*  i  PRESSURE  1  i  i  i .  ij  1  \  IK m ... m D  "_!  two  JA two ieot>  i.. i  I  1\  two  i  tm  1  1.  tm  \  .1  M M  L  1 ."  M M  A R R I V A L TIMf M I C R O S E C O N D S )  Figure 2-2-4. Comparison of N and N , * arrival-time spectra showing the effect of the formation of N j * from N . +  +  P i c r u r p 1 4* A r r i v a l - t . i m p s n p r t r a n f N~t~ wit.fi r n n v p r s j i n n r p a r t i n n t o f o r m t b p  Chapter 1. Experimental and Theoretical Backgrounds  7  There are two additional matters worth mentioning here. The first is "end effects " that arise at the ion source or the detector end of the apparatus. These end effects include issues such as the nonuniformity of the electric field at the entrance to the drift tube, and the finite time required for ions to travel from the drift-tube exit aperture to the detector in a drift-tube mass spectrometer. Experimentally, there are techniques that have been developed to treat end effects. The second is related to the detailed shape of the arrival spectrum. By solving a partial differential equation(PDE) governing the drift and diffusion of the ions in the apparatus with a particular geometry, the shape of the arrivaltime spectrum with a close match to the measured arrival-time spectrum can be predicted. In Fig. 1.2, the P D E governing the diffusion of K  +  in nitrogen  gives a beautiful prediction of the shape of the arrival-time spectrum. In practice, drift-tubes can be used to study kinetics, such as the determination of longitudinal diffusion coefficients, and of reaction rate coefficients from arrival-time spectra. In many experiments, drift-tubes are combined with mass analysis in order to increase the ability to identify ions. Without the help of mass analysis, it is impossible to distinguish the behavior of individual species since when ions drift in the tube, they could react with the gas to produce an unexpected assortment of ions.  1.2  Applications of Ion Mobility Spectrometry(IMS) in Chemistry  The term ion mobility spectrometry (IMS) refers to the principles, practice, and instrumentation for characterizing chemical substances through gas phase ion  Chapter 1. Experimental and Theoretical Backgrounds  8  mobilities. In modern analytical IMS methods, ion mobilities are determined from ion velocities that are measured in a drift tube with supporting electronics. Ion mobilities are characteristic of substances and can provide a means for detection and identification. IMS is a useful method for detecting trace concentrations of organic compounds, and has traditionally found many applications in environmental sensing, and in detection of chemical agents. H.W. Ellis, L. A . Vieland, E.A. Mason, and coworkers have collected plenty of data on mobility and diffusion coefficients for different ions drifting in different bath gases at different field strengths.[7, 8 , 9, 10]. Jarrold and coworkers [3, 11] have used IMS to study atomic clusters of silicon, and carbon clusters. Improved apparatus and better electronic-control systems are used to produce purer ion sources and obtain better experimental data. A variation of IMS, FAIMS(High Field Asymmetric Waveform Ion Mobility Spectrometry) [12] is a new technology capable of separating of gas-phase ions at atmospheric pressure (760 torr) and at room temperature. FAIMS can be operated over a wide range of pressures and has been tested above 1500 torr. FAIMS will operate at lower and higher temperatures. FAIMS works by applying strong oscillating electric fields perpendicular to the drifting field. In this way, ions are forced to wiggle at high field strengths as they drift. FAIMS can separate ions from each other[13], focus ions[13] and trap ions in three dimensions [14].  1.3  The Kinetic Theory of Ion Mobility  Drift tube is an open system. The bath gas is in equilibrium. For the collision rate of ion-bath gas is lower than that of gas-gas, the ions is macroscopically in steady state, but microscopically in non-equilibrium. Due to the electric field, the ion system is an anisotropic non-equilibrium system microscopically. The behaviour of ions in drift tubes fits firmly in the framework of non-equilibrium  Chapter 1. Experimental and Theoretical Backgrounds  9  statistical mechanics. In this section, some concepts and basic terms in the transport properties of ions in gases will be clarified and physical descriptions will be expanded upon. One solution to the Boltzmann equation will be explained in some detail. Before beginning the description of drifting ions, certain simplifications will be made. The ideal system is a localized collection of ions of a single type in a gas of uniform temperature and pressure, where the number density, n, of the ions is low enough that their mutual Coulomb forces of repulsion can be ignored. Essentially, when no external field is present, the system is governed by diffusion processes where the net spatial transport of the ions is produced by the gradient of relative concentrations. A diffusive flow exists in the opposite direction of the gradient. The ionic flux density J is given by Fick's law of diffusion, J =  -DVn,  (1.1)  in which D is the scalar diffusion coefficient. If a weak uniform electric field is now applied throughout the gas, a steady flow of the ions along the field lines will form, and be superimposed on the faster random motion from diffusion. The velocity of the center of mass of the ion swarm is called the drift velocity, V D , and is directly proportional to the electric field intensity, E, that is  v  D  = KE,  (1.2)  in which the constant of proportionality K is called the scalar mobility of the ions. In the weak field limit, a beautiful and simple relation, known as the Einstein equation, exists between the mobility and diffusion coefficient, namely,  Chapter 1. Experimental and Theoretical Backgrounds  10  in which k is the Boltzmann constant, T i s the bath gas temperature, and q is the charge of the ion. The electric field can be raised to a level at which the ions acquire an average energy that greatly exceeds the thermal energy of the gas molecules. In the case, the high electric field intensity causes the situation to be more complicated. The thermal energy becomes less important, and two large components of motion are produced by the drift field: a directed component along the field lines and a random component representing energy acquired from the drift field but converted into random form by collisions with gas molecules. So, in this limit, the mobility can vary rapidly as a function of E/N. The kinetic energy distribution of the ions becomes distinctly non-Maxwellian, and ion diffusion becomes inhomogeneous. Therefore, the diffusion coefficient becomes a tensor with the form, I'D  T  D=  0 \0  0  0 ^ (1.4)  0  D  T  0  DJ L  in which DT is the transverse diffusion coefficient that describes diffusion in directions perpendicular to E , and DL is the longitudinal diffusion coefficient characterizing diffusion in the field direction. At higher field strengths, the Einstein relation of Eq. (1.3) does not hold. However, a generalized Einstein relation(GER) can be defined of the form,  q  _  kT dv L  D  kT ^.^dlnK. L  K  (15)  Chapter 1. Experimental and Theoretical Backgrounds  11  in which TT is the transverse (perpendicular) temperature and TL the longitudinal (parallel) temperature. The constant, steady-state drift velocity in Eq. (1.2) is achieved from a balance between the accelerations from the electric field and the decelerations from collisions with the bath gas. We will model the ideal case with the linear Boltzmann equation in the following part. Before doing this, let us discuss the case where the weak electric field condition doesn't hold. The field energy is negligible[15] compared with the thermal energy if (1.6)  in which M and m are the molecular and ionic masses, respectively, qE\ is the energy gained by an ion in moving a mean free path A in the field direction. Using the relation A = 1/(NQ) where yV is the number density of the bath gas and Q is the ion-molecule collision cross section, the inequality in Eq. (1.6) can be written as  This relation defines the weak field limit. It is natural to define the high field limit by the relation in Eq. (1.7) with the inequality reversed. Note that if the number density, electric field strength or temperature changes, a given experiment may change from the low to the high field limit. The kinetic theory of ion mobility and diffusion is complicated if one wants to describe a system at its most detailed level. The goal here is that we want to relate the phenomenological coefficients K and D to the properties of the ions and neutral gas molecules, and any reactions between them. Like so many kinetic problems in physics, everything can be reduced to solving the Boltzmann equation. Solving the Boltzmann equation is a formidable exercise in mathe-  Chapter 1. Experimental and Theoretical Backgrounds  12  matics. Almost no nontrivial exact solutions are known after more than 100 years. First the Boltzmann equation for ions drifting in a gas will be quickly constructed. In the ideal case, diffusion arises from the gradient of the concentration, and the ion swarm should be subject to the continuity equation, dn „ , — + V•J  +  an = 0, (1.8)  d t  J ( r , t) = v n ( r , t) - D • V n ( r , t), D  in which J is the ionic flux density, n the ion number density, a the frequency of the loss of the ions by chemical reactions, D the diffusion tensor and t the time parameter. From the knowledge of kinetics in physical chemistry, we know that the distribution function is a complete quantity to describe our whole system. All the macroscopic information about the system can be calculated once the distribution function is known. The Boltzmann equation is the master equation for the distribution function, f(v,r,t).  By the distribution function, any quantity of  interest can be found by integration. For example, the ion density is  n(r,t) = J f(v,r,t)dv.  In fact, it is the normalization condition for f(v,r,t).  (1.9)  Other quantities of in-  terest are averages of the ion velocity and energy, namely, (v) = n 1 -m(v ) = n 2  1  / v/dv, (  1  L  1  0  )  / (-mv J/dv. 2  Generally, for any function of the ion velocity, ?/>(v,r,t) , the expectation of it  13  Chapter 1. Experimental and Theoretical Backgrounds  can be expressed as ( V ( r , i ) ) = n - yV(v,r,t)/dv.  (1.11)  1  Basically, the Boltzmann equation is just an equation of continuity for / in the six-dimensional space of v and r, except that the flux is in terms of probability density, that is[15]  = E / /[/(v'.r.tJF^Vj.r.tJ-^v.r.^-CV^r.t)] j J J v ' v ' gain term  x v jaj(6,v j)dftjdVj r  r  loss term  = T(/),  (1.12) where a = qE/m,  (1.13)  = \y-V \,  (1.14)  v  rj  j  dn = sin 6d0d4>,  (1.15)  and / ( v , r , t ) and Fj(Vj,r,t) are the distribution functions of the ions and bath gas respectively, and v and V j are the velocities of the ion and bath gas respectively. Unprimed variables represent the initial values before a collision, and those with primes, the final values after collisions. The index j counts over the different species of neutral gas molecules if there are more than one species present. The equation can be rewritten in the following way for a better physical picture, that is  d f= t  ^ J z i streaming term  ~  tL^ii acceleration term  -  £ 1 2 collision term  •  (  L  1  6  )  \  Chapter 1. Experimental and Theoretical Backgrounds  14  The streaming term describes the change due to the free motion of the ions, in other words, it gives a picture of equilibrium for ions without the external fields, some ions move into the vicinity of r and equivalent others move out. The acceleration term describes the part of the change due to an external force altering v. When Boltzmann first derived this equation, he considered a dilute gas in which only two gas molecules at a time are involved in the collisions, that is, binary collisions. This implies that the collision term is a quadratic operator. Collisions involving three particles are very rare, so the binary collision assumption to some extent is good enough for the description of our system provided that the mean free path is large compared to the collision diameter, but is smaller than the diameter of the drift tube. In the collision term, r(/), the loss term describes the loss of ions having velocities in the vicinity of v due to collisions with gas molecules of velocity V . The gain term describes the gain of ions into the velocity region around v' due to collisions of ions with gas molecules of velocity V . The Boltzmann equation written in Eq. (1.12) has some restrictions: 1. There is no term for ion-ion collisions (no term involving a product of two ion distribution functions), since it has been assumed that the density of ions is very low. The treatment of ion-ion collisions by the Boltzmann equation involves fundamental difficulties due to the long-range behaviour of the Coulomb potential. It should be noted that the neglect of ionion collisions means that the motions of different ions in a drift tube are completely independent. This restriction is very important as it allows simulations to run independently, and a lot of ion simulation boxes can be used simultaneously. 2. The number of ions is constant during the drifting time, which means no ion-gas chemical reactions occur.  Chapter 1. Experimental and Theoretical Backgrounds  15  3. The fact that only positions and velocities appears as variable means that no internal degrees of freedom of the ions or gas are considered. Hence, all collisions are elastic. It is possible to release the last two restrictions, but not the first one. It is difficult to find general solutions of the Boltzmann equation. However there is a useful mathematical tool that is worth being introduced  the  Moment Method. The philosophy of the Moment Method[15] is to give up the effort of finding the complete distribution function but to try to obtain the first few moments (or averages) of the distribution. In statistics, there is a very famous formula called the Edgeworth expansion which can used to break down any distribution function into a series of moments of Hermitian polynomials. Principally, the Moment Method applied here is very similar to this expansion. In drift tubes, we are concerned mostly with steady-state phenomena, or at worst with quasi-steady-state phenomena in which the time variation is slow enough to be treated as a sort of perturbation. In short, the idea here is that there may be some chance of success in finding useful solutions of the Boltzmann equation if we are not excessively ambitious. Multiply the Boltzmann equation from the left by any function ip(v) of the ion velocity and integrate over ion velocities to give the moment equation  (1.17)  Integrating by parts and using the inverse collision property of the Boltzmann collision operator,  v jOj(v j)dQjdVjdv r  r  =  v' jOj(v' j)d£l'jdV'jdv' r  r  (1.18)  Chapter 1. Experimental and Theoretical Backgrounds  16  then gives the final equation as,  + V • (n<W>» -  • <V V> = -n^NjiW), j  r  V  (1.19)  in which the linear collision operator Jj is defined by  J.V = iV- 1 j 1  F,#(v)  - V(v')]t;r^idnjdVj.  (1.20)  The linearity of Jj is a result of the fact that the bath gas distribution Fj is an equilibrium Maxwellian distribution, in turn a result of the assumption of the trace ion concentration. Prom the moment equation a variety of behaviours of the ions can be understood. With reasonable physical assumptions, Eq. (1.17) can be reduced to something more tractable. Begin by expanding the distribution function in terms of a complete basis set. By using a perturbation method, the corresponding equation for the term of the corresponding order can be obtained. Viehland[16] applied the principle to the case of diatomic ions. Here, the distribution was expanded in the following form based on the assumption that the dependence of the ion distribution function, / ( r , v , j , i ) , upon position r, velocity, v , angular momentum j and time t is through its dependence upon spatial gradients of the ions number density n(r, t),  /(r,v,j,t)  in which the  /^(vj)  = E/  [  ,  1  W ) - ( - ^ ' n ( r , ( ) ,  (1.21)  are tensors of rank i and the center dot indicates an i-fold  scalar product. After substituting Eq. (1.21) into Eq. (1.12) and expanding it in terms of its spatial gradient, all the terms with the same order of the gradient of n(r, t) are collected to form the equations of different orders. The first two  Chapter  1. Experimental and Theoretical Backgrounds  17  equations obtained in this way are  (a.A+j)/[°)(v,j)=0,  (1.22)  (a-A + j ) /W(v,j) = (v - v )/I°](v,j). D  (1.23)  By solving the equations of different orders, distribution function approximations of a corresponding order can be found. The real solution is substituted by an asymptotic solution which can give a description of the real system qualitatively and quantitatively .  Chapter 2. Simulation Methodology  18  Chapter 2  Simulation Methodology Solving the Boltzmann equation is not trivial, but with the help of computers, the problem can be overcome by simulating the dynamical system. There are several different roads to reach a firm understanding of the behaviour of the drifting ions. Solving the Boltzmann equation directly with numerical methods means spending a lot time determining the collision operator. Monte Carlo [17, 18] and Molecular Dynamics(MD)[19] techniques are other choices for studying ion mobility. In this work, molecular dynamics techniques were used because they provide a means of obtaining important time-dependent information as well as the data necessary to calculate macroscopic properties of the ions. In this chapter, the details of the MD techniques will be described In this research project, previous work[20, 21, 22] on linear molecular ions is extended to the general three dimensional case. In this chapter, two main topics will be discussed, the dynamic description of the drifting ions and the numerical methods applied in the simulation.  2.1  Frames of Reference  In the simulation, three different frames of reference are used: the laboratory frame, the molecular frame, and the body-fixed frame. The laboratory frame is the frame that is fixed in the space of the experiment. Information such as the distribution function, and drift velocity is defined in the laboratory frame. The  19  Chapter 2. Simulation Methodology  molecular frame is fixed on the ions and gives the relative positions of each of the atoms in the ion. The molecular frame can be chosen arbitrarily, and in some cases, may be selected to be the same as the body-fixed frame. The body-fixed frame is the frame that diagonalizes the moment of inertia tensor for the ion. The molecular frame is usually chosen for ease of input of the ion information, and the body-fixed frame is then determined from this information. For the sake of simulating the ion motion, the principal axes of the ion must be found, by solving the secular equation [23],  ( I - / 1 ) - R = 0,  (2.1)  in which the moment of inertia tensor is given by  d xx  ±  I=  lyx ylzx  in which  I  jk  =  'E m (r?6 i  i  jk  - njn^),  is the Kronecker function, and  Ixy  j \ xz 1  lyz Izy  Izz I  j,k = x,y,z,  rf =  r\  x  + r\  y  + r\ ,  6  z  is the radial vector pointing to the i  th  jk  atom  in the ion. The summation i is over all atoms of mass mj making up the ion. This is a typical eigenvalue problem. The vector R represents the eigenvector of the Eq. (2.1). The three eigenvectors from Eq. (2.1) give the directions of the principal axes of the body-fixed frame. The three eigenvalues I , h, and I are a  c  the principal moments of inertia, and are ordered so that I < lb < I a  c  The body-fixed frame depends upon which eigenvector is assigned as the X axis, Y axis and Z axis according to right-hand orthogonality. There is a convention for the assignment of the X,Y and Z axis. Table 2.1 gives a classification for the Ii according to the symmetry of molecules. Equation (2.1) was solved numerically using the lapack mathematical library  20  Chapter 2. Simulation Methodology  Table 2.1: Classification of ions according to principal moments of inertia. Symmetry linear spherical top symmetric top(prolate) symmetric top (oblate) asymmetric top  Principal moments of inertia  Molecular Example  c o , ocs 2  la = h = h h * h  = Ia  /a = /ft # Ic l a t h * Ic  CH CH C1 Benzene Water 4  3  obtained from www.netlib.org. The propagation is performed in the lab and body-fixed frames, so it is necessary to transform between these frames using the transformation matrix R.  2.2  D y n a m i c Simulation  The dynamics simulation is fully based on classical mechanics. No quantum effects are included since classical mechanics is adequate for describing our system at the temperatures of interest. The Correspondence Principle in quantum mechanics states that quantum laws must reduce to classical laws in the limit of large quantum numbers. At high temperature, most states are populated at the large quantum numbers. This part is the kernel of the whole project, and will be broken down into two parts, a description of the forces and torques, and the ion trajectory evolution.  2.2.1  Forces and Torques  In classical mechanics, the forces can be derived from a scalar potential function V if the system is conservative, namely,  Fi = -ViV.  (2.3)  Chapter 2. Simulation Methodology  21  In our system, all the forces, both the interaction between the ions and gas atoms and the interaction with the electric field, are conservative. The potential energy surface (PES) of the ion-gas system will be calculated in order to give the force between the ion and the gas atoms. For most molecules with cyclic symmetry, spherical coordinates are a good choice for describing the PES. The laboratory frame is set in Cartesian coordinates for the sake of convention and convenience. The conversion from spherical coordinates to Cartesian coordinates is given by x = r sin 9 cos <j>,  (2.4)  y = r sin 9 sin </>,  (2-5)  z = rcos0,  (2.6)  while the conversion from Cartesian coordinates to spherical coordinates is given by f = i sin 9 cos <j> 4- ism9sm(f> +  fccosf?,  9 = i cos 9 cos <f> + j cos 0 sin <> / -  fcsinf?,  ci = — i sin <j> + j cos (j>,  (2.7) (2.8) (2.9)  in which the definitions of f, 9, <j> are shown in the Fig. 2.1. In spherical coordinates, the gradient is given by  V = or P  + r 89  r sm 9 d<j>  (2-10)  Care has to be taken when 9 = 0 or n, since the denominator of the component of <j) becomes zero. In these specific cases, the difficulty is overcome by setting the components of 9 and  to be zero.  Chapter 2. Simulation Methodology  22  Figure 2.1: Schematic view of spherical coordinates and definitions of r, 0, <j>. Explicit expressions for the total force acting on the ion in the lab frame is given by F  =  dr  1 dV _sv _isv _ 4> + qEz. r 89 rsin# d(f> f  §  (2.11)  in which the electric field is chosen to lie along the z direction. The torque arises from both the ion gas interaction(—r x V V ) and ion-field interaction(/x x E) and is given in the lab frame by (recalling that f x 6 = <j>, 6 x $ — f and $ x f = 6) 1 6V  &  dV  1  (2.12)  The choice of a good coordinate system to represent the PES partly depends on the symmetry of the ion, if there are no other important coordinate-related requirements. Numerically, the gradient of the PES is calculated from a discrete set of potential values using cubic splines [24]. The Spline Toolbox 3.0 in Matlab 6.5[25] was used because it is very convenient and controllable. The force is discretized on a set of grids on which interpolation can be used. By balancing computation time and proper accuracy these grids can be made fine enough so  Chapter 2. Simulation Methodology  23  that linear interpolation is accurate. In practice, a big array is used to save the discrete force grids for dV/dr, dV/dO and dV/d<f>. When the force at a specific position is needed, the neighbouring grid points are located, and simple linear interpolation is used to calculate the force.  2.2.2  Ion Trajectory Evolution  This section is devoted to explaining how the positions and velocities of the ion evolve with time. A mixed introduction to the mechanics and corresponding numerical methods is given. Chasles's Theorem in classical mechanics states that "the most general displacement of a rigid body is a translation plus a rotation." [23]. For the translational motion, the evolution equation is straightforward, dr =vdt  F  dv =—dt m  (2.13) (2.14)  in which r is the position vector of ion, and its v velocity. In following sections, the superscript "s" represents the space frame/laboratory frame, and the superscript "b" the body-fixed frame. For the rotational motion, the situation is more complicated. The evolution equation in the lab frame is given by dL  s  = t ? , i = x,y,z,  (2.15)  in which L\ is the lab frame angular momentum, and r* the lab frame torque. Equation (2.15) is transformed to the body-fixed frame, since in this frame there is a simple correspondence between angular momentum and angular velocity, mass and moment of inertia, which will make the description of rotations  Chapter 2. Simulation Methodology  24  much cleaner. Otherwise, we have to deal with the moment of inertia matrix. Euler's equations give the description of rotations in the body-fixed frame, namely, + e wL b  ijk  dt  b  =r,  (2.16)  b  k  in which L\ is the angular momentum in the body-fixed frame, r\ the torque in the body-fixed frame,  the angular velocity in the body-fixed frame, and e^j,  the Levi-Civita symbol whose definition is 0, ijk = < +1,  e  if any two labels are the same; if i,j,k are an even permutation of 1,2,3;  — 1,  (2-17)  if ij,k are an odd permutation of 1,2,3.  In the body-fixed frame, it can be shown that  L\ = Iiul  i = x,y,z,  (2.18)  in which Ii is one of the three principal moments of inertia. Combining Eqs. (2.16) and (2.18) then gives I u -ujy (I -I )  = T>,  I u -uj uj (I -I )  =T ,  I 0J -U U (I -I )  = T.  b  x  x  b  y  b  y  v  z  b  z  b  z  z  b  x  z  x  b  x  y  x  y  b  y  (2.19)  b  The evolution of the orientation of the body-fixed frame also needs to be  Chapter 2. Simulation Methodology  25  4  Figure 2.2: Definitions for the Euler angles ci, if), and 6. described. The rotational transformation matrix cos tp cos (j> — cos 6 sin <f> sin tp A  = | cos tp sin <j) + cos 0 cos ci sin ^  — sin ip cos ^ — cos 0 sin ci cos tp  — sin r/> sin <p + cos 0 cos <j> cos ^  sin 6 sin V>  sin 0 cos ip  sin 0 sin <> / — sin 6 sin ci  cos f? (2.20) "  functions as a conversion between the laboratory frame and the body-fixed frame, that is L  s  L  6  =AL ,  (2.21)  =A- L . X  S  The evolution of the Euler angles depends on the angular velocity as follows u\  = ci sin 9 sin ip + 9 cos tp,  0J2 = <p sin 9 cos ip — 6'sin %p,  (2.22)  W3 = <pcos9 + ip,  in which ci, xp and 0 are the Euler angles, defined in Fig. 2.2. These angles are not to be confused with the spherical polar angles in Fig. 2.1. Euler angles have the disadvantage of being susceptible to "Gimbal lock " [26], where attempts to rotate an object fail due to the order in which the rotations are performed. For example, when a rotation about the y axis by n/2 is made,  Chapter 2. Simulation Methodology  26  you will find the problem, namely the rotation matrix becomes singular. This problem is overcome by using quaternions. Quaternions extend the concept of rotation in three dimensions to rotation in four dimensions. Instead of rotating an object through a series of successive rotations (Euler angles), a quaternion allows one to rotate an object through a single arbitrary rotation axis. Quaternions offer another advantage in that they can be interpolated. This allows for smooth and predictable rotational effects. The definition of the quaternions, qi, in terms of Euler angles is,  ^sin(0/2) cos((V> - 0)/2^ 92  sin(0/2) sin((V> - <W/2)  93  cos(6>/2) sin((V> 4-c/>)/2)  w  (2.23)  \a>s(0/2) cos((V> + <t>)/2) J  in which 5Zi=i 9i = 1With this definition, the following relationship between the quaternions and the rotational transformation matrix can be obtained, namely, ^ql + qj-  A=  1/2  9i92 - 9394 ^ 9193 + 9294  9192 + 9394 92 + 9l - I /  9193 - 9294 ^  9293 - 9i94  (2.24)  9293 + 9i94  2  93 + ll ~ V J 2  The evolution equations for the quaternions are 1. (2.25)  Chapter 2. Simulation Methodology  27  in which  M X = l"l u  92  «?],  93  -92  94  9i  -92  -9i  94  -93  92  93  94  (2.26)  y  The numerical propagator is the core of our simulation. Its accuracy will determine the performance of the simulation. This propagator is based on the Verlet algorithm combined with periodic boundary conditions and the cell method. [20, 21] For the translational motion, Eq. (2.14) is discretized using a finite difference scheme[21] to give,  v (t + At) = v (t) + ^At[a (t) + a ( i + At)}, s  3  s  s  (2.27)  in which At is the length of time for which a step is being made. For the rotational motion, discretizing Eq. (2.15) gives,  L {t + At) = L (t) + ^At[r (t) s  s  s  + r (t + At)]. s  (2.28)  A new modified second-order Verlet algorithm is developed in order to make the calculation efficient and computationally simple. The following is the recipe of how to do it. 1. The values of a" are T at time t + At need to be determined in order to S  use Eqs. (2.27) and (2.28). Therefore, the change of the position and the orientation of the body-fixed frame needs to be numerically described to obtain these values. The new position at t + At is determined solely from  28  Chapter 2. Simulation Methodology information at the current time, namely  a(i).  (2.29)  The values of the quaternions in the new orientation at t + At are also determined solely from information at the current time, namely q(i 4- Ai) = q(t) 4- Aiq(i) 4- - A i  2  (2.30)  q(t).  in which q and q are determined from u\ and cjf with the aid of the Eq. (2.25). The quaternions q(i 4- At) have to be normalized. 2. After obtaining values of the new position and orientation at t 4- At, the values of a" and r  s  at time t + At can be calculated in the same man-  ner as was done at time t, by summing all the interactions (forces and torques) from the gas atoms and external electric field using Eqs. (2.11) and (2.12). The new velocity at the next time step can then be calculated using Eq. (2.27). The angular momentum at the next time step can be calculated using Eq. (2.28). 3. Using the updated values of the quaternions q(i 4- A i ) , then allows the updated transformation matrix A to be calculated with Eq. (2.24). This can subsequently be used to transform the lab frame angular momentum and torque into the body-fixed frame, giving L ( i + Ai) and T (t + At). 6  b  Values of u)\(t 4- Ai) and 0J^(t + Ai)can then be obtained using Eqs. (2.18) and (2.19), respectively. With these updated quantities, q(i 4- A i ) and q(i 4- Ai) are then determined from Eq. (2.25). 4. The step is finished up here. The dynamic evolution of the system continues by repeating steps 1, 2 and 3.  29  Chapter 2. Simulation Methodology  Until now, details about the dynamic simulation have been presented. The algorithm above determines trajectories for all particles and is repeated thousands of times during the simulation.  2.3  Simulation Boxes  In the drift tube, ions move fast through the buffer gas. With the external electric field on, the drift tube is an open system in steady state because the ensemble of bath gas atoms is at thermal equilibrium while the ions are not. The ion energy is continuously dissipated through collisions with the bath gas which in turn transfers the heat to a thermal bath. In the M D simulations, the heat is removed using a dissipation mechanism based on the method of Koutselos.[21, 27] In this mechanism, two simulation boxes[21, 28] are introduced. One box simulates the pure gas bath with a fixed temperature, the other box simulates an ion drifting in a supporting gas. The pure gas bath box is simulated using a standard constant temperature equilibrium M D method.  The interaction  potential between the gas atoms is taken to be a Lennard-Jones(12-6) potential, with a proper cutoff distance chosen for the interaction potential. After the gas simulation box relaxes to equilibrium, the ion simulation box is initialized. The ion simulation box is initialized by placing one ion randomly within it, and assigning initial translational and rotational velocities consistent with the equilibrium of the bath gas. Each ion is surrounded by an interaction sphere whose radius is set to the distance at which the ion-bath gas interaction becomes negligible. The bath gas and ion simulation boxes are then superimposed, and the positions of the bath gas atoms relative to the ion interaction sphere are recorded both before and after each time step. If during the time step a bath gas from the bath gas simulation box enters the interaction sphere of the ion, a  Chapter 2. Simulation Methodology  30  copy of its position and velocity is sent to the ion simulation box. The key is that the pure gas simulation box continues according to its dynamics, and the copy gas atoms evolve in the ion simulation box until they leave the sphere of interaction and are discarded. This is the mechanism by which bath gas energy from collisions with the ions is dissipated. Physically, this means that after interaction between the gas atoms and ions, the gas atoms transfer the energy from collisions to other bath gas atoms far enough away that no further interaction with the ions is possible. In other words, the ions always drift through equilibrium gas bath [21]. The ions evolve according to dynamics discussed in the Sec. 2.2. Because the simulation box is finite in size, simulation techniques have to be introduced in order to approximate a macroscopic system, namely the PES must be truncated, and periodic boundary conditions must be used. Cutoff of the P E S The PES should be continuous. As the distance between the bath gas and ion gets smaller, the potential energy approaches infinity. In a numerical simulation, only a finite potential value can be requested so this limits the closeness of the bath gas-ion pair that can be simulated. From statistical mechanics, we know that for the gas to approach the ion very closely, a very large energy will be needed. The bath gas temperature T in the simulation is usually 300K and according to the Boltzmann distribution, the most probable kinetic energy is on the order of kT. Therefore, the smallest bath gas-ion distance in the simulation must have a corresponding potential energy that is much larger(typically more than an order of magnitude) than ffcT. In the simulation code, a counter tracks gas atoms that move much faster than the average and cause the numerical propagator to sample the PES at distances smaller than the inner cutoff. From the simulations of the N O - H e system, we know that the occurrence of these +  Chapter 2. Simulation Methodology  31  events is very low. The simulation also uses an outer cutoff for the bath gas-ion potential chosen to be the distance at which this interactions become negligible. The outer cutoff helps improve the efficiency of the simulation by selecting only these gas atoms whose interaction are important. Periodic Boundary Condition  Periodic boundary conditions [29] are important in molecular dynamics simulations that make a simulation that consists of only a few hundred atoms behave as if it was infinite in size. By using periodic boundary conditions, the simulation box can be reduced to a much smaller size suitable for computer simulation. Figure 2.3 illustrates the concept of periodic boundary conditions. The shaded box represents the system being simulated, while the surrounding boxes are exact copies in every detail - every particle in the simulation box has an exact duplicate in each of the surrounding cells. Even the velocities (indicated by the arrows) are the same. This arrangement is imagined to fill the whole of space. A result of this is that whenever an atom leaves the simulation cell, it is replaced by another with exactly the same velocity, entering from the opposite cell face. So, the number of atoms in the cell is conserved. Furthermore, no atom feels any surface forces, as these are now completely removed. The cutoff radius is always chosen so that an atom can interact with only one image of any given atom. This means that rcut cannot be greater than half the width of the cell.  Parallel Computing Algorithm  In statistical mechanics, ergodic theory[30] states that time and ensemble averages are equal. If one C P U is used to do the simulation experiment, it may take a long time for the code to converge to statistical averages, especially de-  Chapter 2. Simulation Methodology  32  Figure 2.3: Schematic view of periodic boundary conditions. The parameter rcut is the cutoff radius that is normally applied when calculating the force between two atoms. As can be seen, an atom may interact with one in the neighbouring cell (which is an image of one of the atoms in the simulation cell) because it is within the cutoff radius. tailed ones that are functions of velocity and angular momentum. Therefore, a number of the CPUS are used in parallel to perform the simulation. Each C P U simulates a set of ions, thereby increasing the size of the ensemble in the calculation. This larger ensemble then means that statistical averages converge more quickly in time. The code utilizes a mixed mode in averaging, controlling both the ensemble and temporal aspects. A pool of nodes corresponding to the number of sets of ion simulation boxes are prepared in the calculation. One node is very special, called the root node from which different setup copies of the simulations are generated. At the end of the calculation, the root node has the duty of collecting all the computational results from the other nodes and averaging them. The  Chapter  2.  Simulation  Methodology  33  parallelism is realized through the Message Passing Interface(MPI) [31]. A l l the copies evolve independently because of the stochastic property of the system. Brief details of the calculation are given below. Below is the header file, parallel.h used by the code to implement part of the parallel algorithm. C ******************** MPI DECLARATIONS ********************** include  '/usr/local/include/mpif.h'  Integer  max_procs  Parameter ( max_procs = 44) ! maximum number of cpus Integer  root  symbol f o r t h e r o o t node  Parameter ( r o o t = 0)  rank of t h e main p r o c e s s  I n t e g e r myid  rank of t h e p r o c e s s  Integer  r e t u r n code from MPI c a l l 0 ->o..  ierr  I n t e g e r ncpu  number of cpus i n communicator  Integer t a g  t a g f o r t h e message  Integer  rank of d e s t i n a t i o n and source p  dest.src  I n t e g e r stat(MPI_STATUS_SIZE) Integer  buff_Il,  status of receive e t c .  buff_I2  Integer*8 buff_LI1,buff_LI2 COMMON / p a r a l l e l /  myid.ncpu  These declarations are needed to setup parallel communication between processes running on different CPUs. First, three routines, M P I J N I T (initialization), MPI_COMM_SIZE (define the size), and M P I _ C O M M _ R A N K (define the ranks of the nodes) prepare the pool of nodes, and set up the communication among them. The root node is set by the routines. In each node of the pool, the same simulation will be run with different initial conditions by using different random seeds. Each  Chapter 2. Simulation Methodology  34  node has its own bath gas simulation box, and the same number of ion simulation boxes. The routine called M P I _ B A R R I E R is used to make sure all the nodes are at the same calculation stage so that the root node can be ready to collect the data from the rest nodes. There are several different ways for sending and receiving data among the nodes. The routine, M P I J R . E D U C E , is used here. Generally speaking, the routine can sum all the variables in each node and send them to the root node, then the average of the variables can be calculated from them. The M P I calculation greatly increases the convergence rate of the calculation. More specific details about the M P I can be found in Ref. [31]  2.4  Approximate Distribution Functions  This section provides a description of the parameters that characterize drifting ions, as well an introduction to some approximate ion distribution functions. First, consider the effective temperatures T and T for translational and 4  r  rotational motions, respectively. In a steady state, a dilute ion swarm can be characterized by effective translational temperatures in the directions parallel, Xj|, and perpendicular, Tj_, to the field, that is  (2.31)  in which T* is the total translational temperature, and the electric field vector is applied along the laboratory reference frame z-axis. Because the system is not at equilibrium, T", T\_ and Tj' can all have different values, and these can also differ from the bath gas temperature, T. These differing temperatures simply reflect differing partitions of kinetic energy in different directions. Figure 2.4 plots the drift velocities as a function of the applied electric field for N O drifting +  Chapter  2. Simulation  25  50  75  35  Methodology  100  125  150  175 200  Reduced Field Strength, E/N (Td)  Figure 2.4: Drift velocities, as determined by the MD method, as a function of E/N for NO+ drifting in helium at 300 K . [22]  in helium. Generally, drift velocities increase in a linear manner at small E/N, but become non-linear at large E/N. The rotational motion in a steady state can be characterized by the average angular momentum in the laboratory reference frame (J), and average rotational energy. Cylindrical symmetry rules that the average angular momentum should satisfy (J) = 0, since collisions can align but not orient a diatomic molecule in an axially-symmetricflow[32,33]. The effective rotational temperatures in the directions parallel, Tfi, and perpendicular, Tj_, to the field can be defined as  1  -  2kl  +  2kl~  3  T  ±  +  3 " ' 3  ( 2  "  3 2 )  in which I is the moment of inertia of the ion(taken to be linear in this case),  Chapter 2. Simulation Methodology  36  and T is the total rotational temperature. r  The general partition theory in thermodynamics fails here because the ion swarm is an open and non-equilibrium system. According to the physical argument of Viehland[34, 35], energy is transferred into the internal degrees of freedom of the ions only by the collisions with the bath gas atoms, and transferred out of the internal and translational degrees of freedom of the ions only through the translational motion of the bath gas atoms. Hence the characteristic temperatures T ,T must be equal when the system is in a steady state, r  l  that is T* = T .  (2.33)  The values of T and T ' obtained in the simulations are plotted in Fig. 2.5a. r  and the agreement with Eq.(2.33) is remarkable over the entire range of field strengths. Figure 2.5b plots the partitioning of the rotational and translational energy between the parallel and perpendicular directions, in which the important relation T{ = Tf ,  (2.34)  is seen. This indicates that, in addition to Eq.(2.33), energy exchange between the ion and bath in the direction perpendicular to the electric field is on average equally partitioned between translational and rotational modes. This is expected since collisions perpendicular to the field are characterized by the equilibrium distribution of the bath gas, and hence in the steady state lead to this equal partitioning. In the following, some approximate forms for the rotationally averaged velocity distribution functions[21, 22] of the ions in steady state are introduced.  Chapter 2. Simulation Methodology  37  CD Z> CO V—  CD CL  E  CD CD .> O CD  1000 h  LU  Drift Velocity, v (nm/ps) D  Figure 2.5: (a) Effective total temperatures as a function of drift velocity for N O drifting in helium at 300 K. (b) Partitioning of the rotational and translational energy between perpendicular and parallel directions.[22] +  This first form is given by  (2.35)  / > J L , « , , ) = / ™ > , , ) ^ — exp - = £ t -  in which f (v\\) ww  is the Whealton—Woo form[36], 1/2  /  W W  («||)  =  (2.36  2^8 ((v -v ) ) 2  2  l{  D  1/2  x exp •  24/3,52 <  _ y  {v\\ ~v ) D  >  VD  hi  1/2  erfc I  -e '  1 2  2 /3<52 < 4  _ )2 VD  >  («|| - V ) ~ ^ + D  |  Chapter 2. Simulation Methodology  38  i n w h i c h erfc is t h e c o m p l e m e n t a r y error f u n c t i o n , a n d  _ 2S  2  (2.37)  is defined i n terms o f t h e skewness parameter, <5, w i t h  S=)\  11  (2-38)  ({v\\ -  VD) ) 2  1  a n d v±_ a n d v\\ are the velocities p e r p e n d i c u l a r a n d p a r a l l e l t o the field d i r e c t i o n , respectively. T h e v e l o c i t y - c o m p o n e n t c o u p l i n g is i n c o r p o r a t e d t h r o u g h t h e t e r m Uj_(u||).  D e t a i l e d definitions o f these quantities are g i v e n i n R e f . [20]. T h e  fy(v±,v\\)  distribution  was c o n s t r u c t e d t o reproduce t h e drift velocity, a n d  t r a n s l a t i o n a l t e m p e r a t u r e s o f the t r u e d i s t r i b u t i o n . A second a p p r o x i m a t e form[37] is g i v e n b y  tf (v ,v )  =fi (W ,W\\) | l + |a||W||(2W| - 3 )  c  ±  T  {l  ±  - l)(4W|f - 12WJ + 3) + | ( / 3 - 1)  +  X  (2.39)  (Wi - 4Wl + 2) + 2  7 l  (W|  2  -1)W\\  + ^(72-l)(W|-l)(2W|f-l)},  in which  f « - ( i s f i 5 ? H ^ - * Wi =  m  \2kT{  W = n  2kTZ  40)  v< = 2  (2.41)  1/2  m  '  (2  (v\\ -v ) D  v\\ - v = V2{{v\\-v yy/ ' D  2  D  Chapter 2.  39  Simulation Methodology  The other parameters appearing in E q . (2.39) are related to averages over combinations of W± and W\\, namely  «||=  <Wfi>,  (2.42)  0||=  |<W||>.  (2-43)  P±=  \(Wl),  (2.44)  7i =  (WlWtf,  (2.45)  72=  2<W|W?),  (2.46)  in which the average for any function (g) is given by CO  / When f (v±,v\\) v  /-OO  dv\\g(vA.,v\\)f (v±,v\\).  27n;_ ch;_L / V—  -oo  (2.47)  v  L  oo  = /^ (ux,W||), «|| = 7i = 0 and /?y = (3± = <y = 1 • Unlike T  2  the form in E q . (2.35), this form incorporates the velocity-component coupling through fully averaged moments of v± and v\\ via a Gram-Charlier expansion of the distrbution[16, 38].  The distribution f^ (v±,v^) c  is constructed to re-  produce precisely the drift velocity, translation temperatures, as well as all the moments in Eqs. (2.42) —(2.46).  The essence of this approximate form is to  correct the Gaussian distribution f$ (v±,v\\)  by adding higher-order moments.  T  In a study of the N O - H e system[20], the full velocity-angular momentum +  distribution was accurately approximated as  tnf (v±,v^)N(vx,v\\)  F(vj_,v\\,J,y)  v  X  6  X  P  2/fcTfV,^)  [ J  \ \  2IkT? (v ,v )_ df  ±  ll  (2-48)  Chapter 2. Simulation Methodology  40  in which the velocity-dependent temperatures, T" (v±,v\\) and Tff *(v±,v\\), are determined from the velocity-dependent angular momentum averages J\ (v±, v and J]_(v±,v\\) using a procedure described in Ref. [20]. The normalization factor is given by Eq.(28) of Ref. [20]. The coupling of angular moment and velocity are incorporated through the angular tempertures  T^ f,(v±,v^). d  Chapter 3. Results  41  Chapter 3  Results There are two main sets of results to be presented. The first is for the system N O - A r [37], where distribution functions are calculated and the theories [20, +  21, 22] discussed in the last section of Chapter 2 are tested against the results of simulation. The second set of results is for the system H 0 - H e . This is a +  2  new project in which the method for treating linear ions[21, 22, 28] is extended to treat general, rigid, three-dimensional ions. Please refer to Chapter 2 for the methodology of this extension. Results to be presented include calculated drift velocity, alignment parameters, and distribution functions.  3.1  NO+ - A r  3.1.1  General Description of the Calculation  The motivation for studying the N O - A r system was for comparison with pre+  vious results from the NO -He system. One difference between the N O - A r +  +  and N O - H e systems is the mass ratio of the ion to buffer gas. For NO -He, +  +  the ratio is approximately 7.5, which is much larger than 1.0 so that less momentum could be transferred to gas atoms. For N O - A r , the ratio is 0.75 which +  is close to 1.0. The velocity distribution is particularly anisotropic when the mass ratio is close to 1.0. For N O - A r , the kinetics allow more momentum +  tranfer through collisions. Furthermore, this limit provides the critical tests for any approximate distributions. The detailed molecular dynamic simulations of  Chapter 3. Results N0  +  42  drifting through Ar over a range of field strength will be presented.  A number of ab initio calculations [39, 40, 41, 42] have been reported for the ground-state potential energy surface of N O - A r . The PES in Ref. [41] has +  been used in this simulation. The N O ion is treated as a rigid linear rotor +  with a bond distance of 2.0091 ao, and a dipole moment of /z = 0.1360a.u. Simulations are made at reduced field strengths of 75, 112, and 250 Td using a bath gas number density of 0.1 atoms/nm at a constant bath gas temperature 3  of 300K, with a time step of 5 fs. The Ar-Ar bath gas interaction was modelled with a Lennard-Jones function with e = 124K and a = 0.3418nm. With the electric field on, ion-induced dipole effect could exist in PES. But the dipole moment havs an effect on cross setion of collision. N O - A r will present different +  behaviour from the N O - H e due to Ar is more easilier to be induced than He. +  3.1.2  Testing Approximate Distributions  Tables 3.1 and 3.2 list the values of a number of calculated moments and physical quantities[37] at a series of reduced field strengths for N O drifting in Ar and +  He. For a given field strength, the ion drift velocity, VD, and fully averaged (2)  quadrupolar alignment parameter, OQ , are much smaller in the argon bath as compared with the helium bath. Here, the quadrupolar alignment parameter is defined as [43] 4  2)  =2(P (cosG)) . 2  (3.1)  This parameter has a numerical value between +2 (when J lies along z), and -1 (when J lies perpendicular to z). When it is equal to 0, it means that the whole system is isotropic macroscopically. The relations between the effective translational temperatures, given by Eqs. (2.33) and (2.34) are found to be satisfied to within numerical error in the NO -He case[22]. From Table 3.1, +  the relations are not satisfied very well for N O - A r , especially at the highest +  Chapter 3. Results  43  Table 3.1: Values of calculated parameters for N O - A r at 300 K for a range of reduced field strengths[37]. +  Parameter Ti[K]  n  [K] T* [K] ri[K]  r,[ [K] (v'i) [10 m /s ] (vi) [10 m /s ] v [10 m/s] (Jl) [10 amu m /s] (Jp [10 amu m /s] (J ) [10 amu m /s] e  12  2  2  4  4  3  D  e  2  (2)  %  2  2  6  2  2  6  2  2  Reduced Field 0 75 300.0 1088 1602 300.0 300.0 1259 300.0 1323 1272 300.0 300.0 1306 0.166 0.603 0.055 0.925 0.0 0.759 0.280 1.238 0.140 0.595 0.420 1.833 0.0 -0.010  Strength 112 3653 5250 4185 4898 4688 4828 2.02 10.3 1.608 4.584 2.194 6.778 -0.017  (Td) 250 5981 8393 6785 8094 7821 8003 3.31 27.5 2.090 7.575 3.660 11.235 -0.016  field strength. The averages of various translational and rotational temperature are also listed in Table 3.1. It is obvious that the effective temperature of the ions is higher than the temperature of the gas bath. With the field strength increasing, the translational and rotational effective temperatures also increase. For a Gaussian distribution, (i>j_) — 2(v ) , so that deviations from this relation 2  2  L  can tell us whether the behaviour of the system is non-Gaussian. The values of the parameters defined in the Eqs. (2.42) —(2.46), are listed in Table 3.3. All the parameters at 75 Td are biggest among all the field strengths. After 112 Td, the changes of all the parameters are very small. Theories introduced in Sec. 2.4 are tested in the N O - A r case. First, the dis+  tribution as a function of the parallel velocity only, f (v\\) = f v  Q  f (v±,v^)dv± v  are compared. These distributions result after integrating the corresponding fv(v\\,vj_) over the whole range of perpendicular velocities, and for the approx-  Chapter 3. Results  44  Table 3.2: Values of calculated parameters for NO -He at 300 K for a range of reduced field strengths [20]. +  Parameter Ti[K]  n [K]  T [K] T]_ [K] l  T,[ [K] < u > [10 m /s ] <v]_ > [10 m /s ] v [10 m/s] < J'i > [10 amu m /s' ] < J | > [10 amu m /s ] < J > [10 amu m / s ] 2  e  12  2  2  4  4  3  D  e  2  (2)  %  2  2  6  2  2  6  2  2  a  0 299.8 299.7 299.7 295.5 294.8 295.3 0.1661 0.055 0.0 0.2814 0.1404 0.4218 -0.000(1)  10 340.4 350.9 343.9 356.3 338.2 350.3 0.1883 0.071 0.588 0.3395 0.1610 0.5005 -0.017(3)  Reduced Field Strength (Td) 25 50 100 530.0 991.0 2094 767.2 1570 3388 609.2 1184 2525 650.6 1290 2684 543.6 1003 2049 614.6 1204 2472 0.2930 0.5464 1.155 0.176 0.630 2.85 1.413 2.491 4.058 0.6089 1.2072 2.5116 0.2544 0.4696 0.9588 0.8622 1.6742 3.4647 -0.060(3) -0.085(2) -0.090(2)  Table 3.3: Values of calculated parameters for N O - A r at 300 K for a range of reduced field strengths[37]. +  Reduced Field Strength (Td) Parameter 0 75 112 250 a 0.0 0.220 0.193 0.185 ll 1.0 1.202 1.094 1.055 Al 1.0 1.270 1.254 1.255 P± 0.0 0.170 0.149 0.146 7i 1.0 1.225 1.146 1.146 72  imate forms of Eqs. (2.35), (2.39), and (2.40) give  /„>||)  =/^II),  (3.2)  Chapter 3. Results  45  0.4 p  Figure 3.1: Plots of f (v\\) for N O drifting in argon at 300 K with a reduced field strength of 112 Td. The results from molecular dynamics simulations are given by open circles, while those from Eqs. (3.2)-(3.4) are given by the dashed, dotted, and solid lines, respectively. +  v  /„ (f||) = / GC  3 T  J(/3|| - 1)  h ) {1 + |«||W|,(2W|f - 3) + 1  6  8  (3.3)  x(4W -12W, +3)}, 4  2  |  ^  =  ( 2 ^ )  '  6 X P (  -^I " }  <-> 3  4  The accuracy of the three forms is tested in Fig. 3.1. The simulation distribution is asymmetric, having a maximum at a velocity less than vr>, and a high-velocity tail that decays more slowly than its low-velocity one. Obviously, the Gaussian form, fl {v\\), is not representing this shape very well. The T  distributions fy {v\\) and f„{v\\) have better shapes since the high-order moC  ment information is included in them. In particular, the high-velocity tail is  Chapter 3. Results  46  Table 3.4: Values of fitted parameters for N O - A r at 300 K for a range of reduced field strengths. Most of the coefficients have estimated errors of ± 3 in the last reported digit.[37] +  Reduced Field Parameter 0 75 0.0 0.247 Pi 0.0 0.074 P2 0.0 -0.023 Pz 0.0 0.001 PA 0.0 0.015 "oi 0.0 0.106 "02 0.0 -0.035 "03 0.0 0.004 "04 0.0 0.142 "ib 0.0 0.012 "li 0.0 -0.051 "12 0.0 0.016 "13 0.0 -0.001 "14 0.0 -0.032 "20 Jl 0.0 -0.086 "oi "02  a "03 a 11  11  "04 "io  a" "n a "12 a" 11  "13 J\ "l4 "20  0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0  0.068 0.000 -0.001 0.186 -0.009 -0.063 0.024 -0.003 -0.043  Strength (Td) 112 250 0.299 0.411 0.083 0.024 -0.058 -0.118 0.009 0.037 0.062 0.042 -0.021 -0.018 0.002 0.022 0.006 -0.005 0.025 0.003 0.022 0.038 0.010 0.00 -0.011 -0.026 0.00 0.013 0.002 0.010 -0.096 0.055 0.034 -0.021 0.064 0.031 -0.026 -0.020 0.016 -0.011  -0.093 0.041 0.008 -0.006 0.032 0.025 -0.010 0.015 -0.007 0.00  better reproduced by the these two approximate forms. But both of them fail to reproduce the low-velocity tail, and predict maxima that are at too high a velocity. From previous studies[20] , it is known that vx. and v\\ are correlated. One way to examine this correlation is to plot the variation of v\ as a function of Let's take a look at the behaviour of the coupling term  Uj_(«||)  in Eq. (2.35).  Chapter 3. Results  47  In the study of the N O - H e system[20], the parallel-velocity dependence of the +  perpendicular moment was expanded as oo  ^(«n) = <«i>{i + E  -  <*r>».  (-) 35  m=l  in which z\ = (v\\ —VD)/VD-  The corresponding parameters obtained from fitting  the M D results are listed in the Table 3.4. The value ft\ is larger than 0.20 at the reduced field strengths of 75, 112 and 250 Td, which shows that v\ has quite strong dependence on v\\. However, the fit and 0 are also very significant 3  in the expansion. The dependence of ^ ( u y ) at 75, 112 and 250 T d is shown in Fig. 3.2. There is a minimum value at v\\ ns 0. For the N O - H e system[20], only +  one parameter, /3i, was needed to fully account for the temperature variation. However, in the N O - A r case, terms up to fourth order are needed, as seen by +  the values of f} listed in Table 3.4. The curve obtained from the fit at 112 T d m  is shown in Fig. 3.2 where it can be seen to give a good but by no means perfect representation of the simulation curve. Another way to examine the coupling between the velocity components is by plotting distributions at several fixed values of  . These slices are defined  as («x) = 7 0 0 7  J  0  T~Tt  \ •  (-) 3  6  2irvxdvxf (vx,v\\) v  Three different slices, g ^ (vx) from Eqs. (2.35), (2.39), and (2.40) are compared v  in Fig. 3.3 for five different values of z\ spanning from negative to positive. The approximation of Eq. (2.35) is reasonably accurate except when v\\ lies between zero and VD, in which case its accuracy is not so good. In this velocity range, g (vx) has a shape that deviates appreciably from Gaussian. The GramV[{  Charlier approximation of Eq. (2.39) predicts values almost identical to those of Eq. (2.35) except for large negative values of z\, when it has severe oscillations  Chapter 3. Results  -  3  -  2  -  1 (  0 V  I I -  48  1 V  D )  /  V  2  D  Figure 3.2: Plots of v ±(v\\)/(v j_) for NO+ drifting in argon at 300K for three reduced field strengths. The results from molecular dynamics simulations at 75, 112, and 250 Td are given by the open triangles, circles, and squares, respectively. Also shown is the fitted curve (dashed line) at 112 Td predicted by Eq. (3.5) with the values of Pm listed in Table 3.4. The solid line denotes simulation results for N O drifting in helium at 300 K for a reduced field strength of 100 Td. 2  2  +  and predicts negative values for the distribution(these are off the scale of the plot). This oscillatory behaviour is caused by the terms within the braces of Eq. (2.39). The first two terms, a \W\\(2W| - 3) and (3\\ - l)(4Wft - 12Wfi + 3) {  interfere destructively, resulting in a very small value. The oscillatory term (8± - l ) ( W l  _  4 W l + 2)/2 (which is independent of v\\), then dominates the  expression. In addition to velocity correlation, there are also correlations between velocity and angular momentum. To study the dependence of the the rotational temperatures, Jj_(v±,v\\) and  J|(ux>  v\\), on velocity, they are expanded accord-  Chapter 3. Results  0  0.5  1  49  1.5 Vj_/<  2 1/2 Vj_>  Figure 3.3: Plots of 2nvj_g ^ (vj_) for N O drifting in argon at 300 K with a reduced field strength of 112 Td plotted at the five values of v\\ indicated. The results from molecular dynamics simulations are given by open circles, while those from Eqs. (2.35), (2.39), and (2.40) are given by the dashed, dotted, and solid lines, respectively. +  v  50  Chapter 3. Results 0.5 A  0.45  V  0.4  i i i i i i i i i i i i i i i i i i i | i i i i | i i i i | i i i i  " H 0.35 ^  0.3'  i i i i i i i i i i i i i i i i i i i i i i i i i i i i i i i i:  Figure 3.4: Plots of J]_(v±)/(J ) (lower panel) and J (v±)/(J )(upper panel) as a function of v± for N O drifting in argon at 300 K from molecular dynamics simulations. Results for reduced field strengths of 75, 112, and 250 Td are drawn as open triangles, circles, and squares, respectively. The lines drawn through the symbols are meant to aid the eye only 2  2  2  +  ing to [44] ^1IIK^H)=(J1II)  l,  i + 5Xo n=l  i \ \ X, /  1  oo  + E °oiK  - <*r>] (3.7)  m=l 0 0  0  „=lm=l  0  / „,2n  \  \\ J-/  /  V  The corresponding parameters are listed in the Table 3.4. The behaviour of the  Chapter 3. Results  -  3  -  2  -  1  ( v „ - v  51  0 D  ) / v  1  2  3  D  Figure 3.5: Plots of J]_(v\\) / (J ) (lower panel) and J|(v||)/(</ ) (upper panel) as a function of v for N O drifting in argon at 300 K from molecular dynamics simulations. Results for reduced field strengths of 75, 112, and 250 Td are drawn as open triangles, circles, and squares, respectively. The lines drawn through the symbols are meant to aid the eye only. 2  2  +  rotational moments of Eq. (3.7) is shown in Figs. 3.4 and 3.5. As a function of the perpendicular velocity, the moments show simple linear behaviour, as evidenced in Fig. 3.4, except for the lowest field strength where some curvature is seen. This behaviour is the same as that seen with the N O — H e system , +  and implies that terms linear in v\ dominate the expansions of Eq. (3.7). In contrast, as a function of the parallel velocity, the moments have quite complex behaviours that change as a function of field strength. In particular, J (u||) has a minimum close to v\\ — 0, and increases for velocities moving x  away from zero. This behaviour is very similar to that seen in Fig. 3.2. On  Chapter 3. Results  52  argon at 300 K from molecular dynamics simulations at a reduced field strength of 112 Td. the other hand, J (v\\) has a local maximum at approximately 2  = 0, with an  overall double minimum behaviour. This qualitative structure appears at the two higher reduced fields but seems to change to a simple, single-minimum one at the lowest reduced field strength. Clearly, many high-order terms in v\\ will need to be retained in Eq. (3.7) if this behaviour is to be reproduced. This is in sharp contrast to the NO -He system in which only the linear and quadratic +  terms were found to be significant[20]. The overall behaviour in the full-velocity space can be examined in Fig. 3.6 in which values of the velocity-dependent total rotational temperature,  T (v±,v\\) r  are plotted. The local maximum in Fig. 3.5 arises from a similar feature in Fig. 3.6 that appears near v± — 0 for low uy. This feature decays quite quickly as v± increases but is pronounced enough to survive upon integration over v±, producing the bump in Fig. 3.5. Finally, consider the plots in the Fig 3.7 of the alignment parameters as  Chapter 3. Results  53  a) •  E CO CO  a. c 0 E c "CO  0  1  2  3  4  J /<J > 2  2  -3 -2 -1  0 1 2 0 1 2 3 4 5 6 7  ( H " D) V  V  7 V  D  V'/<V > 2  1  Figure 3.7: Plots of quadrupole alignment parameters as a function of angular momentum, a '(J), (left most panel), parallel velocity, a ^(u||), (middle panel), and perpendicular velocity, a^\vj_), (right most panel) for N O drifting in argon at 300 K for three reduced field strengths. The results from molecular dynamics simulations at 75, 112, and 250 Td are given by the open triangles, circles, and squares, respectively. Predictions of Eq. (2.48) using Method II[37] at 75, 112, and 250 Td are given by the solid, dotted, and dashed lines, respectively. 2  0  2  0  +  functions of J , wy, and v±. As a function of angular momentum (left most panel of Fig. 3.7), the approximate distribution, Eq. (2.48), predicts essentially linear behaviour, with higher angular momentum states having stronger preference for negative alignments. This behaviour is precisely that found in previous studies of NO -He.[20] The simulation values at the lowest reduced field strength follow +  this trend but those at the higher reduced field strengths deviate substantially from it for large angular momentum, becoming constant or even increasing slightly in this limit. It is clear that this behaviour will not be reproduced by  Chapter 3. Results  54  Eq. (2.48). As a function of the perpendicular velocity component (right most panel of Fig. 3.7), the simulation values of the alignment parameter show the same behaviour seen previously [20], namely a weak linear but increasing dependence upon velocity. In this case, the approximate distribution predicts essentially constant values, and fails to capture the increase at larger velocities. As a function of the parallel velocity component (middle panel of Fig. 3.7), the simulation values of the alignment parameter have rather large negative values at high velocities. As the velocity decreases, the alignment parameter increases, becomes positive, reaches a maximum positive value, and then decreases towards zero again at the lowest velocities. This final decrease, thereby producing a positive maximum in the alignment parameter, has not been reported previously. Previous studies showed that the alignment parameter does indeed become positive for low parallel velocities but there was no indication that these values would begin to decrease at even lower velocities. The approximate distribution Eq. (2.48) predicts values that are in reasonably good agreement with the simulation values, matching both the position and decay of the positive region at low velocities and the slope of the curves at high velocity. The maximum in the alignment parameter occurs close to z\ = — 1.  3.2  H 0 + - He 2  This project is a first attempt to extend the simulation of drifting ions from the linear to the generally-shaped case. The H 0 - H e system was chosen for two +  2  reasons. First, there is a corresponding set of experimental data for comparison with our simulation data. Second, H 0 2  +  has some symmetry that can be used  to reduce the calculational demand for generating an ab initio PES for the system. The simulation code is modified for generally-shaped ions system and  Chapter 3. Results  55  does not explicitly account for any molecular symmetry.  3.2.1  P E S of the H O + - He Complex a  Because no global PES for H 0 - H e was available in the literature, one was +  2  calculated using Gaussian 98[45]. To calculate a good PES for an open-shell many-atom system is not a trivial job. The UMP2 method and a modified basis set in the same class as the aug-cc-pVTZ [46] was used in the calculation of the PES. At the same time, the Basis Set Superposition Error(BSSE) correction was applied to PES. The high symmetry of H 0 2  +  helps reduce calculation cost.  Spherical coordinates are used to represent the PES. The plane of the H 0  +  2  ion lies in the y-z plane and the origin of the coordinate system is identical with the center of mass of the H 0 2  to the H 0 2  +  +  ion. The relative position of He with respect  ion is described by the spherical coordinate triple (r, 8, a) whose  definition is given in the Fig. 2.1 with a being the same as ip. Just one quarter of the whole coordinate space is calculated, with the rest reproduced by symmetry. The grid points used in the PES calculation were 8 from 0° to 180° in steps of 12°, a from 0° to 90° in steps of 10°, r from l.lOA to 2.50A in steps of 0.05A, 2.60A to 4.00A in steps of O.lA, 4.25A to 5.00A in steps of 0.25A, and 6.00A to 11.OA in steps of 1.0A. Several different two-dimensional slices of the PES are displayed in Figs. 3.83.11. In Fig. 3.8 is plotted the ground state and first excited state potential as a function of r and 8 for a = 10°. The first excited state corresponds roughly to moving one electron from a bonding orbital of H 0 2  +  to a half-filled lone pair  orbital on the oxygen. In this case, the ground state and excited state cross at the 8 = 180°, r « 1.40A and 8 = 0°, r « 1.40A. In practice, the distance of the crossing, 1.40A, is so small that it is almost impossible for gas atoms to reach it because the field is not strong enough to make gas atoms reach the  Chapter 3. Results  56  Figure 3.8: Two-dimensional slice of the PES for the H 0 - H e system as function of r and 9 taken with fixed a = 10°. +  2  distance of crossing point. Or it can be said that the possibility for gas atoms to reach that distance is close to zero. Therefore, this crossing should not affect the simulation results. In Fig. 3.9 is plotted the ground state PES as a function of r and 9 for a = 0°. This geometry corresponds to the case in which He moves perpendicular to the H2 0  +  plane. For fixed 9, the radial curves have wells at small distance. A local  minimum is found around r RS 2.50A and 9 fa 84°. The angle corresponds to the helium atom almost hovering directly above the H 2 0 ion. +  In Fig. 3.10 is plotted the ground state PES as a function of r and 9 for a = 90°. This geometry corresponds to the case in which the helium atom is in the H 0 2  +  plane. The surface in this case is similar to that in Fig. 3.9. The  global minimum of the ground state PES is found in this in-plane geometry at around r RS 3.40A and 9 RS 60°. The well depth is about 15 x I O hartrees. This - 4  Figure 3.10: Two-dimensional slice of the PES for the H20 -He system as a function of r and a taken with fixed a = 90° +  Chapter 3. Results  58  Figure 3.11: Two-dimensional slice of the PES for the H 2 0 - H e system as a function of r and a taken with fixed 6 = 60° +  6 angle is close to the angle between the cyclic axis and the OH bond of  H 0  +  2  so that the global minimum corresponds to a hydrogen-bonding geometry. In Fig. 3.11 is plotted the ground state PES as a function of r and a for 6 — 60°. This slice corresponds to a cone whose point is close to the oxygen atom and whose edges pass through the OH bonds of H 2 0  +  .  For all a, wells  exist in the radial curves. These grow from shallowest at a = 0°(perpendicular to the  H2 0  +  plane) to deepest at the global minimum at a = 90° (in the  H 2 0  +  plane).  3.2.2  Other Results  In Fig. 3.12 is plotted a sequence of drift drift velocities at different field strengths from 6 Td to 100 Td compared to experimental data. The figure shows that the calculated drift velocities agree within the error range of the  Chapter 3. Results  59  Field Strength(Td)  Figure 3.12: The drift velocities of H 0 drifting in helium as a function of reduced strength both from experiment [7] and MD calculation. +  2  experimental data. Thus, both the PES and extended MD method are correct, and are accurately reproducing measured data.  Chapter 4. Conclusions  60  Chapter 4  Conclusions Theoretical research on drifting ions provides a way of compensating for the lack of the experimental tools available for studying the details of ions motion in a bath gas. These theoretical studies can provide us a concrete understanding of the ion bath gas system. In this thesis, a classical molecular dynamics method is utilized. The N O - A r system was expected to provide a stringent test for several +  approximate distributions given by Eqs. (2.35), (2.39) and (2.40) because its ion/bath gas mass ratio is close to one. In general, the approximate form correctly predicted the qualitative behaviour of a range of physical observables at a microscopic level, and for many properties, semi-quantitatively correct as well. The agreement though is not as good as that found for systems with ion-bath gas ratios greatly exceeding unity. For example, in the NO -He system, most +  relations involving velocity-angular momentum coupling were quite simple and linear[20, 22]. The smaller mass ratio of the N O - A r system provided a larger window of +  observation in velocity space, and revealed behaviour not seen previously. First, the velocity correlation term u (i;||) in Eq.(2.35) has a minimun at z\ « — 1 (rex  ferring to Fig. 3.2). This is the first time a turning point for the correlation term has been observed. This property is quite likely generic for many systems. The behaviour of the rotational moments of Eq. (3.7) is shown in Figs. 3.4 and 3.5. As a function of the perpendicular velocity, the moments show simple  Chapter 4. Conclusions  61  linear behaviour, as evidenced in Fig. 3.4, except for the lowest field strength where some curvature is seen. This behaviour is the same as that seen with the N O - H e system, and implies that terms linear in v\ dominate the expansions +  of Eq. (3.7). In contrast, as a function of the parallel velocity, the moments have quite complex behaviours that change as a function of field strength. In particular, J\(v\\) has a minimum close to z\ — 0, and increases for velocities moving away from zero. This behaviour is very similar to that seen in Fig. 3.2. Physically, these minima guarantee that the velocity-angular momentum distribution of Eq. (2.48) will decay quickly to zero, since these moments appear in exponential terms. Furthermore, the dependence of the quadrupole alignment parameter on the parallel velocity showed a maximum in the region of positive alignment occurring at low velocities. The dependence of the quadrupole alignment parameter on the angular momentum is not linear as found in the N O - H e system, which turns up at high field strength. +  The code is designed to simulate generally-shaped ions drifting in bath gases. The H2 0 - H e system was picked for testing, even though it is not so generic +  in the sense of shape. An ab initio PES was calculated at the level of UMP2 employing an aug-cc-pVTZ basis set by using Gaussian 98. There is surface crossing in the PES at small distances. A sequence of calculated drift velocities for reduced field strengths from 0 to 100 E/N were compared with the corresponding experimental ones and were found to agree within their error range. The comparison shows that the simulation code works correctly. The present simulation work gives good insight into the microscopic behaviour of drifting ions system, but studies to date are still limited to small, simple ions. Future theoretical studies could focus on bigger, more complicated ions. To do this, some extension to the simulation needs to be done in the next step of research. For example, new mechanisms for energy transfer between  Chapter 4. Conclusions  62  the ions and gas bath, and coupling to the environment may have to be introduced. The more complicated PES's of large ions, including the crossing of the ground state and excited states at reachable distances for simulation, may require the addition of new methods for treating non-adiabatic effects. Future work also hopes to include the internal degrees of freedom of ions, such as vibration. Currently, ions are treated as rigid bodies, but momentum and energy transfer to their internal degrees of freedom could affect mobilities significantly. When more degrees of freedom are introduced, there could be new numerical approximations required to reduce the computational load.  Bibliography  63  Bibliography [1] W. S. Barnes, D. W. Martin, and E. W. McDaniel. Phys. Rev. Lett., 6:110, 1961. [2] K . B. McAfee and D. Edelson. Proc. 6th Int. Con. Ionizaion Phenomena  Gases, volume 1. SERMA, Paris, 1963. [3] M . F. Jarrold. J. Phys. Chem., 99:11, 1995. [4] J. T. Moseley, I. R. Gatland, D. W. Martin, and E. W. McDaniel. Phys. Rev., 178:234, 1969. [5] G. M . Thomson, J. H. Schummers, D. R. James, E. Graham, I. R. Gatland, M . R. Flaneery, and E. W. McDaniel. J. Chem. Phys., 58:2402, 1973. [6] J.-T. Moseley, M . Snuggs, D. W. Martin, and E. W. McDaniel. Phys. Rev., 178:240, 1969. [7] H. W. Ellis, R. Y . Pad, E. W. Mcdaniel, A. L. Viehland, and E. A. Mason. Atomic data and Nuclear data tables, 17:177, 1976.  [8] H. W. Ellis, E. W. Mcdaniel, D. L. Albrition, L. A. Viehland, S. L. Lin, and E.A. Mason. Atomic data and Nuclear data tables, 22:179, 1978.  [9] H . W. Ellis, M . G. Thackston, E. W. Mcdaniel, and E. A. Mason. Atomic data and Nuclear data tables, 31:113, 1984.  Bibliography  64  [10] L . A . Viehland and E . A . Mason. Atomic data and Nuclear data tables,  60:37, 1995. [11] M . F. Mesleh, J . M . Hunter, A. A. Shvartsburg, G. C. Schatz, and M . F . Jarrold. J. Phys. Chem., 100:16082, 1996. [12] I. A. Buryakov, E. V . Krylov, E. G. Nazarov, and U . K . Rasulev. Int. J. Mass Spectrom. Ion Proc, 128:143, 1993. [13] R. Guevremont and R. W. Purves.  Review of Scientific Instruments,  79:1370, 1999. [14] R. Guevremont, R. W. Purves, D. A. Bamett, and L. Ding. Int. J. Mass Spectrom., 193:45, 1999. [15] E . A. Mason and E . W. McDaniel. Transport Properties of Ions in Gases.  Wiley, NewYork, 1988. [16] L. A. Viehland and A. S. Dickinson. Chem. Phys., 193:255, 1995. [17] H . Skullerud. J. Phys. B, 6:728, 1973. [18] S. N . Lin and J. N . Bardsley. J. Chem. Phys., 66:435, 1977. [19] A. D. Koutselos. J. Chem. Phys., 102:7216, 1995. [20] R. Baranowski and M . Thachuk. Phys. Rev. A, 63:32503, 2001. [21] R. Baranowski and M . Thachuk. J. Chem. Phys., 110:11383, 1999. [22] R. Baranowski and M . Thachuk. J. Chem. Phys., 111:10061, 1999. [23] H. Goldstein and J . Safko C. Poole. Classical Mechanics. Addison Wesley, San Francisco, 2002. [24] C. de Boor and J. E. Marsden. Practical Guide to Splines. Springer-Verlag New York, 1978.  Bibliography  65  [25] Matlab. Spline Toolbox For Use with Matlab, User's Guide, Version 3.0.  http://www.mathworksxom/access/helpdesk/help/toolbox/splines/splines.shtml, 2003. [26] http://guardian.curtin.edu.au/cga/faq/angles.html.  [27] A . D. Koutselos. J. Chem. Phys., 102:7216, 1995. [28] R. Baranowski and M . Thachuk. Physics in Canada, 55:261, 1999. [29] M.P. Allen and D.J. Tildesley. Computer Simulations of Liquids. Oxford  University, Oxford, 1994. [30] A. Katok and B. Hasselblatt. An Introduction to the Modern Theory of  Dynamical Systems. Cambridge University Press, Cambridge, 1996. [31] W. Gropp, E. Lusk, A . Skjellum, and R. Thakur. Portable Parallel Programming with the Message Passing Interface (Scientific and Engineering  Computation). MIT Press, Cambridge, second edition, 1999. [32] R. N . Zare. Angular Momentum. Understanding Spatial Aspects in Chem-  istry and Physics. Wiley, New York, 1988. [33] D. P. Pullman, B. Friedrich, and D. R. Herschbach.  J. Chem. Phys.,  93:3224, 1990. [34] L. A. Viehland. Chem. Phys., 101:1, 1986. [35] L. A. Viehland, S. Lin, and E. Mason. Chem. Phys., 54:341, 1981. [36] J. H. Whealton and S.-B. Woo. Phys. Rev. A, 6:2319, 1970. [37] Xin Chen, R. Araghi, R. Baranowski, and M . Thachuk. J. Chem. Phys.,  116:6605, 2002. [38] L. A . Viehland. Chem. Phys., 179:71, 1994.  Bibliography  66  [39] T. G. Wright, V . Spirko, and P. Hobza. J. Chem. Phys., 100:5403, 1994. [40] T. G. Wright. J. Chem. Phys., 105:7579, 1996. [41] A . M . Bush, T. G. Wright, V . Spirko, and M . Jufek. J. Chem. Phys., 106:4531, 1997. [42] M . Robbe, M . Bencheikh, and J.-P. Flament. Chem. Phys. Lett, 210:170, 1993. [43] A. J . Orr-Ewing and R. N . Zare. Annu. Rev. Phys. Chem., 43:315, 1994. [44] R. Baranowski and M . Thachuk. Phys. Rev. A, 64:62713, 2001. [45] D. C. Rapaport. The Art of Molecular Dynamics Simlation. Cambridge  University, Cambridge, 1995. [46] D. Roth, O. Dopfer, and J . P. Maier. Phys. Chem. Chem. Phys., 3:2400, 2001.  

Cite

Citation Scheme:

        

Citations by CSL (citeproc-js)

Usage Statistics

Share

Embed

Customize your widget with the following options, then copy and paste the code below into the HTML of your page to embed this item in your website.
                        
                            <div id="ubcOpenCollectionsWidgetDisplay">
                            <script id="ubcOpenCollectionsWidget"
                            src="{[{embed.src}]}"
                            data-item="{[{embed.item}]}"
                            data-collection="{[{embed.collection}]}"
                            data-metadata="{[{embed.showMetadata}]}"
                            data-width="{[{embed.width}]}"
                            async >
                            </script>
                            </div>
                        
                    
IIIF logo Our image viewer uses the IIIF 2.0 standard. To load this item in other compatible viewers, use this url:
http://iiif.library.ubc.ca/presentation/dsp.831.1-0061288/manifest

Comment

Related Items