UBC Theses and Dissertations

UBC Theses Logo

UBC Theses and Dissertations

Combining statistics and scattering calculations for improved predictions of molecular collision observables Cui, Jie 2015

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

Item Metadata


24-ubc_2015_september_cui_jie.pdf [ 3.1MB ]
JSON: 24-1.0166524.json
JSON-LD: 24-1.0166524-ld.json
RDF/XML (Pretty): 24-1.0166524-rdf.xml
RDF/JSON: 24-1.0166524-rdf.json
Turtle: 24-1.0166524-turtle.txt
N-Triples: 24-1.0166524-rdf-ntriples.txt
Original Record: 24-1.0166524-source.json
Full Text

Full Text

Combining statistics and scattering calculations forimproved predictions of molecular collision observablesbyJie CuiA THESIS SUBMITTED IN PARTIAL FULFILLMENTOF THE REQUIREMENTS FOR THE DEGREE OFDoctor of PhilosophyinTHE FACULTY OF GRADUATE AND POSTDOCTORALSTUDIES(Chemistry)The University of British Columbia(Vancouver)August 2015c© Jie Cui, 2015AbstractThis thesis describes four types of cold molecular collision systems with increas-ing complexity: from simple atom-diatomic molecule to complicated polyatomicmolecule-polyatomic molecule. The early thesis work is concerned with three openquestions pertaining to collision dynamics of cold molecules. We demonstrate thepossibility of controlling collisional decoherence of ultracold molecules by tuningan external magnetic field. We then provide insight into the feasibility of evapora-tive cooling of molecules, and for the first time incorporate the uncertainty analysisof the potential energy surface (PES) into scattering calculations. In addition, weuse classical trajectory methods to study the effects of the interaction strength andthe geometry of rigid polyatomic molecules on the formation of long-lived colli-sion complexes at low collision energies.The second half of the thesis work is focused on combining statistical method-ology and scattering calculations to address two major problems in molecular dy-namics calculations: increasing computational complexity and uncertainties due toinaccuracies of PES. Using a small number of scattering calculations, we show thatwe can build a Gaussian Process (GP) model to statistically approximate collisionoutcomes for complex molecules, and then perform the uncertainty analysis andthe sensitivity analysis. We also demonstrate that trained by a combination of clas-sical and quantum calculations, a GP model can provide an accurate description ofthe quantum scattering cross sections, even near scattering resonances.iiPrefacePart of the material in Chapter 3 was published in the article: J. Cui and R. V.Krems, Collisional decoherence of superposition states in an ultracold gas near aFeshbach resonance, Phys. Rev. A 86, 022703 (2012). The project was identifiedand designed by Roman Krems and the author. The author developed the code andperformed all the numerical calculations and derived the analytical expressionsrelating the collisional decoherence to scattering lengths.Part of the material in Chapter 4 was published in the paper: J. Cui and R. V.Krems, Elastic and inelastic collisions of 2Σ molecules in a magnetic field, Phys.Rev. A 88, 042705 (2013). The project was designed by R. V. Krems and the au-thor. The original code was contributed by Y. V. Suleimanov. The author modifiedthe code and carried out all the numerical calculations and data analysis.Part of the material in Chapter 5 was published in the paper: J. Cui, Z. Liand R. V. Krems, Collision lifetimes of polyatomic molecules at low temperatures:Benzene-benzene vs benzene-rare gas atom collisions, J. Chem. Phys. 141, 164315(2014). Z. Li developed the classical trajectory methodology for solving cold colli-sions of polyatomic molecules, and wrote the code for classical calculations exceptfor the subroutine for constructing the potential energy surface of benzene dimer,which is the author’s work. The author performed all the numerical calculationsand data analysis.Part of the material in Chapter 7 was accepted for publication in Physical Re-view Letters, arXiv:1503.01432. The project was designed by the author underthe guidance of Roman Krems. Z. Li contributed her codes for the classical andquantum dynamics calculations. The author performed all the design and analysisof the scattering calculations.iiiTable of ContentsAbstract . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . iiPreface . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . iiiTable of Contents . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ivList of Figures . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . viiAcknowledgments . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xiv1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11.1 Cooling techniques . . . . . . . . . . . . . . . . . . . . . . . . . 11.2 Collisions of cold molecules . . . . . . . . . . . . . . . . . . . . 31.3 Challenges for scattering calculations . . . . . . . . . . . . . . . 41.4 Thesis outline . . . . . . . . . . . . . . . . . . . . . . . . . . . . 82 Background material . . . . . . . . . . . . . . . . . . . . . . . . . . 102.1 Quantum mechanical formulation . . . . . . . . . . . . . . . . . 102.2 Differential and integral cross sections . . . . . . . . . . . . . . . 112.3 Single-channel scattering . . . . . . . . . . . . . . . . . . . . . . 132.3.1 Partial wave analysis . . . . . . . . . . . . . . . . . . . . 132.3.2 Computation of integral cross sections . . . . . . . . . . . 152.4 Multi-channel scattering . . . . . . . . . . . . . . . . . . . . . . 172.5 Log-derivative method for scattering calculations . . . . . . . . . 19iv3 Tuning collisional decoherence of ultracold molecules by a magneticfield . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 213.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 213.2 Methodology . . . . . . . . . . . . . . . . . . . . . . . . . . . . 233.2.1 Denisty matrix and decoherence . . . . . . . . . . . . . . 233.2.2 Theory of collisional decoherence . . . . . . . . . . . . . 243.2.3 Complex scattering length from S-matrix . . . . . . . . . 283.2.4 Collisions of Mg(1S) with NH (3Σ) in a magnetic field . . 293.3 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 313.4 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 384 Elastic and inelastic collisions of 2Σ molecules in a magnetic field . . 394.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 394.2 Computation details . . . . . . . . . . . . . . . . . . . . . . . . . 424.3 Spin relaxation mechanisms . . . . . . . . . . . . . . . . . . . . 474.4 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 555 Collision lifetimes of polyatomic molecules at low temperatures: ben-zene - benzene vs benzene - rare gas atom collisions . . . . . . . . . 585.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 585.2 Methodology . . . . . . . . . . . . . . . . . . . . . . . . . . . . 605.2.1 Molecular dynamics simulation of rigid molecule-atom col-lisions . . . . . . . . . . . . . . . . . . . . . . . . . . . . 615.2.2 Initialization of the dynamic state of a collision complex . 655.2.3 Potential energy surfaces for C6H6-Rg and C6H6-C6H6 . . 665.2.4 Details of numerical calculations . . . . . . . . . . . . . . 745.3 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 745.4 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 876 Statistical background of Gaussian process models . . . . . . . . . . 906.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 906.2 Relevant statistical theories . . . . . . . . . . . . . . . . . . . . . 916.2.1 Definition and properties of multivariate normal distribution 916.2.2 Random function (“Process”) . . . . . . . . . . . . . . . 93v6.2.3 Gaussian random function . . . . . . . . . . . . . . . . . 966.2.4 Latin hypercube sampling . . . . . . . . . . . . . . . . . 986.2.5 Maximum likelihood estimation . . . . . . . . . . . . . . 1006.3 Fitting a Gaussian process model . . . . . . . . . . . . . . . . . . 1016.3.1 Why use a Gaussian process model? . . . . . . . . . . . . 1026.3.2 Prediction based on Gaussian process model . . . . . . . 1056.3.3 Maximum likelihood estimation of model parameters . . . 1076.4 Two-step Gaussian process models . . . . . . . . . . . . . . . . . 1086.5 Bayesian calibration of numerical models . . . . . . . . . . . . . 1146.5.1 Introduction to Bayesian inference . . . . . . . . . . . . . 1146.5.2 Metropolis-Hastings MCMC algorithm . . . . . . . . . . 1176.5.3 Bayesian modelling with Gaussian processes . . . . . . . 1176.5.4 Bayesian model calibration . . . . . . . . . . . . . . . . . 1207 Gaussian Process Model for Collision Dynamics of Complex Molecules1277.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1287.2 GP model of scattering observables . . . . . . . . . . . . . . . . . 1297.3 Gaussian process for interpolation of scattering data . . . . . . . . 1327.3.1 Mean lifetime dependence analysis for the C6H6 – Ar system1327.3.2 Anaysis of Interaction between Mass and Polarizability . . 1377.3.3 Fitting a High-Dimensional Model . . . . . . . . . . . . . 1397.4 Bayesian calibration of dynamical calculations . . . . . . . . . . 1437.4.1 Fitting experimental data . . . . . . . . . . . . . . . . . . 1437.4.2 Fitting quantum calculations . . . . . . . . . . . . . . . . 1467.5 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1468 Outlook . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 149Bibliography . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 153viList of FiguresFigure 3.1 Energy levels of the NH molecule in a magnetic field. . . . . . 33Figure 3.2 The real (upper panel) and imaginary (lower panel) parts ofthe Mg – NH s-wave scattering lengths: solid curves – state a;dashed curves – state b; dotted curves – state d of Figure 3.1. 34Figure 3.3 Upper left panel: The absolute difference of the scatteringlengths for NH molecules in states a and b (labeled in Fig-ure 3.1) as a function of the magnetic field. Lower left panel:The t = 0 decoherence rate of NH molecules prepared in thesuperposition of states a and b as a function of the magneticfield. Right panels: The time dependence of the density matrixelements ρela,a (dashed lines), ρelb,b (solid lines) and |ρela,b| (dottedlines) at two magnitudes of the magnetic field. . . . . . . . . . 35Figure 3.4 The quantity Q = αa±√Rad/CT 1/2− (βa−βd)2: solid curve– minus sign; dashed curve – plus sign. . . . . . . . . . . . . 37Figure 4.1 Cross sections for elastic (upper panel) and inelastic (lowerpanel) collisions of 2Σ molecules with γSR = 0.0415cm−1 andBe = 4.2766cm−1 in a magnetic field of 100 Gauss. The col-ored bars correspond to the different number of rotational statesincluded for each molecule in the basis set expansion (4.4),ranging (left to right) from 3 to 7. . . . . . . . . . . . . . . . 46viiFigure 4.2 Cross sections for elastic (left panel) and inelastic (right panel)collisions of 2Σ molecules with γSR = 0.0415cm−1 and Be =4.2766cm−1 computed for the collision energy 10−3 cm−1 andthe magnetic field 100 Gauss. The colored bars correspondto the different number of rotational states included for eachmolecule in the basis set expansion (4.4), ranging (left to right)from 3 to 7. The results are averaged over 35 calculations withdifferent interaction potentials. . . . . . . . . . . . . . . . . . 48Figure 4.3 Cross sections for elastic (upper panels) and inelastic (lowerpanels) collisions of molecules with γSR = 0.0415cm−1 andBe = 4.2766cm−1 computed for the collision energy 10−3 cm−1and the magnetic field 100 Gauss. The results in the left pan-els are obtained with a single potential energy surfaces corre-sponding to λ = 1. The results in the right panels are averagedover the cross sections computed with 50 different potentials.The calculations are performed with N ≤ 2 for molecule A inthe basis set expansion (4.4). The colored bars in each panelcorrespond to the different number of rotational states includedfor molecule B in the basis set expansion (4.4), ranging (left toright) from 3 to 11. . . . . . . . . . . . . . . . . . . . . . . . 49Figure 4.4 Cross sections for elastic (upper panel) and inelastic (middlepanel) collisions of molecules with Be = 4.2766 cm−1 com-puted as functions of γSR for the collision energy 10−3 cm−1and the magnetic field 100 Gauss. The lower panel shows theelastic-to-inelastic cross section ratio. The results are averagedover 48 calculations with different interaction potentials. Thevertical bars indicate the 2σ interval of the cross section valuesand their ratios. The insets illustrate the polynomial regressionfits. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 51viiiFigure 4.5 Cross sections for elastic (upper panel) and inelastic (middlepanel) collisions of molecules with γSR = 0.0415 cm−1 com-puted as a function of Be for the collision energy 10−3 cm−1and the magnetic field 100 Gauss. The lower panel shows theelastic-to-inelastic cross section ratio. The results are averagedover 20 calculations with different interaction potentials. Thevertical bars indicate the 2σ interval of the cross section valuesand their ratios. . . . . . . . . . . . . . . . . . . . . . . . . 53Figure 4.6 Cross sections for elastic (upper panel) and inelastic (middlepanel) collisions of molecules with γSR = 0.5 cm−1 computedas a function of Be for the collision energy 10−3 cm−1 and themagnetic field 100 Gauss. The lower panel shows the elastic-to-inelastic cross section ratio. The results are averaged over35 calculations with different interaction potentials. The verti-cal bars indicate the 2σ interval of the cross section values andtheir ratios. The insets illustrate the polynomial regression fits. 54Figure 4.7 The lowest magnitude of the expected elastic-to-inelastic crosssection ratio for selected 2Σ molecules. The results are ob-tained by interpolation of the data in Figures 4.4, 4.5 and 4.6using the spectroscopic constants borrowed from Ref. [100]. . 57Figure 5.1 (a) The internal coordinates used to describe the Rg – C6H6collision complex: R is the relative distance between centers ofmass, θ (the polar angle) is the angle between R and the 6-foldsymmetry axis of benzene, and φ (the azimuthal angle) is theangle between the projection of R onto the molecular plane andthe x axis, which is perpendicular to a C-C bond of benzene.(b) The relative geometries of the Rg–C6H6 collision complexfor the out-of-plane, vertex-in-plane, and side-in-plane config-urations. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 68ixFigure 5.2 Radial dependence of the interaction energy of benzene withHe (black solid line), Ne (red dashed line), Ar (blue dottedline), Kr (green dot-dashed line), and Xe (purple long-dashedline) compared with the corresponding benzene-benzene inter-action averaged over the normal angles of the molecules (ma-genta dot-dashed line). Upper panel: out-of-plane configura-tion; central panel: vertex-in-plane configuration; lower panel:side-in-plane configuration. The geometries are illustrated inFigure 5.1. . . . . . . . . . . . . . . . . . . . . . . . . . . . 69Figure 5.3 (a) Benzene as a combination of two rings: inner carbon ringand outer hydrogen ring. (b) Model for the interaction betweenthe two rings: R is the relative distance between the centers ofmass, θ (the polar angle) is the angle between R and the 6-fold symmetry axis of benzene, φ (the azimuthal angle) is theangle between the projection of R onto the molecular planeand the x axis, which is perpendicular to a certain C-C bondof benzene, and γ (the normal angle) is the angle between the6-fold symmetry axes of the two molecules. (c) The relativegeometries of the benzene dimer for the sandwich, T-shapedand parallel-displaced (PD) configurations. . . . . . . . . . . 71Figure 5.4 The mean lifetimes of the quasi-bound collision complexes ofbenzene with He (black circles), Ne (red squares), Ar (blue di-amonds), Kr (green up-triangles), Xe (purple down-triangles)as functions of the collision energy. The inset shows the distri-butions of short trajectories and long trajectories for each raregas at a collision energy of 2 cm−1. For each collision energy,the results are averaged over the rotational energy of benzenewith the Boltzmann distribution corresponding to the tempera-ture 6 K. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 76Figure 5.5 The mean lifetimes of the quasi-bound collision complexes ofrare gas atoms with benzene as functions of temperature. . . 77xFigure 5.6 Left panel: The mean lifetimes of the quasi-bound collisioncomplexes of benzene with He (black circles), an artificial atomdescribed in the text (green crossed circles), benzene (bluestars), and Xe (purple down-triangles) as functions of colli-sion energy. Right panel: The mean lifetimes of the quasi-bound collision complexes of benzene with He (black circles),benzene (blue stars), Ne (red squares), and Xe (purple down-triangles) as functions of temperature. . . . . . . . . . . . . . 80Figure 5.7 The mean lifetimes of the quasi-bound collision complexes ofbenzene with atoms of varying mass (indicated on the hori-zontal axis) computed with the potential surfaces for the in-teraction of benzene with He (black circles), Ne (red squares),Ar (blue diamonds), Kr (green up-triangles), Xe (purple down-triangles) at a collision energy of 5 cm−1. . . . . . . . . . . . 81Figure 5.8 The time dependence of the relative distance R, the transla-tional energy (black solid line), the rotational energy of themolecule (red dashed line), and the interaction potential (bluedotted line) for a collision of Xe with benzene at a collisionenergy of 2 cm−1. . . . . . . . . . . . . . . . . . . . . . . . 82Figure 5.9 The time dependence of the relative distance R, the polar angleand the azimuthal angle as defined in Figure 5.1 for a collisionof Xe with benzene at a collision energy of 2 cm−1. . . . . . . 83Figure 5.10 The time dependence of the distance R, the translational energy(black solid line) of the relative motion, the rotational energyof the first molecule (red dashed line), the rotational energyof the second molecule (green dash-dotted line), and the in-teraction potential (blue dotted line) for a benzene - benzenecollision with a collision energy of 2 cm−1. . . . . . . . . . . 85Figure 5.11 The time dependence of the distance R, polar angle θ , az-imuthal angle φ and normal angle γ as defined in Figure 3 fora benzene - benzene collision with a collision energy of 2 cm−1. 86xiFigure 5.12 The lifetime of 400 randomly sampled trajectories for benzene- benzene collisions as a function of the difference between themaximum and minimum normal angles γmax − γmin sampledduring the collision. . . . . . . . . . . . . . . . . . . . . . . 88Figure 6.1 Scatter plot of 1000 samples from a bivariate normal distribu-tion. Each circle represents a sampled vector y with elementsy1 and y2 plotted on the X-axis and Y-axis, respectively. . . . 94Figure 6.2 A random function F(·) mapping from a set Ω of outcomes ωin a sample space to deterministic functions f (·). . . . . . . . 95Figure 6.3 Five draws from a random function . . . . . . . . . . . . . . 96Figure 6.4 (a) a 4× 4 Latin square (b) a completed 4×4 Sudoku grid . . . 99Figure 6.5 left panel: 20 functions fi(x) randomly drawn from a ran-dom function F(x); right panel: 20 functions fi(x|yn) passingthrough all sampled data yn (red circle) randomly drawn fromthe same random function. . . . . . . . . . . . . . . . . . . . 103Figure 7.1 (a) and (b): The lifetime dependence on the rotational temper-ature and collision energy for C6H6 – Ar collisions. . . . . . . 133Figure 7.2 The predicted mean lifetime based on a 20-point design anda 50-point design versus true mean lifetime of the equispacedgrid of 100 points when power parameter is set to be 1.95 . . . 134Figure 7.3 The predicted mean lifetime based on a 20-point design anda 50-point design versus true mean lifetime of the equispacedgrid of 100 points when power parameter is set to be 2. . . . . 135Figure 7.4 The surface produced by the GP model. The lines connect thevalues (circles) computed from the classical trajectories withthe values predicted by the GP model. . . . . . . . . . . . . . 137Figure 7.5 The 3D global response surface produced by the GP model forC6H6 – Rg collision lifetimes with the atomic mass and thePES depth plotted on the x- and y-axis at Tr = 4 K and E = 4cm−1. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 139xiiFigure 7.6 The contour plot produced by the GP model for C6H6 – Rgcollision lifetimes with the atomic mass and the PES depthplotted on the x- and y-axis at Tr = 4 K and E = 4 cm−1. . . 140Figure 7.7 Accuracy of the GP model with variable PES parameters forthe prediction of the collision lifetimes. The scatter plot com-pares the predicted values with the computed values. The errorof the GP model is the deviation of the points from the diagonalline. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 141Figure 7.8 Relative effect of the variation of Tr, E and the PES parameterson the collision lifetimes. The blue area of the bars showsthe uncorrelated contribution of the corresponding variable andthe green area – the effect that depends on one or more othervariables. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 142Figure 7.9 Energy dependence of the collision lifetime for Ar - C6H6 withthe error interval obtained by varying all the individual PESparameters by ±3 %. . . . . . . . . . . . . . . . . . . . . . . 143Figure 7.10 Randomly generated pseudo-observation data (red circles) vsthe results obtained from the CT calculations. The off-set ofthe data points from the dashed line represents the inherentinaccuracy of the scattering calculations. Inset: Accuracy ofthe GP model. . . . . . . . . . . . . . . . . . . . . . . . . . . 145Figure 7.11 The total 15000 MCMC simulations for θ2(left panel), last10000 simulations (out of 15000) of all eight potential param-eters. The burn-in period is 5000. . . . . . . . . . . . . . . . 146Figure 7.12 GP model (solid curve) of quantum scattering cross sections(squares and circles) for C6H6 - He collisions with CT calcu-lations (not shown) serving as an intermediate model. The 60squares are the quantum results used to train the models. Theinset shows the GP model (solid curve) trained by the quantumresults without CT calculations. . . . . . . . . . . . . . . . . 147xiiiAcknowledgmentsForemost, I would like to express my sincere gratitude to my supervisor Prof. Ro-man Krems for his continuous support of my Ph.D study and research. His pa-tience, motivation, enthusiasm and broad knowledge, make him a great mentor.I would like thank Prof. Gren Patey and Prof. Alex Brown for proof-reading mythesis. I also would like to thank the rest of my thesis committee for their insightfulcomments and inspiring questions to my thesis.A special thanks goes to Zhying Li for her generous help and fruitful collab-oration throughout my PhD study. Her advice on both research as well as on myfuture career is very valuable. I also want to thank all my previous and currentgroup members for the great time I spent in our group. I really enjoy the atmo-sphere in our office and appreciate their friendship and support.Last but not least, I must acknowledge my beloved husband Ping. He hasbrought a lot of happiness and fun to my life. Without his love, encouragement and24-hour computer technical support, I would not have finished this thesis.xivTo my husband Ping and mygrandparentsxvChapter 1IntroductionSince the first demonstration of laser cooling and trapping of neutral atoms in the1980s [1], the field of ultracold atoms has achieved substantial progress, leadingto the creation of Bose-Einstein condensates (BECs) [2, 3] and Fermi degener-ate gases of ultracold atoms [4], the observation of the Bardeen-Cooper-Schrieffer(BCS) superfluid transition [5], as well as many other interesting discoveries [6, 7].Inspired by the achievements with ultracold atoms, effort has been made to cre-ate ultracold ensembles of molecules. The rich internal structure of moleculespromises new physics absent in cold atoms and a wide range of applications, suchas precision measurement [8], quantum computing [9] and controlled chemistry[10]. Despite the enormous success achieved in experiments with atoms, laser cool-ing does not work for the majority of molecules due to the complexity of molecularstructure [11]. In order to create ultracold molecules, researchers have developedseveral alternative methods. In this chapter, we give a brief review of current exper-imental and theoretical work on the creation of dense ensembles of ultracold polarmolecules, and summarize several questions that remain unanswered. The chapterconcludes with a summary of motivation for the work presented in this thesis.1.1 Cooling techniquesCurrent techniques for producing cold and ultracold molecules fall into two cat-egories: indirect methods and direct methods. There are many excellent reviews1of the experimental developments [12, 13]. This section only offers a snapshot ofsome main techniques related to the work in this thesis.The indirect cooling methods start with ultracold atomic gases with high phase-space density and then convert atomic pairs into molecules [14]. Feshbach reso-nance linking [15–17] and photoassociation [18–20] are two of the most widelyused approaches for this purpose. Indirect methods have so far achieved temper-atures as low as 300 nK [21]. The biggest drawback of the indirect techniques isthat they are limited to molecular species composed of atoms amenable to lasercooling.Direct cooling methods start with hot molecules and cool them to low temper-atures (T ∼ 10 mK to 1 K). Examples of such methods include buffer-gas cooling[22], sympathetic cooling [23] and evaporative cooling [24]. The buffer-gas cool-ing technique relies on the thermalization of molecules in a cold buffer gas [22].The thermalization process is driven by momentum transfer in elastic collisions ofmolecules with the buffer-gas atoms, so it appears to be general and applicable toany molecular system. It has been successfully used for cooling CaH [22], CaF[25], NH [26, 27], PbO [28], CrH and MnH [29] to T ∼ 1 K. Recently, buffer-gascooling was extended to complex polyatomic molecules naphthalene (C10H8) [30]and benzonitrile (C7H5N) [31].Buffer gas cooling methods can only cool molecules to the cold regime T ∼1–10 mK, where molecules can be confined in an external field trap. To reach theultracold regime, further cooling can potentially be achieved by sympathetic cool-ing or evaporative cooling. Both techniques require that elastic collisions, whichthermalize molecules, far exceed inelastic or reactive collisions, which lead to traploss. Sympathetic cooling was initially developed for cooling singly-charged ions[23, 32] and molecular ions, such as MgH+[33], NH+4 and C2H+5 [34]. Recently,there were several experimental observations [35] and proposals for sympatheticcooling of neutral polar molecules [36–39]. Evaporative cooling is similar to theprocess of cooling a cup of hot coffee: the hottest water molecules escape fromthe cup as steam, leading to a colder temperature of molecules left behind in thecup. Evaporative cooling is commonly done with alkali metal atoms in magnetictraps [24] for the creation of BECs. By now, the only experimental demonstrationof evaporative cooling of molecules was reported for an ensemble of OH(2Π) rad-2icals. In this experiment, the temperature of the molecular ensemble was reducedfrom 51 to 5.1 milliKelvin [40]. It was shown that inelastic OH–OH collisions canbe suppressed in this temperature interval due to the peculiar structure of moleculesin the 2Π electronic state.1.2 Collisions of cold moleculesIn this section, we discuss four types of collision systems of increasing com-plexity: atom–diatomic molecule, diatomic molecule–diatomic molecule, atom–polyatomic molecule and polyatomic molecule–polyatomic molecule. These fourtypes of collision systems also reflect the progress of the research projects reportedin this thesis.The experiments on cooling molecules and the promise of new applicationshave stimulated intense research of the effects of external fields on ultracold molec-ular collisions [41], including the possibility of using external fields for tuning elas-tic and inelastic cross sections and controlling chemical interactions of molecules.The rich internal structure also makes ultracold molecules a promising platformfor the implementation of a quantum computer. In quantum computing, a qubit,the quantum analogue of the classical bit, is a coherent superposition of quantumstates. Qubits can be stored in the long-lived internal states of molecules. More-over, the long-range dipole-dipole interaction between molecules makes it possibleto engineer fast logic gates between remote qubits. However, the main practicalobstacle is decoherence. Molecular collisions are an important source of deco-herence in experiments with atoms and molecules in the gas phase. This raisesthe first question considered in this thesis: is it possible to suppress collisionaldecoherence by applying external fields?Motivated by experimental demonstration of buffer-gas cooling of a great vari-ety of diatomic molecules in a magnetic trap, current studies aim at further coolingmolecular ensembles to ultracold temperatures. The evaporative cooling, as a po-tential way to bridging the milliKelvin – microKelvin temperature gap, relies onelastic molecule-molecule collisions. If molecules are trapped in excited states,the evaporative cooling process is hindered by inelastic collisions which releaseenergy and expunge the colliding species from the trap. Early theoretical studies3of atom - molecule collisions in a magnetic field found that inelastic scatteringis driven by an interplay of intra-molecular spin-dependent interactions and theanisotropy of the atom–molecule interaction potentials [42–44]. The interactionanisotropy in molecule - molecule collisions is generally much larger than that inatom - molecule collisions. This leads to the second question addressed in thisthesis: can molecules be evaporatively cooled? what are the criteria for moleculesto be good candidates for evaporative cooling experiments?In addition to getting temperatures lower, another trend of current researchin molecular physics is getting ultracold molecules larger, i.e., extending cool-ing techniques to complex polyatomic molecules [45, 46]. As a possible coolingmethod for this purpose, buffer-gas cooling may be impeded by various unwantedprocesses, such as resonant scattering and the ensuing formation of atom - moleculeclusters. The probability of cluster formation is clearly related to the complexity ofatom-polyatomic molecule system. Li et al. [47, 48] performed classical trajectorycalculations of collision lifetimes for polyatomic molecule complexes in a buffer-gas cooling experiments at a temperature of T = 6.5 K to elucidate the effects of themolecular size and geometry on the formation of long-lived collision complexes.However, their investigation was limited to collisions of polyatomic molecules withhelium. The molecular interaction strength does not substantially increase withthe increasing complexity of molecules. At the same time, the strength of atom-molecule interactions is closely related with the polarizability of the atom. Thedynamics of more complex polyatomic molecule–polyatomic molecule collisionsis still completely unexplored. This leads to our third question: what is the role ofthe interaction strength and the geometry of polyatomic molecules on the formationof long-lived collision complexes at low collision energies?1.3 Challenges for scattering calculationsMy early PhD thesis work was motivated by the above three open questions. Dur-ing that work, I developed a love-hate relationship with cold molecules in a mag-netic field. On the one hand, the rich internal structure of molecules leads to manyinteresting discoveries in both experimental and theoretical studies as well as awide variety of applications. In addition, the use of external fields, together with4multiple scattering resonances, does allow for tuning collisional decoherence byvarying an external magnetic field near a Feshbach resonance, which will be dis-cussed in detail in Chapter 3. On the other hand, those unique properties of com-plex molecules are extremely hard to obtain from scattering calculations.The first challenge comes from the computational cost. The computational costof the time-independent quantum scattering calculations is proportional to the cubeof the number of basis functions used to represent the total wave functions. Evenfor simple atom-diatomic molecule collisions, the size of the basis set requiredfor converged results can be large, especially for systems with large masses andhighly anisotropic and deep potentials. The number of channels becomes muchgreater for molecule–molecule collision systems. The presence of external fieldsmakes the problem even more difficult. In the absence of external fields, the to-tal angular momentum J is conserved, and the Hamiltonian matrix has a block-diagonal structure, so the scattering calculations can be performed for each valueof J separately. The interaction of molecules with magnetic fields disturbs thesymmetry of the scattering problem. This makes quantum calculations of crosssections for molecule - molecule collisions in a magnetic field prohibitively diffi-cult. If the problem is formulated in the fully uncoupled basis set representation[42], fully converged calculations are impossible [49, 50]. For example, Janssenet al. [51] investigated cold and ultracold NH–NH collisions at zero field. Theircomputation based on a fully uncoupled basis set including up to 6 rotational statesand 7 partial waves, involved the numerical integration of 2382 coupled differen-tial equations. In the presence of an external magnetic field, the same calculation[49, 50] would require the integration of 38170 equations. In order to evaluatethe feasibility of evaporative cooling of certain 2Σ molecules in a magnetic field,we performed rigorous quantum calculations based on a total angular momentumrepresentation in the body-fixed frame, which allows for computations with muchlarger basis sets (Chapter 4). Where it comes to the third question involving atom-polyatomic molecule or polyatomic molecule-molecule collisions, the complexityof polyatomic molecules precludes full quantum calculations due to a large numberof degrees of freedom. Wherever possible, we used the classical trajectory methodto simulate the collisions of polyatomic molecules (Chapter 5).The other challenge is related to the uncertainty in the interaction potential sur-5faces. Despite the well-established methodology and affordable computation cost,the quantitative results from scattering calculations for atom-diatomic moleculesystems are very limited due to the lack of accurate ab initio potential energy sur-faces. Even when the accurate potential energy surfaces are available, binary colli-sions at low temperatures are very sensitive to details of the interaction potentials.A small variation of the interaction potential may result in large changes in thecollision cross sections. Many theoretical studies on the sensitivity of collision dy-namics to potential energy surfaces were performed by simply scaling the originalab initio potential energy surfaces, or equivalently the reduced mass of the collisioncomplex. It was shown that the shift of the resonance positions caused by the scal-ing may lead to dramatic change in the cross sections [37, 49–51]. Moreover, slightvariations of the interaction potential due to different fitting approaches may alsosignificantly alter the collision dynamics at ultracold temperatures [44]. There-fore, calculations based on a single potential surface would yield highly unreliableresults and conclusions. Compared with atom-molecule collisions, the strongerinteractions between polar molecules lead to denser resonance structures and con-sequently wider variability in collision outcomes due to uncertainties in potentialsurfaces.It is more meaningful to provide a probability distribution of collision out-comes by incorporating the uncertainty in the interaction potential surfaces intoscattering calculations for molecule-molecule collisions. That is, it would be desir-able to determine the distribution of collision outcomes η(V) of interest, such asthe elastic or inelastic scattering cross sections and elastic-to-inelastic cross sectionratios, when the the potential surface V is varied according to a certain distributionV ∼ pi(·). A Monte Carlo sampling would be a natural method for this purpose.Based on a Monte Carlo simulation, one can first draw a sufficiently large sampleof potential surfaces v1, v2, · · · , vn from the distribution of V∼ pi(V), and then carryout scattering calculations on each of those potential surfaces to obtain a sampleof collision outcomes η(v1), η(v2), · · · , η(vn). Using the law of large numbers,we can approximate the mean and variance of η(V) as the sample average and thesample variance of η(v1), η(v2), · · · , η(vn). In practice, the potential surface isgenerally defined by multiple terms in the angular expansion, so V is a functionof a high-dimensional random vector. For complex molecule-molecule collision6systems, the high-dimensional Monte Carlo sampling is very time-consuming. In-stead, in Chapter 5, we reduce the multi-dimensional random vector V to a scalarλ , and sample interaction potentials by simple scaling of an accurate ab initio po-tential surface, leading to a manifold of surfaces with different interaction strengthsand anisotropies.This scaling procedure increases the magnitude of both the lowest-order (isotropic)term and the higher-order terms representing the anisotropy of the potential, butpreserves the ratio of these terms. In order to find the distribution of collisionoutcomes nonambiguously covering the entire class of molecules considered, weneed a more sophisticated procedure involving an ensemble of potential surfacesgenerated by varying each term in the angular expansion independently. For com-plex systems, one can only afford calculations for a finite set of data points D ={(θ 1,η(θ 1)),(θ 2,η(θ 2)), · · · ,(θ n,η(θ n))}, where θ i represents any experimen-tally controllable parameters (such as the collision energy, temperature or externalfields) as well as the potential energy surface parameters. This leads to the lastquestions considered in this thesis: is it possible to make inference about the col-lision outcomes based on a finite set of model runs, or a more challenging one: isit possible to develop an analytical function of collision observables on the under-lying parameters by fitting a limited number of calculated results? If the answeris yes, previous problems, such as the computational cost and uncertainty analy-sis, can be resolved. In addition, this would provide a general solution to manyproblems encountered in the studies of collision dynamics of atoms or molecules,including parameter fitting by solving the inverse scattering problem, the predic-tion of collision properties of a specific molecular system based on the informationfor another molecule, the efficient calculation of thermally averaged observables.My later PhD thesis work was mainly concerned with the last questions. Tosolve this problem, I drew some great inspiration from statistics, and combinedstatistics and scattering calculations for improved predictions of molecular colli-sion observables. For any non-chaotic collision system, using a small number ofscattering calculations, we show that we can build a Gaussian Process model tostatistically approximate collision outcomes. The Gaussian Process model can belater used as an efficient surrogate for computationally intensive calculations inthe prediction of collision outcomes and uncertainty analysis. In addition, we ex-7plore the feasibility of using Gaussian process models to calibrate the scatteringcalculations in order to fit experimental observations or more complex scatteringcalculations.1.4 Thesis outlineChapter 2 describes the theoretical framework of rigorous quantum calculations. Itbriefly introduces the quantum mechanical formulation of the scattering problem.Chapter 3 demonstrates that the rate of collisional decoherence can be tuned byvarying an external magnetic field near a Feshbach resonance, and the collisionaldecoherence can be suppressed when the scattering lengths for the quantum statesin the coherent superposition are identical. In addition, a method for measuring thescattering length of ultracold particles is proposed.In Chapter 4, we perform rigorous quantum calculations for the cross sectionsin binary collisions of 2Σ molecules in a magnetic trap. Based on the randomsampling of the potential energy strengths, we provide the expectation intervalsfor the cross sections and the elastic-to-inelastic cross section ratios. We concludethat the majority of 2Σ molecules can be evaporatively cooled under conditionscorresponding to trapped molecular ensembles at T ∼ 10−3 K.In Chapter 5, we use classical trajectory calculations to study the effects of theinteraction strength and the geometry of rigid polyatomic molecules on the for-mation of long-lived collision complexes at low collision energies. We briefly de-scribes the basic ingredients of the molecular dynamics simulation of rigid molecule.We illustrates that the mean lifetimes of the collision complexes increase mono-tonically with the strength of the atom – molecule interaction, while the meanlifetimes of polyatomic molecule-molecule collision complexes are significantlyreduced due to non-ergodic effects.Chapter 6 provides the statistical background for the use of Gaussian Processmodels to emulate scattering calculations. In addition to interpolation, it discussesthe methodology of Bayesian calibration for fitting scattering calculations to realobservations in experiments.In Chapter 7, we demonstrate that a Gaussian Process model can be com-bined with a small number of scattering calculations to provide an accurate multi-8dimensional dependence of scattering observables on the collisions of polyatomicmolecules previously discussed in Chapter 5. We further explore possible applica-tions of Gaussian process model, such as calibration of our scattering calculationsto fit experimental observations or more complex scattering calculations.9Chapter 2Background materialThe purpose of this chapter is to summarize the details of quantum scattering the-ory for use in later chapters. We give a brief introduction to quantum scatteringtheory and present the general methodology for calculating non-reactive scatter-ing cross sections based on the work of Taylor [52], Arthurs and Dalgardo [53] andKrems, Stwalley and Friedrich [54]. We begin by discussing simple single-channelcollisions of two structureless particles, and the partial wave analysis to derive ex-pressions for the integral scattering cross sections. After that, we take into accountthe internal structure of the colliding particles, leading to a set of coupled differen-tial equations, and outline the numerical method for calculating elastic and inelasticcross sections in the multi-channel scattering events.2.1 Quantum mechanical formulationIn this thesis, we only deal with binary collisions. It is convenient to reduce therelative motion of two particles of mass mA and mB to the dynamical problem of avirtual particle with the reduced mass µ = mAmBmA+mB subject to an interaction potentialV (r,θ ,φ), where r,θ ,φ are the spherical polar coordinates.Before the scattering event, the particle can be considered to move along the zaxis with momentum h¯~k. The wave function of the particle is a plane waveψinc = Aeikz (2.1)10where k is the wave number, and A is a normalization factor. We will use atomicunits throughout this thesis. In atomic units, h¯ = 1, and α ≈ 1/137, where h¯ is thereduced Planck constant, and α is the fine-structure constant. In these units, theenergy of the particle is given byE = k22µ (2.2)The scattered part of the wave function is a spherical wave [55],ψsc = A f (θ ,φ)eikrr(2.3)where f (θ ,φ) is the scattering amplitude.2.2 Differential and integral cross sectionsWe assume that a uniform beam of particles with current density Jinc is directedalong the z-axis, and then scattered by a target fixed at the origin O to differentdirections with outgoing density Jsc(θ ,φ). To measure the strength of scattering ina given direction (θ ,φ), we introduce the solid angle dΩ, covering the ranges ofangles [θ ,θ +dθ ] and [φ ,φ +dφ ]. The area subtended by the solid angle isdS = r2 sinθdθdφ = r2dΩ (2.4)where r is the distance between the detector used for counting the scattered parti-cles and the target.The number of scattered particles passing through dS per unit time isdN = Jsc(θ ,φ)dS. (2.5)It is proportional to the incoming current density Jinc,dN = r2Jsc(θ ,φ)dΩ= Jincdσ (2.6)The ratio dσdΩ is the differential cross section, that is the cross section per unit solidangle.11The integral cross section is obtained by integrating the differential cross sec-tions over all solid angles,σ =∫4pi(dσdΩ)dΩ=∫4pir2(Jsc(θ ,φ)Jinc)dΩ (2.7)The integral cross section σ has the dimension of area, and represents the effectivearea of the incident beam scattered in all directions.The current density resulting from a wavefunction ψ is given by [56]J = 12iµ [ψ∗(∇ψ)− (∇ψ)∗ψ]= 1µ Im[ψ∗(∇ψ)] (2.8)At r→ ∞, the gradient operator can be approximated as∇= rˆ ∂∂ r (2.9)Inserting the expressions for ψinc in Eq (2.1) and ψsc in Eq (2.3) into Eq (2.8), wegetJinc = |A|2kµ (2.10)Jsc = |A|2kµr2 | f (θ ,φ)|2 (2.11)The replacement of Jinc and Jsc in Eq (2.7) by Eq (2.10) and Eq (2.11) leads to theexpression for the integral cross sections in terms of the scattering amplitudeσ =∫4pi| f (θ ,φ)|2dΩ (2.12)122.3 Single-channel scattering2.3.1 Partial wave analysisThe relative motion of two structureless particles interacting through a potentialV (R) is described by a time-independent Schro¨dinger equation:Hˆψ = Eψ (2.13)withHˆ =− 12µ∇2 +Vˆ (R) (2.14)where Hˆ is the total Hamiltonian, a sum of the kinetic energy operator and thepotential energy operator, ψ is the total wave function, R is defined by both theradial coordinates and the angular coordinates, and E = k22µ is the total energy ofthe colliding system.In spherical polar coordinates, the time-independent Schro¨dinger equation iswritten as [55][− 12µ1r2∂∂ r r2 ∂∂ r +Lˆ2(θ ,φ)2µr2 +Vˆ (r,θ ,φ)]ψ(k,r,θ ,φ) = Eψ(k,r,θ ,φ) (2.15)or[1r2∂∂ r r2 ∂∂ r −Lˆ2(θ ,φ)r2−2µVˆ (r,θ ,φ)+ k2]ψ(k,r,θ ,φ) = 0 (2.16)where Lˆ is the orbital angular momentum operator.Asymptotically, the total wave function of the colliding pair is a superpositionof the incoming wave function and the scattered wave function,ψ(k,R) R→∞−−−→ A(eikz + f (θ ,φ)eikRR)(2.17)13On the other hand, when Vˆ (r,θ ,φ) → 0, Eq (2.16) reduces to the Schro¨dingerequation for a free particle,[1r2∂∂ r r2 ∂∂ r −Lˆ2(θ ,φ)r2+ k2]ψ(k,r,θ ,φ) = 0 (2.18)The total wave function can be expanded in terms of spherical harmonics Ylm(θ ,φ)as [55]ψ(k,r,θ ,φ) =∑l∑mFlm(k,r)Ylm(θ ,φ) (2.19)For spherically symmetric and cylindrically symmetric interaction problems, thewave functions have no dependence on the azimuthal angle φ . In this case, thetotal wave function can be expanded in terms of Legendre polynomials [55],ψ(k,r,θ) =∑lFl(k,r)Pl(cosθ) (2.20)Each term in this expansion is called a partial wave, with l = 0,1,2, · · · referred toas the s, p,d, · · · waves.Substituting Eq (2.20) into Eq (2.18), and using the orthogonality property ofLegendre polynomials, we obtain the radial equation[d2dr2+ 2rddr+ k2− l(l +1)r2]Fl(k,r) = 0 (2.21)Replacing r with xk , we transform the radial equation into a spherical Bessel differ-ential equation[d2dx2+ 2xddx+(1− l(l +1)x2)]Fl(x) = 0 (2.22)There are two types of independent solutions to a spherical Bessel differential equa-tion: spherical Bessel function Jl(x) and spherical Neumann functions Nl(x). Thegeneral solution is a linear combination of Jl(x) and Nl(x) with unknown coeffi-14cients Bl and ClFl(x) = BlJl(x)+ClNl(x) (2.23)whereJl(x) x→∞−−−→1xsin(x− lpi/2) (2.24)Nl(x) x→∞−−−→ −1xcos(x− lpi/2) (2.25)Thus, the radial functions asymptotically becomeFl(x) x→∞−−−→1x[Bl sin(x− lpi/2)−Cl cos(x− lpi/2)](2.26)= Alxsin(x− lpi/2+δl) (2.27)whereAl =√|Bl|2 + |Cl|2 (2.28)δl =−arctan(Cl/Bl) (2.29)and δl is the scattering phase shift.2.3.2 Computation of integral cross sectionsThe radial function can be expressed in terms of the incoming (e−ix) and outgoing(eix) waves asFl(x) x→∞−−−→ Alx(−i)le−iδl2i(e−ix−Sleix)(2.30)where Sl is an S-matrix element,Sl = e2iδl (2.31)Writing the incoming wave function e−ikz in terms of the Legendre polynomials15in Eq (2.17), we obtainψ(k,R) R→∞−−−→ A[∑l(2l +1)il sin(kr− lpi/2)krPl(cosθ)+ f (θ)eikrr](2.32)By matching Eq (2.32) with the total wave function for free particles expanded inpartial wavesψ(k,R) R→∞−−−→∑lAl sin(kr− lpi/2+δl)rPl(cosθ) (2.33)we can extract the expression for the scattering amplitude f (θ)f (θ) = 12ik∑l(2l +1)(Sl−1)Pl(cosθ) (2.34)The integral cross section can be computed asσ =∑lσl =pik2 ∑l(2l +1)|Sl−1|2 (2.35)At very low collision energies, σl>0s vanish , and the integral cross section is de-termined by the s-wave scattering cross section σ0:σ0 =4pik2sin2 δ0(k) (2.36)We introduce another frequently used quantity – scattering length, which is definedas [54]a0 =− limk→0tanδ0k(2.37)The zero-energy cross section can be written asσ0 = 4pia20 (2.38)162.4 Multi-channel scatteringIn this section, we discuss collisions of particles with internal structure. Particleswith internal structure may undergo elastic or inelastic scattering.The total Hamiltonian of two particles (A and B) with internal structure is thesum of the Hamiltonian describing the relative motion of the colliding pair and theHamiltonian of the separated particlesHˆ =− 12µ1r2∂∂ r r2 ∂∂ r +Lˆ2(θ ,φ)2µr2 +Vˆ (r,τ )+ Hˆas(τ ) (2.39)where τ denotes collectively the internal coordinates, and Hˆas(τ ) represents theHamiltonian yielding the asymptotic states of the particles.The total wave function can be expanded in terms of the products of radialwave functions, spherical harmonics and the eigenfunctions of Hˆasψ(r,τ ) = 1r ∑α ∑l∑mFα,l,m(r)φα(τ )Ylm(θ ,ϕ)= 1r ∑iFi(r)Φi(τ ) (2.40)where i is a combination of three quantum numbers i ≡ (α, l,m), representing apossible scattering channel. Channels with threshold energies smaller than E arecalled open channels, otherwise they are closed channels. At ultracold tempera-tures, we consider transitions between open channels.Substituting Eq (2.40) into the Schro¨dinger equation Hˆψ(r,τ ) = Eψ(r,τ ), weobtain[∂ 2∂ r2 +2µE−Lˆ2r2−2µVˆ (r,τ )−2µHˆas(τ )]∑iFi(r)Φi((τ )) = 0 (2.41)Multiplying the resulting equation by Φ∗j yields a set of coupled-channel differen-tial equations:d2dr2Fj(r)+∑iW jiFi(r) = 0 (2.42)17whereW ji =∫τΦ∗j[2µE− Lˆ2r2−2µVˆ (r,τ ),−2µHˆas(τ )]Φi(τ )dτ= 〈Φ j|2µE−Lˆ2r2−2µVˆ (r,τ ),−2µHˆas(τ )|Φi〉 (2.43)Eq (2.42) can be written in matrix notation as[d2dr2+W(r,τ )]F(r) = 0 (2.44)where F(R) is a square matrix with each column being a linearly independent so-lution.Using numerical approaches, which will be discussed later, one can computethe S-matrix and the corresponding T -matrix. The matrix element Si j can be inter-preted as the probability amplitudes of the transitions from an incoming channeli≡ (α, l,m) to an outgoing channel j ≡ (α ′, l′,m′).The matrix T is defined byTi j = Si j−δi j (2.45)and the scattering amplitude can be written in terms of the T -matrix elements as[42]fα→α ′ =∑lm∑l′m′2pikαil+1−l′Y ∗lm(~k)Ylm(R′)Tαlm→α ′l′m′ (2.46)where kα = 2µ(E− εα), and ε is the eigenenergy of the asymptotic state α .Thus, the expressions for the differential and integral cross sections for thetransitions between the asymptotic states α and α ′ aredσα→α ′dΩ = | fα→α ′ |2 (2.47)18σ = pik2α∑lm∑l′m′|δαα ′δll′δmm′−Sαlm→α ′l′m′ |2= pik2α∑lm∑l′m′|Tαlm→α ′l′m′ |2 (2.48)2.5 Log-derivative method for scattering calculationsWe use an efficient algorithm for numerically solving the radial matrix Schro¨dingerequation based on the theory presented by Johnson [57] and Manolopoulos, Jamiesonand Pradhan [58].Instead of solving Eq (2.44) directly, we define a log-derivative matrixY(r) = F′(r)F−1(r) (2.49)Differentiating Y(x) and substituting the second derivative term with−W(r,τ )F(r),we obtain the first-order Ricatti equationY′(r)+W(r,τ )+Y2(r) = 0 (2.50)Within the classical forbidden region r→ 0, F(r)→ 0, leading to Y(r)→ ∞. Inthe calculations, we start from a point rmin deep inside the classically forbidden re-gion (where the total energy is less than the repulsive potential). At this point, weset Y(rmin) to a finite but very large number, such as Y(rmin) = 1030 I, and propa-gate it using the Johnson’s log-derivative propagator [57][58] with the appropriatestep size until the asymptotic region where the radial matrix function satisfies theboundary condition:F(r) r→∞−−−→ J(r)+N(r)K (2.51)J(r) and N(r) are diagonal matrices whose elements are made up of Riccati-Besselfunctions and modified spherical functions. The matrix K is related to Y(rmax) bydifferentiating the equation with respect to r and then multiplying it by its inverse,K =−[Y(rmax)N(rmax)−N′(rmax)]−1[Y(rmax)J(rmax)−J′(rmax)] (2.52)19The matrix K is partitioned into four blocksK =(Koo KocKco Kcc)(2.53)where Koo, Kcc, Koc and Kco corresponds open-open, closed-closed,open-closedand closed-open submatrices, respectively.The T-matrix can be expressed in terms of KooT = 2iKooI− iKoo(2.54)By extracting Ti j from the above matrix, we can compute the integral cross sectionsfor transitions between channels i and channels j.20Chapter 3Tuning collisional decoherence ofultracold molecules by a magneticfieldUltracold molecules in a superposition of two internal states are a promising plat-form for implementation of quantum computing and precision measurement. How-ever, collisional decoherence is a practical obstacle. In this chapter, we demonstratethat collisional decoherence of ultracold atoms or molecules in a coherent super-position of non-degenerate quantum states is suppressed when both the real andimaginary parts of the scattering lengths for the states in the coherent superposi-tion are equal. We show that the rate of collisional decoherence can be enhancedor suppressed by varying an external magnetic field near a Feshbach resonance.For some resonances, the suppression is very dramatic. We propose a method formeasuring the scattering length of ultracold particles in excited quantum states ex-hibiting Feshbach resonances.3.1 IntroductionQuantum coherence between internal states of atoms and molecules is exploitedfor many diverse applications, including Ramsey interferometry [59], quantumcomputing [60], quantum simulation of spin-lattice Hamiltonians [61] and exci-21ton energy transfer [62], coherent control of intra- and inter-molecular dynamics[63], magnetometry [64], laser-field alignment of molecules [65] and stereochem-istry [66]. These applications are always limited by decoherence. Decoherencearises from random forces exerted on atoms or molecules in coherent superposi-tion states. Decoherence underlies the quantum-classical system transition [67]and understanding decoherence mechanisms is key to unraveling the relationshipbetween quantum and classical mechanics. Molecular collisions are an importantsource of decoherence in experiments with atoms and molecules in the gas phase.Due to the stochastic nature of molecular collisions, collisional decoherence is usu-ally uncontrollable. This limits the experimental studies of collisional decoherenceand other decoherence sources, especially when several decoherence mechanismshave similar timescales. On the other hand, collisional decoherence can be used asa probe of interparticle interactions in a gas [68, 69].There is currently great interest in experimental work with atoms and moleculesat ultralow temperatures. Cooling molecular gases to subKelvin temperatures elim-inates many complications of the experimental measurements at ambient temper-atures and opens up the possibility of controlling intermolecular interactions [13,54]. This suggests a possibility of designing a system with tunable collisional deco-herence. The most widely used technique for controlling interparticle interactionsin ultracold atomic gases is based on tuning a Feshbach resonance with an exter-nal magnetic field [70]. The collision properties of ultracold atoms in a particularquantum state are determined by two parameters: the real and imaginary parts ofthe scattering length. Feshbach resonances modify the scattering length, allowingfor the magnetic field control of atomic collisions [70]. In the present work, weconsider an ensemble of ultracold particles (atoms or molecules) prepared in a co-herent superposition of two quantum states and explore the possibility of using aFeshbach resonance to suppress collisional decoherence. With two states in thecoherent superposition, the decoherence is determined by the real and imaginaryparts of two scattering lengths. We present numerical calculations that demonstratethe possibility of tuning these four parameters to the regime where the damping ofcoherent oscillations must be dramatically suppressed.The second problem we consider here is the possibility of using an experimen-tal measurement of collisional decoherence for determining the scattering length of22ultracold atoms or molecules. The scattering length of ultracold atoms can be deter-mined by measuring the thermalization rate of an ultracold gas driven out of ther-mal equilibrium. This method can only be used for atoms in the absolute groundstate immune to inelastic collisions. An alternative approach is based on mappingout the positions and widths of Feshbach resonances by measuring three-body re-combination and fitting the experimental data by a model based on multichannelscattering calculations [70]. This method can provide the scattering lengths ofatoms in different internal states, but the extension of this approach to molecularsystems is anticipated to be extremely difficult. Clearly, there is a need for the de-velopment of a versatile experimental technique for measuring the scattering lengthof ultracold atoms and molecules. We present an analysis demonstrating that, if thescattering length for any state is known, the magnetic field dependence of the scat-tering length of any other state can be determined by measuring the depletion ofcoherence.3.2 Methodology3.2.1 Denisty matrix and decoherenceThe density matrix is a useful tool for the description of decoherence [71]. Wewould like to illustrate this concept using a two-level system. If a system is origi-nally prepared in a coherent superposition of states |a〉 and |b〉|ψ〉= c1|a〉+ c2|b〉the corresponding density matrix written in the basis {|a〉, |b〉} reads asρˆ = |ψ〉〈ψ|= c1c∗1|a〉〈a|+ c1c∗2|a〉〈b|+ c2c∗1|b〉〈a|+ c2c∗2|b〉〈b|In the matrix representation:ρˆ =(c1c∗1 c1c∗2c2c∗1 c2c∗2)23The off-diagonal terms in the matrix representation describe the quantum coher-ence between states |a〉 and |b〉. The decoherence is the decay of the off-diagonalterms. When there is complete loss of coherence between the two states, all off-diagonal matrix elements vanish, and the initial pure state degrades to a statisticalmixtureρˆ = c1c∗1|a〉〈a|+ c2c∗2|b〉〈b|3.2.2 Theory of collisional decoherenceThe microscopic theory of collisional decoherence due to molecular collisions isbased on the work presented by Hemming and Krems [72], and it is general forany atomic or molecular system. For concreteness, we consider an ensemble ofultracold molecules prepared in a coherent superposition of two non-degeneratestates, labeled a and b, and placed in an ultracold buffer gas of structureless atoms,such as Mg(1S). The buffer gas represents the environment that destroys coher-ence between molecular states. The coupling to the environment is enabled by theatom - molecule collisions. We assume that the buffer-gas is dilute enough to elim-inate three or more-body collisions, and the minimum energy difference betweeninternal states is larger than the trap depth. As a consequence, inelastic collisionscan release enough energy to eject both the molecules and the buffer-gas atomsinvolved, leading to the depletion of the molecular ensemble.Collisions can entangle the states of molecules with the buffer-gas atoms. How-ever, we are only interested in describing the ensemble of molecules. The reducedone-particle density operator ρˆ ′ for the molecular internal states after collisions isobtained by tracing over the translational degrees of freedom of the buffer gas atomρˆ ′ = trgas{ρˆ pair′} (3.1)where ρˆ pair′ is the density operator decribing the atom-molecule mixed states aftercollisions.Since the incoming and outgoing scattering wave functions are related by a24two-particle scattering operator Sˆ [52]|ψout〉= Sˆ|ψin〉, (3.2)The density operator ρˆ pair′ after collisions is given byρˆ pair′ = Sˆρˆ pairSˆ†= Sˆ(ρˆ⊗ ρˆgas)Sˆ† (3.3)where ρˆ and ρˆgas are the initial density operators for the molecules and buffer gasatoms, respectively.Eq (3.1) can be rewritten asρˆ ′ = trgas{Sˆ(ρˆ⊗ ρˆgas)Sˆ†} (3.4)Following Ref. [72], we describe decoherence by two measurable quantities, thedecay of the off-diagonal elements of the density matrix (the total coherence signal)|ρela,b(t)| =√ρela,b(t)ρela,b(t) (3.5)and the relative coherence describing the coherence of molecules remaining in thegas after inelastic collisions,ηab =(ρela,bρelb,aρela,aρelb,b)1/2, (3.6)where the matrix elements ρeli j are defined asρeli j =∫trgas[Sˆii(ρˆi j⊗ ρˆgas)Sˆ†j j]d3P, (3.7)Sˆii = 〈i|S|i〉 is an operator acting on the spaces for the relative translational motionof the colliding pair and the motion of the center of mass, ρˆi j denotes the matrix25elements 〈i|ρˆ| j〉 of the initial density operator for the molecules, and the integrationis over the translational momenta of the molecules. At first sight, it may appear thatthe quantities defined by Eqs. (3.5) and (3.6) are unaffected by inelastic collisions.However, due to the unitarity property of the Sˆ operator, the presence of inelasticcollisions divides the Sˆ operator into elastic and inelastic parts,Sˆ = Sˆel + Sˆin= ∑iSˆii +∑i 6= jSˆi j (3.8)Thus, the rates of inelastic scattering affect Sˆii’s and hence the density matrix ele-ments defined by Eq. (3.7).The time dependence of ηab is given by [72]ηab(t) = ηab(t = 0)[1+T 1/2tξ1 +T(ξ2,1t +t22ξ2,2)+ · · ·], (3.9)where T is the temperature of the buffer gas and ξ are the expansion coefficientsdefined below. The decoherence rate isddtηab(t) = ηab(0)[T 1/2ξ1 +T (ξ2,1 + tξ2,2)+ · · ·]. (3.10)At low temperatures and short times, the decoherence rate is thus determined bythe coefficient ξ1, while the coefficients ξ2,1 and ξ2,2 may become important atlarger times and higher temperatures. The s-wave scattering amplitude of ultracoldmolecules in state i ( fii) can be expanded in powers of linear momentum (p = h¯k)as follows:fii(p)≈−ai +bi p+ ci p2 (3.11)Ref. [72] presents the expressions for the coefficients ξ in terms of the expansioncoefficients ai and bi of the scattering amplitude. It is more convenient for thepresent study to express the coefficients ξ in terms of the real and imaginary partsof the scattering length for molecules in states a and b. To do this, according to ref.26[73], we rewrite the s-wave scattering amplitudes fii asfii(p) =1gi(p2)− ip(3.12)where gi(p2) is an analytical even function of the momentum p.A Taylor series expansion of gi(p2) has the formgi(p2) =−1ai+ 12r0 p2 + · · · (3.13)where ai is the scattering length and r0 is the effective range [73]. In the presenceof inelastic collisions, the scattering length is a complexai = (αi− iβi) (3.14)Substituting Eq (3.13) and Eq (3.14) into Eq (3.12), we expand fii(p) in powersof p asfii(p)≈−ai +[2αiβi +(α2i −β 2i)i]p+ ci p2 (3.15)By matching Eq (3.11) with Eq (3.15), we can express the coefficients ξ1, ξ2,1 andξ2,2 as follows:ξ1 = −25/2pi1/2ngask1/2Bm∗1/2[(αb−αa)2 +(βb−βa)2], (3.16)ξ2,1 = 12pingaskBr3/2(βb +βa)[(αb−αa)2 +(βb−βa)2]/h¯, (3.17)ξ2,2 ={32pi3n2gaskBm∗+8pikBn2gasm[3(2r+1)1/2 + 1+2r+3r2rsin−1(rr+1)−4(1+ r)]}×|ab−aa|4−64pikBn2gasm[3(2r+1)1/2 + 1+2r+3r2rsin−1(rr+1)−4(1+ r)]×(β 2b +β 2a +βbβa)|ab−aa|2. (3.18)Here, ngas is the density of the buffer gas, kB is the Boltzmann constant, m is the27mass of the buffer gas atom, m∗ is the reduced mass of the atom - molecule collid-ing pair, r is the ratio of the buffer gas atom mass and the molecule mass, αi andβi denote the real and imaginary parts, correspondingly, of the scattering length aidescribing ultracold collisions of molecules in state i with buffer gas atoms.Eqs. (3.16) - (3.18) lead to two observations: (i) The decoherence rate vanishesif both the real and the imaginary parts of the scattering lengths for the two states aand b have the same magnitude and sign. (ii) If one of the states a or b is the abso-lute ground state, the corresponding scattering length is real and the decoherencerate can never vanish, even if the real parts of the two scattering lengths are equal.This means that the decoherence of coherent superpositions of states a and b canbe suppressed to zero only if both of the molecular states are excited states.3.2.3 Complex scattering length from S-matrixAccording to Eqs. (3.16) - (3.18), the problem for evaluating the relative coherenceηab and the decoherence rate reduces to finding the the real and the imaginary partsof the scattering lengths for the two states a and b in the superposition. The diag-onal S-matrix element in the incoming channel i is defined in terms of a complexphase shift δi with a positive imaginary partSii = e2iδi (3.19)Since a scattering length is the limit of− tanδ(k)k as k→ 0, we can relate Sii with thescattering length ai bySii =1− iaik1+ iaik(3.20)Since Sii is also complex, it is convenient to split it into its real and imaginarypart as Sii = Re(Sii)+ iIm(Sii) and express the real and the imaginary parts of thescattering length as follows:αi =−2Im(Sii)k [(1+Re(Sii))2 + Im(Sii)2](3.21)βi =1−Re(Sii)2− Im(Sii)2k [(1+Re(Sii))2 + Im(Sii)2](3.22)28In the numerical calculation, we first calculate the K-matrix by using the bound-ary conditions for asymptotic radial wave functions, and then express the real andimaginary parts of the T -matrix in terms of open-open submatrix of the K-matrixKoo by using Eqs (2.54) as:Im(T) = 2Koo(I+K2oo)−1 (3.23)Re(T) = −2K2oo(I+K2oo)−1 (3.24)Finally, we can obtain the real and imaginary parts of the S-matrix by using therelation between the T -matrix and the S-matrix (Eq (2.45))3.2.4 Collisions of Mg(1S) with NH (3Σ) in a magnetic fieldThe total Hamiltonian of the Mg–NH system in a magnetic field can be representedas Eq(2.39) with the asymptotic Hamiltonian of the form:Hˆas = Hˆrot + Hˆspin−rot + Hˆspin−spin + HˆZeeman, (3.25)where Hˆrot, Hˆspin−rot, Hˆspin−spin and HˆZeeman are the rotational, spin-rotation, spin-spin and Zeeman Hamiltonians of the 3Σ-molecule, respectively.Following the formalism for quantum-mechanical theory of molecular scatter-ing in an external field developed by Krems and Dalgarno [42], we expand thechannel wave functions φi in a set of space-fixed uncoupled basis functions.φi = |NMN〉|SMS〉|LML〉 (3.26)where N and S are the rotational and spin angular momenta of the 3Σ-molecule, Lis the orbital angular momentum describing the end-over-end rotation, MN , MS andML are the projections of N, S and L on the quantization axis, which is set to be theZ-axis, and i≡ (L,ML,N,MN ,S,MS) is a set of the six quantum numbers.The matrices of Lˆ2, Hˆrot and HˆZeeman are diagonal in basis (3.26), and are given29by following expressions:〈φi| Lˆ2 |φi′〉 = δii′ L(L+1) (3.27)〈φi| Hˆrot |φi′〉 = δii′ BeN(N +1) (3.28)〈φi| HˆZeeman |φi′〉 = δii′ 2µ0BMS (3.29)where Be is the rotation constant of the molecule, µ0 is the Bohr magneton, and Bis the magnitude of a magnetic field.The diagonal and off-diagonal matrix elements of Hˆspin−rot have the followingexpressions:〈φi| Hˆspin−rot |φi〉 = γMNMS (3.30)〈φi| Hˆspin−rot |φi′〉 = δL,L′δML,M′LδN,N′δMN ,M′N±1δMS,M′S∓1 γCii′ (3.31)where Cii′ = 12 [N(N+1)−M′N(M′N±1)]1/2[S(S+1)−M′S(M′S∓1)]1/2, and γ is thespin-rotation constant of the molecule.The matrix elements of Hˆspin−spin are defined as〈φi| Hˆspin−spin |φi′〉= δL,L′δML,M′LλSS(2/3)(4pi/5)1/2√6∑q(−1)q〈NMN |Y2,−q(rˆ)|N′M′N〉〈SMS|[S⊗S](2)q |SM′S〉= δL,L′δML,M′LλSS(2/3)√30∑q(−1)q(−1)−MN [(2N +1)(2N′+1)]1/2×(N 2 N′−MN −q M′N)(N 2 N′0 0 0)(−1)S−MS(S 2 S−MS q M′S)(3.32)where λSS is the spin-spin constant of the molecule.The atom-molecule interaction potential is expanded in spherical harmonics,and its matrix elements in the basis (3.26) are evaluated using the Wigner-Eckart30theorem [74]〈φi|Vˆ (r,τ) |φi′〉= δMS,M′S ×∑λ(L λ L′0 0 0)(N λ N′0 0 0)[(2L+1)(2L′+1)(2N +1)(2N′+1)]1/2×∑Mλ(−1)Mλ−ML−MN(L λ L′−ML −Mλ M′L)(N λ N′−MN Mλ M′N)(3.33)3.3 ResultsThe scattering length of ultracold atoms and molecules can be tuned by an exter-nal magnetic field in the presence of scattering Feshbach resonances [13]. Near theresonances, both the real and imaginary parts of the scattering length undergo rapidvariation as functions of the external field. Atoms and molecules in different inter-nal states usually exhibit scattering resonances at significantly different magnitudesof external fields. If the scattering lengths of the states a and b were both real, itwould be generally easy to find an interval of an external field where the scatteringlength of state a would vary rapidly and that of state b would be independent ofthe external field. It would therefore be easy to find an external field magnitude, atwhich the two scattering lengths would be equal and the decoherence rate wouldvanish. However, it is not clear if the two scattering lengths can be made the same,if they are both complex.To explore this, we calculate the decoherence rate for NH molecules preparedin a coherent superposition of different rotational states in a buffer gas of ultracoldMg atoms. The experimental work on cooling NH molecules using buffer gascooling is currently actively pursued [26]. It was also shown theoretically that NHmolecules can be effectively cooled if placed in an ultracold gas of Mg atoms [37].The cross sections for Mg - NH collisions at ultracold temperatures exhibit multipleFeshbach resonances [37]. We use the following parameters for our calculations:T = 10−6 K, ngas = 1011 cm−3, m = 24.3 g/mol, and r = 1.62. The scatteringlengths are computed using the Mg - NH interaction potential presented in Ref.[38], which is generated by using highly correlated ab intio electronic structure31calculation. The scattering calculations are performed as described in Ref. [42].The energy levels of the NH molecule are separated into manifolds that canbe labeled by a (nearly conserved) quantum number N of the rotational angularmomentum. In the presence of a magnetic field, each N manifold is split intoseveral fine structure and Zeeman sublevels (see Figure 3.1).Within each manifold, the energy levels can be labeled by the quantum numberof the total angular momentum ( j) and its projection on the magnetic field axis(m j). The total angular momentum is a good quantum number at zero magneticfield. The label j is used for states that adiabatically correlate with a particularj-level at zero field. We consider NH molecules prepared in a coherent superpo-sition of the state ( j = 4,m j = −4) in the N = 3 manifold (state a) and the state( j = 3,m j = −3) in the N = 2 manifold (state b). We assume that the moleculesare initially in a pure state so that ηa,b(t = 0) = 1. Figure 3.2 presents the scatter-ing lengths for ultracold collisions of Mg atoms with NH molecules in these states.Figure 3.3 shows that both the real and imaginary parts of the scattering lengths be-come similar at certain magnitudes of the magnetic field and demonstrates that thedecoherence rate is dramatically suppressed at certain magnitudes of the magneticfield. The calculation shows that the decoherence rate can be modified over therange spanning six orders of magnitude by varying the magnetic field. In order toverify that the suppression of collisional decoherence observed in Figure 3.3 is notcoincidental, we repeated the calculations for several different pairs of molecularstates, labeled in Figure 3.1. For each combination of states, there is an intervalof the magnetic field, where the decoherence rate is suppressed by more than threeorders of magnitude.Eqs. (3.16) - (3.18) show that the decoherence rate is determined by the differ-ence of the scattering lengths in states a and b. This suggests that a measurementof the decoherence rates can be used to determine the scattering length of state b ifthe scattering length of state a is known. To do this, it is necessary to measure thetime evolution of both |ρela,b(t)| and ηa,b(t). At low temperatures, the time evolutionof |ρela,b(t)| is given by [72]|ρela,b(t)| =√ρela,b(t)ρelb,a(t) = |ρa,b(0)|eζ0t[1+T 1/2ζ1t + · · ·], (3.34)32193194195196197198199969798991000 0.4 0.8 1.2 1.6 2.0Magnetic Field (T)303132333435Energy (cm-1)state astate bstate dstate c321243120 -4 3-3-2mjN=3jN=2N=1state e 4Figure 3.1: Energy levels of the NH molecule in a magnetic field.33-20020α (a.u.)state bstate astate d0.4 0.5 0.6 0.7 0.8 0.9 1Magnetic Field (T)0204060β (a.u.)Figure 3.2: The real (upper panel) and imaginary (lower panel) parts of theMg – NH s-wave scattering lengths: solid curves – state a; dashedcurves – state b; dotted curves – state d of Figure 3.1.3410-410-2100102104|aa-ab|2  (a.u.)ρ aa(t) / ρ aa(0)ρ bb(t) / ρ bb(0)|ρ ab(t)| /  |ρ ab(0)| Field (T)10-1010-810-610-410-2Decoherence Rate ( s-1)012345Time (s) 3.3: Upper left panel: The absolute difference of the scattering lengthsfor NH molecules in states a and b (labeled in Figure 3.1) as a functionof the magnetic field. Lower left panel: The t = 0 decoherence rateof NH molecules prepared in the superposition of states a and b as afunction of the magnetic field. Right panels: The time dependence ofthe density matrix elements ρela,a (dashed lines), ρelb,b (solid lines) and|ρela,b| (dotted lines) at two magnitudes of the magnetic field.35whereζ0 = −4pi h¯ngasm∗(βa +βb)2, (3.35)ζ1 =25/2pi1/2k1/2B ngasm∗1/2[(βa +βb)2− (αa−αb)2]. (3.36)In the limit of low T , the decay of |ρela,b(t)| is entirely determined by Eq. (11),which provides the relation between βa and βb. The magnitudes of βa and βb canalso be determined directly by measuring the decay of the populations of states aand b. To obtain αb from αa, one can use the decoherence rate (3.10)αb = αa±√RabCT 1/2− (βb−βa)2 (3.37)where Rab = dηab(t→ 0)/dt and C = 25/2pi1/2ngask1/2B /m∗1/2. Unfortunately, thisequation is not unambiguous and provides two possible values for αb. This prob-lem arises because the decoherence rate is determined by the absolute differenceof the scattering lengths.To overcome this problem, we propose first to create a coherent superpositionof state a and another state d, whose scattering length has no resonances and isindependent of the magnetic field B in a given interval of magnetic fields ∆B, wherestate a exhibits, at least, one Feshbach resonance. Most atomic and molecularsystems possess such states. In ultracold alkali metal atoms, hyperfine states withthe largest value of the total angular momentum f and the angular momentumprojection m f = f are resonance-free. In the molecular system considered here, themaximally stretched state labeled d in Figure 3.1 exhibits no resonances (see Figure3.2). The real part of the scattering length αd is also, in principle, undeterminedαd =αa±√Rad/CT 1/2− (βa−βd)2. However, because αa is an analytic functionof the magnetic field [75, 76] and Rad is a continuous function of the magnetic field,only one of the signs gives the value of αd that is independent of the magneticfield. This is illustrated in Figure 3.4. Thus, given αa(B) and the decoherencerate Rad(B), it is possible to determine αd non-ambiguously. Once αd and βd areknown, the scattering length of any other state x exhibiting a Feshbach resonancecan be determined non-ambigously by measuring the decoherence rate Rax as a360.4 0.5 0.6 0.7 0.8 0.9 1Magnetic Field (T)-30-20-1001020Q (a.u.)Figure 3.4: The quantity Q = αa±√Rad/CT 1/2− (βa−βd)2: solid curve –minus sign; dashed curve – plus sign.function of the magnetic field.In conclusion of this section, we need to point out the challenges of measuringthe decoherence rate at ultralow temperatures. We note that inelastic collisions leadto the decay of all the density matrix elements ρela,a, ρelb,b and |ρela,b| due to trap losseven if the scattering lengths of the states a and b are equal and there is no decoher-ence of trapped molecules. As demonstrated in Figure 3.3, when the decoherencerate dηab/dt is suppressed, the matrix elements ρela,a, ρelb,b and |ρela,b| decay at thesame rate. This decay limits the lifetime of the experiments. The decay rates of thediagonal matrix elements are determined by the imaginary parts of the scatteringlengths. The decay rate of the off-diagonal matrix element |ρela,b| is determined byboth the real and imaginary parts of the scattering length. Eqs. 3.10 and 3.34 showthat the effect of elastic scattering (the real parts of the scattering length) is con-37tained in the coefficients proportional to T 1/2 or higher powers of temperature. Inorder to measure the effect of elastic scattering on collisional decoherence, it maytherefore be necessary to raise the temperature of the gas. However, the collisionsystem must remain in the s-wave scattering regime. The optimal temperature ofthe experiments must therefore be just below the temperature, at which the p-wavecontribution to the scattering cross sections becomes important. This temperatureis different for different collision systems and depends on the reduced mass ofthe colliding particles, the strength of the inter-particle interaction, the strength ofs-wave scattering and the magnitude of the resonant enhancement of the s-wavescattering cross sections.3.4 ConclusionIn this chapter, we have answered the first question raised in Chapter 1 by demon-strating that collisional decoherence of ultracold molecules prepared in coherentsuperpositions of internal states can be controlled by tuning an external magneticfield near a scattering resonance. Our results show that the decoherence rate issuppressed when both the real and imaginary parts of the scattering lengths for thestates in the coherent superposition are equal. The scattering resonance modifiesthe scattering length of one of the states in the coherent superposition, which canbe exploited to reduce or enhance the rate of collisional decoherence by varyingthe external field.The ability to tune collisional decoherence in a molecular gas by varying anexternal field can be used for numerous practical and fundamental applications.We have demonstrated that measuring the decay of the total and relative coherencesignals can be used to determine the scattering length of molecules in any quantumstate x, if the scattering length of one state is known. The proposed method relieson measuring the decay of coherence between a state that exhibits resonances anda state that is resonance-free. Tuning collisional decoherence and measuring thecoherence times can also be exploited for calibrating other less controllable sourcesof decoherence.38Chapter 4Elastic and inelastic collisions of2Σ molecules in a magnetic fieldIn this chapter, we calculate the cross sections for elastic scattering and Zeeman re-laxation in binary collisions of molecules in the ro-vibrational ground state of a 2Σelectronic state and the Zeeman state with the electron spin projection MS = 1/2on the magnetic field axis. This is the lowest-energy state of 2Σ molecules con-fined in a magnetic trap. To take into account the uncertainty in the potential en-ergy surface, we average calculation results over multiple molecule - moleculeinteraction potentials, and provide the expectation intervals for the cross sectionsand the elastic-to-inelastic cross section ratios. We find that the elastic-to-inelasticcross section ratios under conditions corresponding to trapped molecular ensem-bles at T ∼ 10−3 K exceed 100 for the majority of 2Σ molecules. The range of 2Σmolecules expected to be collisionally unstable in magnetic traps at T < 10−3 Kis limited to molecules with the spin-rotation interaction constant γSR > 0.5 cm−1and the rotational constant Be < 4 cm−1.4.1 IntroductionA major focus of current research in molecular physics is the preparation of ultra-cold molecules [13, 54]. While several experimental techniques have been devel-oped for cooling molecules to temperatures of about 10−3 Kelvin [13, 22, 54, 77], it39is desirable to obtain dense molecular samples at temperatures below 10−6 Kelvin[13]. Molecules at such ultracold temperatures can be produced by photoassocia-tion of ultracold atoms [20, 78]. However, this technique is limited to molecularspecies composed of atoms amenable to laser cooling. A number of promising ap-plications – such as the development of quantum simulators of condensed-matterHamiltonians [61], precision measurements of fundamental constants [79], andcontrolled chemistry [80] – need a greater variety of ultracold molecules. There-fore, it is necessary to develop methods for cooling molecules from 10−3 Kelvin toultracold temperatures. Despite this goal, bridging the milliKelvin – microKelvintemperature gap remains a major challenge.Once cooled to temperatures of T ∼ 1− 10 mK, molecules can be confinedin an external field trap. Further cooling of a molecular ensemble can potentiallybe achieved by evaporation of the most energetic molecules, as is commonly donewith alkali metal atoms in magnetic traps [24]. The evaporative cooling relies onenergy thermalization enabled by momentum transfer in molecule - molecule col-lisions conserving the internal energy of the colliding species. If molecules aretrapped in excited states, the evaporative cooling process is hindered by inelasticcollisions which release energy and expunge the colliding species from the trap.For example, dc magnetic or electric traps confine molecules in excited Zeeman orStark states that can undergo collisional energy relaxation [81], removing particlesfrom the most dense region of the trap. This leads to loss of low-energy moleculesand results in heating. The prospect for evaporative cooling of molecules in mag-netic traps thus hinges on the relative probability of elastic and inelastic scatteringin molecule - molecule collisions in a magnetic field.The experiments on magnetic trapping of CaH(2Σ+) [22] CrH(6Σ+) [82], MnH(7Σ+)[82] and NH(3Σ−) [26, 27] in a cold buffer gas of He atoms stimulated multiple the-oretical studies of atom - molecule collisions in a magnetic field [42–44, 76, 83–87]. It was found that Zeeman relaxation in atom - molecule collisions is inducedby an interplay of intra-molecular spin-dependent interactions and the anisotropyof the atom - molecule interaction potentials [42–44]. A series of experimental[22, 25, 84] and theoretical [43, 44, 76] studies showed that collisions of mag-netically trapped molecules with weakly interacting He atoms are predominantlyelastic, as a consequence of weak interaction anisotropy. At the same time, it was40demonstrated that the rates of Zeeman relaxation in collisions of molecules withmore polarizable atoms characterized by a stronger interaction anisotropy are muchgreater [36, 37, 76, 85]. The interaction anisotropy in molecule - molecule colli-sions is generally much larger than in atom - molecule collisions. This raises thequestion, can molecules be evaporatively cooled?The evaporative cooling of molecules in a magnetic trap is certainly possibleunder specific conditions. Recently, Stuhl et al. [40] demonstrated the evaporativecooling of an ensemble of OH(2Π) molecules confined in a magnetic trap from 51to 5.1 milliKelvin. It was shown that inelastic Zeeman relaxation in OH-OH col-lisions can be suppressed in this temperature interval due to the peculiar structureof molecules in the 2Π electronic state preventing cold molecules from excursionsinto the strongly interacting regime. The rates of Zeeman relaxation must also besuppressed in the limit of very low temperatures and weak magnetic fields. Thissuppression is a result of the centrifugal barriers which reduce the collision flux intothe inelastic scattering channels [83, 88]. As a result of this mechanism, moleculescooled to 10−6 Kelvin or below should be amenable to evaporative cooling in shal-low magnetic traps [89]. However, in order to trap molecules at 10−3 Kelvin, it isnecessary to apply magnetic fields of finite strength (& 50 Gauss). The possibil-ity of evaporative cooling of molecules at such magnetic fields and temperaturesremains an open question.To answer this question, it is necessary to perform rigorous quantum calcula-tions of molecular scattering properties in the presence of a magnetic field. Thetheoretical analysis of cold molecular collisions in a magnetic field is, however,complicated by two factors [49, 50, 89, 90]. First, molecule - molecule collisionsat low temperatures are affected by scattering resonances, which are sensitive todetails of intermolecular interaction potentials. A small variation of the interac-tion potential may therefore lead to large changes of the collision cross sections.Second, the interaction of molecules with magnetic fields disturbs the symmetryof the scattering problem. This makes quantum calculations of cross sections formolecule - molecule collisions in a magnetic field prohibitively difficult. If theproblem is formulated in the fully uncoupled basis set representation [42], fullyconverged calculations are impossible [49, 50]. As a result, there are very limitedreliable theoretical data for molecule - molecule scattering cross sections in a mag-41netic field. In particular, nothing is known about the relative probabilities of elasticand inelastic scattering of 2Σ radicals, the simplest kind of open-shell moleculeswith a non-zero magnetic moment.The goal of the present study is to provide predictions of the limits of elas-tic and inelastic scattering cross sections of 2Σ molecules in a magnetic field. Toreduce the uncertainty stemming from the limited accuracy of the interaction po-tentials and scattering resonances, we average the results of quantum scatteringcalculations over data obtained with multiple interaction surfaces. This producesstatistical expectation intervals of cross sections. We employ the time-independentquantum scattering formalism using the total angular momentum representationthat allows large basis sets and dramatically reduces the basis set truncation er-ror [90]. The calculations are performed as functions of the molecular structureparameters in order to provide useful reference points for future experiments.4.2 Computation detailsWe consider collisions of molecules prepared in the ro-vibrational ground state of a2Σ electronic state and the Zeeman state with the electron spin projection MS = 1/2on the magnetic field axis. This is the lowest-energy state of 2Σ molecules con-fined in a magnetic trap [22]. Evaporative cooling relies on momentum trans-fer in molecule - molecule collisions that preserve the spin alignment. How-ever, molecule - molecule collisions may lead to inelastic relaxation populatingthe lower-energy Zeeman state characterized by MS = −1/2. Molecules in thestate with MS = −1/2 are expunged from the trap by the magnetic field gradi-ent. The ratio of the cross sections for elastic collisions and the Zeeman relaxationMS = 1/2→MS =−1/2 is thus a critical parameter determining the possibility ofevaporative cooling of 2Σ molecules in a magnetic trap.We use the numerical technique developed by Tscherbul and Dalgarno [91] andTscherbul, Suleimanov and Krems [90]. The details of the theory were describedpreviously [42, 90, 91]. Traditionally, the time-independent quantum calculationsof molecular scattering properties were based on representing the total wave func-tion of the collision complex by a basis set expansion in terms of the eigenstatesof Jˆ2 and JˆZ , where Jˆ is the total angular momentum of the collision complex and42JˆZ is the Z-component of Jˆ in the space-fixed coordinate frame. The substitutionof this expansion in the time-independent Schro¨dinger equation leads to a systemof coupled differential equations for the expansion coefficients [53]. In the ab-sence of external fields, the total angular momentum is conserved and the totalangular momentum basis leads to a significant reduction of the number of coupledequations. The interaction of molecules with an external field induces couplingsbetween states of different total angular momenta. Therefore, it was originally sug-gested that the collision theory of molecules in external fields is best formulatedusing a fully uncoupled space-fixed basis representation [42]. The uncoupled basissimplifies the evaluation of the Hamiltonian matrix elements so it was adopted bymany authors [92–95].If evaluated in the total angular momentum basis, the matrix of the molecule- field interactions is tridiagonal. Tscherbul and Dalgarno [91] and Tscherbul [96]showed that the tridiagonal character of the Hamiltonian matrix can be exploitedfor quantum scattering calculations of molecular properties in fields. As describedin Refs. [90, 91, 96], the scattering cross sections for molecules in external fieldscan be calculated using the total angular momentum representation with multipleJ states included in the basis. The convergence of the scattering cross sectionscan then be sought with respect to the number of J states. The results of Refs.[90, 91] show that a few coupled J states are usually sufficient for full convergenceof low-energy collision properties in a magnetic field, which makes the total an-gular momentum basis much more efficient than the fully uncoupled basis. If thetotal angular momentum basis is used, the collision theory can be formulated usingeither a space-fixed or body-fixed representation [90]. Here, we use the body-fixedformulation described in detail in Refs. [90, 91].We assume that the space-fixed (SF) quantization axis Z is directed along themagnetic field vector and the body-fixed (BF) quantization axis z is directed alongthe vector R joining the centers of mass of molecules A and B. The Hamiltonianfor two identical 2Σ molecules in a magnetic field can be written (in atomic units)asHˆ =− 12µR∂ 2∂R2 R+lˆ22µR2 +Vˆ (R,θA,θB,φ)+Vˆdd(R, SˆA, SˆB)+ Hˆ(A)as + Hˆ(B)as (4.1)43where µ is the reduced mass of the collision complex, lˆ = Jˆ− NˆA− NˆB− SˆA− SˆB isthe orbital angular momentum expressed in terms of the total angular momentumJˆ, the rotational angular momenta NˆA and NˆB, and the electronic spin angular mo-menta SˆA and SˆB of the two molecules. The operator Vˆ (R,θA,θB,φ) is parametrizedby the interaction potential between molecules A and B depending on the polar an-gles θi of each molecular axis relative to the intermolecular axis R as well as thedihedral angle φ = φA−φB.The interaction between the unpaired electrons of the molecules gives rise tothe intermolecular magnetic dipole - dipole interactionVˆdd(R, SˆA, SˆB) =−g2sµ20(24pi5)1/2 α2R3 ∑q(−)qY2,−q(Rˆ)[SˆA⊗ SˆB](2)q , (4.2)where gs is the electron g-factor, α is the fine-structure constant, µ0 is the Bohrmagneton and [SˆA⊗ SˆB](2)q is the tensor product of rank-1 tensors SˆA and SˆB [74].The asymptotic Hamiltonian Hˆ(i)as of each 2Σ molecule in a magnetic field is givenbyHˆ(i)as = BeNˆ2i + γSRNˆi · Sˆi +gsµ0BSˆZi (4.3)where Be is the rotational constant, B is the magnitude of the magnetic field, SˆZi isthe Z-component of Sˆi, and γSR is the spin-rotation interaction constant.The eigenstates of the total Hamiltonian are expanded as follows [90]|Ψ〉= 1R ∑αA,αB∑J,ΩFMαAαBJΩ(R)|NAKNA〉|SAΣA〉|NBKNB〉|SBΣB〉|JMΩ〉 (4.4)where αi denotes collectively the quantum numbers {Ni,KNi ,Si,Σi} of molecule i,KNi ,Σi and Ω are the projections of Nˆi, Sˆi and Jˆ on the BF quantization axis z andM is the projection of Jˆ on the SF quantization axis Z. In the presence of externalfields, only M is a good quantum number. This basis set is properly symmetrizedto account for the exchange symmetry of identical molecules, as described in Ref.[90].The substitution of the symmetrized expansion into the Schro¨dinger equationyields a set of coupled differential equations. For all computations discussed in44the next section, we integrate these equations on a grid of R from 4.5 A˚ to 500 A˚using the log-derivative algorithm [57, 58] with a step size of 0.05 A˚. The resultspresented in Figures 4.2 and 4.3 are obtained with the grid of R from 4.5 A˚ to100 A˚. The numerical solutions yield the log-derivative matrix that is converted tothe scattering S-matrix and the matrix of cross sections using standard equations[42]. We perform the calculation with the interaction potential operator definedas Vˆ = λVNH−NH, where VNH−NH is the 4D interaction potential surface for themaximum spin state of the NH - NH dimer recently computed by Janssen et al.[50, 51]. Each result is averaged over multiple calculations with different λ in theinterval λ = 0.5−4. This interval of λ values generates a wide range of interactionpotentials. When λ = 2.0, the interaction potential depth is similar to than of theRbCs - RbCs interaction potential [97], which is dramatically different from thatof the NH dimer. The cross sections computed with any given value of λ may beaffected by resonances and cannot be regarded as meaningful. However, multiplecalculations with different values of λ provide the reliable expectation intervals ofthe cross sections. The results of the following sections are presented with one-σstandard deviations representing ∼68% confidence intervals.In order to select the proper basis set, we computed the cross sections for elas-tic and inelastic collisions of molecules with γSR = 0.0415 cm−1 and the rota-tional constant Be = 4.2766 cm−1. These parameters correspond to the moleculeCaH(2Σ+). Figure 4.1 shows the results computed with different basis sets corre-sponding to the maximum quantum number Nmax = 2–6 of the rotational angularmomentum for each molecule in the expansion (4.4). For each calculation, the ba-sis included four total angular momentum blocks corresponding to the fixed totalangular momentum projection M = 1. Figure 4.1 illustrates that the elastic scatter-ing cross sections are almost independent of Nmax, while the cross sections for theZeeman relaxation significantly increase with Nmax, indicating that it is necessaryto include, at least 6 rotational states (N = 0–5) in the calculations.As explained above, elastic and inelastic processes at the collision energies∼ 10−3 K are particularly important for evaporative cooling. The lower panel ofFigure 4.1 shows that the inelastic cross sections at E = 10−3 cm−1 ' 1.4 mK arevery sensitive to the basis set, exhibiting a significant increase when the maximumnumber of rotational states in the basis is increased from 6 to 7. We have confirmed4510010210410610−6 10−5 10−4 10−3 10−2Collision energy (cm−1)Elastic cross section  (A°2)10−110010110210310−6 10−5 10−4 10−3 10−2Collision energy (cm−1)Inelastic cross section  (A°2 )23 45623456Figure 4.1: Cross sections for elastic (upper panel) and inelastic (lowerpanel) collisions of 2Σ molecules with γSR = 0.0415cm−1 and Be =4.2766cm−1 in a magnetic field of 100 Gauss. The colored bars cor-respond to the different number of rotational states included for eachmolecule in the basis set expansion (4.4), ranging (left to right) from 3to 7.46by the analysis of the energy dependence of the cross sections that this difference isdue to a scattering resonance. The results of Figure 4.1 are obtained with a singlepotential energy surface corresponding to λ = 1. By contrast, Figure 4.2 displaysthe elastic and inelastic scattering cross sections averaged over the results of 35calculations with different λ drawn randomly from a uniform distribution in theinterval λ = 0.5−4.Figure 4.2 demonstrates that the averaging over the data computed with differ-ent interaction potentials reduces the basis set truncation error. To illustrate thispoint better, we repeated the calculation with the basis set that included the stateswith N = 0,1 and 2 for molecule A and a variable number of rotational states formolecule B. This allowed us to examine the role of high-energy rotational stateswith the rotational angular momentum up to N = 11. The calculations with a singlepotential energy surface showed significant sensitivity to the presence of the rota-tional states with N < 12. When averaged over the results of 50 calculations withdifferent surfaces, the cross sections appear to be well converged when N ≥ 5, seeFigure 4.3. The calculations presented in the next section are performed with sixrotational states N = 0−5 for each molecule in the basis set expansion (4.4).4.3 Spin relaxation mechanismsThere are two terms in the Hamiltonian (4.1) that induce couplings between thestates with MS = +1/2 and MS = −1/2. The magnetic dipole - dipole interaction(4.2) provides direct couplings between the product states |MSA = 1/2〉|MSB = 1/2〉and |MSA = −1/2〉|MSB = −1/2〉. In addition, the states with different MS arecoupled by the spin-rotation interaction γSRNˆ · Sˆ in each molecule. When bothmolecules are in the ground rotational state N = 0, the MS = 1/2↔ MS = −1/2couplings are induced by a combination of the N = 0↔N > 0 couplings due to theanisotropy of the intermolecular interaction potential and the spin-rotation interac-tion in the states with N > 0 [42]. The cross sections for inelastic collisions leadingto spin relaxation MS = +1/2→MS = −1/2 in one or both molecules are deter-mined by the competition of these two mechanisms. While the magnetic dipole- dipole interaction is fixed, the dynamical processes induced by the spin-rotationinteraction depend on the magnitude of the spin-rotation interaction constant and4726695 26302 27641 27701 279151001021042 3 4 5 6Elastic cross section  (A°2)4.327.8 35.2462.4 509.11001011021032 3 4 5 6Inelastic cross section  (A°2)NmaxFigure 4.2: Cross sections for elastic (left panel) and inelastic (rightpanel) collisions of 2Σ molecules with γSR = 0.0415cm−1 and Be =4.2766cm−1 computed for the collision energy 10−3 cm−1 and the mag-netic field 100 Gauss. The colored bars correspond to the different num-ber of rotational states included for each molecule in the basis set expan-sion (4.4), ranging (left to right) from 3 to 7. The results are averagedover 35 calculations with different interaction potentials.481001011021031042 3 4 5 6 7 8 9 10Elastic cross section  (A°2 )1001011021031042 3 4 5 6 7 8 9 10Elastic cross section  (A°2 )1001011021032 3 4 5 6 7 8 9 10Inelastic cross section  (A°2 )1001011021032 3 4 5 6 7 8 9 10Inelastic cross section  (A°2 )Nmax NmaxFigure 4.3: Cross sections for elastic (upper panels) and inelastic (lowerpanels) collisions of molecules with γSR = 0.0415cm−1 and Be =4.2766cm−1 computed for the collision energy 10−3 cm−1 and the mag-netic field 100 Gauss. The results in the left panels are obtained witha single potential energy surfaces corresponding to λ = 1. The resultsin the right panels are averaged over the cross sections computed with50 different potentials. The calculations are performed with N ≤ 2 formolecule A in the basis set expansion (4.4). The colored bars in eachpanel correspond to the different number of rotational states includedfor molecule B in the basis set expansion (4.4), ranging (left to right)from 3 to 11.49the rotational constant of the molecule [42].In order to explore the dependence of the cross sections for elastic and inelas-tic scattering on γSR, we average the calculations over the results obtained with 48interaction potentials. The calculations are performed for molecules with the rota-tional constant Be = 4.2766 cm−1, the collision energy 10−3 cm−1 and a magneticfield magnitude of 100 Gauss. This particular magnetic field magnitude is chosenbecause the potential energy shift of 2Σ molecules in the |N = 0,MS = 1/2〉 statedue to the Zeeman effect at 100 Gauss is about 4.7×10−3 cm−1 = 6.7×10−3 K. Amagnetic trap with a field gradient up to 100 Gauss can thus confine molecular en-sembles with temperatures T < 5×10−3 K. In the limit of vanishing magnetic field,the low-energy cross sections for the Zeeman relaxation MS = 1/2→MS =−1/2tend to zero [43]. As the field increases from zero, the inelastic cross sections in theabsence of resonant scattering generally increase and the elastic cross sections areindependent of the magnetic field [93]. The results at B = 100 Gauss thus representthe highest probability of inelastic scattering in a magnetic trap with molecules atT < 5× 10−3 K. As the temperature of the molecular ensemble is reduced, thedepth of the magnetic trap can be lowered, leading to smaller cross sections forinelastic scattering and higher ratios of elastic-to-inelastic cross section [93].We note that in addition to spin relaxation, 2Σ molecules in a magnetic trapmay also undergo chemical reactions through non-adiabatic transitions to the sin-glet spin state of the collision complex [80]. These non-adiabatic transitions areinduced by the same couplings as the inelastic Zeeman relaxation. Such chemicalprocesses have been explored for the case of 3Σ molecules [98]. Unlike Zeemanrelaxation, chemical reactions produce molecular fragments with high kinetic en-ergy. Therefore, the suppression mechanism mentioned above does not apply tochemical reactions and the reaction rates are largely independent of the magneticfield. Because the chemical reactions and the Zeeman relaxation are mediated bythe same couplings, the rates of these two processes are very similar at high fields.The calculations in Ref. [98] showed that the rates for the chemical reactions andthe spin relaxation in NH - NH collisions are very similar at 100 Gauss. For heaviermolecules, the Zeeman relaxation rates must approach the reaction rates at lowermagnetic fields.Figure 4.4 shows the averaged results for the elastic and Zeeman relaxation501.3 × 1041.5 × 1041.7 × 1041.9 × 1040 0.2 0.4 0.6 0.8 1γSR (cm−1)Elastic cross section  (A°2 )●● ● ●● ● ●● ●● ●04008001200160020000 0.2 0.4 0.6 0.8 1γSR (cm−1)Inelastic cross section  (A°2 )1003005007009000 0.2 0.4 0.6 0.8 1γSR (cm−1)Ratio  (σ elσ in)! ! ! ! ! !! !! ! !!in = a + b  "SR2! ! ! ! ! ! ! ! ! ! !Figure 4.4: Cross sections for elastic (upper panel) and inelastic (middlepanel) collisions of molecules with Be = 4.2766 cm−1 computed asfunctions of γSR for the collision energy 10−3 cm−1 and the magneticfield 100 Gauss. The lower panel shows the elastic-to-inelastic crosssection ratio. The results are averaged over 48 calculations with differ-ent interaction potentials. The vertical bars indicate the 2σ interval ofthe cross section values and their ratios. The insets illustrate the poly-nomial regression fits.51cross sections with the ±σ intervals and illustrates two important observations.First, the elastic-to-inelastic ratio is very large at low magnitudes of γSR. This indi-cates that the magnetic dipole - dipole interaction is relatively ineffective. Moleculeswith small spin-rotation interaction constants are thus likely to be amenable toevaporative cooling. This is in agreement with a recent calculation of the Zee-man relaxation cross-sections in collisions of CaH(2Σ+) molecules with Li (2S)atoms [99]. In this paper, Tscherbul and coworkers showed that collision-inducedspin-relaxation may be very inefficient even in systems with very large interac-tion anisotropy, providing it is mediated by weak spin-dependent interactions. Thespin-dependent interactions serve as a bottleneck suppressing the Zeeman relax-ation even in the limit of infinitely strong interaction anisotropy [99].Second, Figure 4.4 shows that the inelastic cross sections can be well approxi-mated by a quadratic function of γSR. This indicates that the spin-relaxation processis largely driven by second-order interactions involving the matrix elements of thespin-rotation interaction [42]. The role of the magnetic dipole-dipole interactionis apparent at low values of γSR, where the dependence of the Zeeman relaxationcross sections on γSR is weak. The strong dependence of the inelastic cross sectionson γSR at γSR > 0.4 cm−1 suggests that the mechanism of the Zeeman relaxation atthese strengths of the spin-rotation interactions is dominated by the γSR-dependentcouplings. The elastic-to-inelastic ratio is inversely proportional to γSR and ap-pears to be above 100, even for molecules with the largest possible spin-rotationalconstants ∼ 1 cm−1. This value of the elastic-to-inelastic ratio is considered to becritical for the efficacy of evaporative cooling and is often used as a limit belowwhich the cooling experiments become unfeasible [100].In the case of collisions with atoms, the Zeeman relaxation of molecules inthe state |N = 0,MS = 1/2〉 is very sensitive to the energy separation between therotational states [25, 101]. In order to examine the sensitivity of the Zeeman re-laxation in molecule - molecule collisions to the rotational energy splittings, wepresent in Figures 4.5 and 4.6 the elastic and inelastic cross sections as functionsof the rotational constant. The calculation in Figure 4.5 is for molecules with a verysmall value of γSR = 0.0415 cm−1 and the results in Figure 4.6 are for moleculeswith γSR = 0.5 cm−1. When the spin-rotation interaction is weak, the elastic andinelastic cross sections exhibit a weak dependence on the rotational constant of the521 × 1041.4 × 1041.8 × 1042.2 × 1041 2 3 4 5 6 7 8 9 10Be  (cm−1)Elastic cross section  (A°2 )●● ● ● ●●●● ●●0501001502001 2 3 4 5 6 7 8 9 10Be  (cm−1)Inelastic cross section  (A°2 )300600900120015001 2 3 4 5 6 7 8 9 10Be  (cm−1)Ratio  (σ elσ in)Figure 4.5: Cross sections for elastic (upper panel) and inelastic (middlepanel) collisions of molecules with γSR = 0.0415 cm−1 computed asa function of Be for the collision energy 10−3 cm−1 and the magneticfield 100 Gauss. The lower panel shows the elastic-to-inelastic crosssection ratio. The results are averaged over 20 calculations with differ-ent interaction potentials. The vertical bars indicate the 2σ interval ofthe cross section values and their ratios.531 × 1041.2 × 1041.4 × 1041.6 × 1041.8 × 1041 2 3 4 5 6 7 8 9 10Be (cm−1)Elastic cross section  (A°2 )●●●● ● ● ● ● ● ●0500100015002000250030001 2 3 4 5 6 7 8 9 10Be (cm−1)Inelastic cross section  (A°2 )100300600900120015001 2 3 4 5 6 7 8 9 10Be (cm−1)Ratio  (σ elσ in)! ! !! ! ! ! ! ! !! ! ! !! ! !!!!Figure 4.6: Cross sections for elastic (upper panel) and inelastic (middlepanel) collisions of molecules with γSR = 0.5 cm−1 computed as a func-tion of Be for the collision energy 10−3 cm−1 and the magnetic field 100Gauss. The lower panel shows the elastic-to-inelastic cross section ratio.The results are averaged over 35 calculations with different interactionpotentials. The vertical bars indicate the 2σ interval of the cross sectionvalues and their ratios. The insets illustrate the polynomial regressionfits.54molecule. This indicates that, at this magnitude of the spin-rotation interaction,the Zeeman relaxation is largely driven by the magnetic dipole - dipole interac-tion. By contrast, the inelastic cross sections displayed in Figure 4.6 decreasemonotonically as the rotational constant Be increases. This is another evidence thatthe mechanism of the Zeeman relaxation induced by the spin-rotation interactionclearly dominates at these values of γSR. The results presented in Figure 4.6 showthat the elastic-to-inelastic ratios may exceed 100 even for molecules with the spin-rotation interaction constant as large as γSR = 0.5 cm−1, providing the rotationalconstant is greater than 4 cm−1. The range of 2Σ molecules expected to be colli-sionally unstable in magnetic traps at T < 10−3 K is thus restricted to moleculeswith γSR > 0.5 cm−1 and Be < 4 cm−1.4.4 ConclusionDiatomic molecules in the 2Σ electronic state represent the simplest class of molec-ular radicals that can be confined in a magnetic trap. While the magnetic trappingof 2Σ molecules has been achieved, the feasibility of evaporative cooling of 2Σmolecules to ultracold temperatures depends on the magnitudes of the cross sec-tions for elastic scattering and inelastic Zeeman relaxation and is yet to be experi-mentally demonstrated. The present work provides insight into the relative magni-tudes of these cross sections. In order to draw general conclusions, we performedthe computations with multiple potential energy surfaces, yielding the expectationintervals for both the elastic and inelastic scattering cross sections.The first important result of this work is illustrated by Figures 4.1, 4.2 and 4.3.These calculations show that the cross sections averaged over multiple potential en-ergy surfaces are much less sensitive to the basis set truncation error. This indicatesthat statistically meaningful results can be obtained by averaging over multiple cal-culations with fairly restricted basis sets. This also made the calculations presentedin Figures 4.4, 4.5 and 4.6 feasible. Reducing the number of rotational states in thebasis set from 7 to 6 reduces the computation time of each data point on a fast CPUfrom about 1 month to about one week. It thus took about 48 weeks of CPU timeto compute each interval depicted in Figure 4.4. Averaging over multiple interac-tion potentials also allows for the analysis of the expectation values of the cross55sections as functions of the molecular parameters.We analyzed the results as functions of the molecular parameters (the spin-rotation interaction strength and the rotational constant magnitude) for the col-lision energy 10−3 cm−1 and the magnetic field 100 Gauss. This collision en-ergy corresponds approximately to the energy of molecules at about 1.4×10−3 K.This magnetic field magnitude is representative of the field strength experiencedby molecules in magnetic traps with the trap depth ∼ 5× 10−3 K. The inelasticcross sections are expected to be smaller and the elastic-to-inelastic ratios higherat lower magnetic fields.Our results show that the elastic-to-inelastic ratios are consistently greater than100, except for molecules with very large constants of the spin-rotation interac-tion and small rotational constants (Be < 4) cm−1. This indicates a good prospectfor evaporative cooling of a great variety of 2Σ molecules to temperatures be-low 10−3 K. Molecules with small constants of the spin-rotation interaction ap-pear to be the best candidates for evaporative cooling, irrespective of the magni-tude of the rotational constant. This challenges the previous conclusion [42] thatonly molecules with large rotational constants should be amenable to collisionalcooling. Molecules with the spin-rotation interaction constants > 0.4 cm−1 havefavourable collision properties only if their rotational constant Be > 5 cm−1. Refs.[100] and [102] provide a list of the spectroscopic constants for a variety of di-atomic molecules in the 2Σ electronic state. Figure 4.7 displays the lower limit ofthe expectation values of the elastic-to-inelastic ratios for these molecules obtainedby the interpolation of our results. To obtain these data, we used the lowest pointof a 2σ deviation of the elastic-to-inelastic ratio from the mean value, correspond-ing to the ∼ 95 % confidence interval. The figure illustrates that the majority ofthe selected 2Σ radicals should be amenable to evaporative cooling starting at mil-liKelvin temperatures. It is important to emphasize that the results presented inFigure 4.7 should be interpreted as expectation values accurate to within ∼95 %.In any specific case, the collisional interaction between the molecules may be af-fected by details of the interaction potentials, leading to lower elastic-to-inelasticcross section ratios. Nevertheless, the results of Figure 4.7 are very encouragingfor the prospect of evaporative cooling of 2Σ molecules as most of the data appearto be well above 100.56●●● ● ●●●● ●●●●●●●● ● ●PdHRhCPdD HgD HgHBSBaFBaDSiNCPCdHBaHZnDCaFZnHCN BOCaH0100200300400500Ratio  (σ elσ in)Figure 4.7: The lowest magnitude of the expected elastic-to-inelastic crosssection ratio for selected 2Σ molecules. The results are obtained by in-terpolation of the data in Figures 4.4, 4.5 and 4.6 using the spectroscopicconstants borrowed from Ref. [100].57Chapter 5Collision lifetimes of polyatomicmolecules at low temperatures:benzene - benzene vs benzene -rare gas atom collisionsIn this chapter, we use classical trajectory calculations to study the effects of theinteraction strength and the geometry of rigid polyatomic molecules on the for-mation of long-lived collision complexes at low collision energies. We illustratesthat the mean lifetimes of the collision complexes increase monotonically withthe strength of the atom – molecule interaction, while the mean lifetimes of poly-atomic molecule-molecule collision complexes are significantly reduced due tonon-ergodic effects.5.1 IntroductionThe rapid development of experimental techniques for the production of cold (< 5K) and ultracold (< 10−3 K) molecular gases has opened new opportunities forfundamental research in molecular physics[70, 103], chemical dynamics [80, 104,105] and beyond [13, 106, 107]. It is now possible to prepare large ensemblesof alkali metal dimers at microKelvin temperatures [108–111] and a variety of di-58atomic molecules at temperatures near 1 Kelvin [22, 112, 113]. Fuelled by thepromise of new applications [45, 46], there is currently growing interest in extend-ing the cooling techniques to complex polyatomic molecules. At present, thereare two viable routes to cooling polyatomic molecules: deceleration of molecularbeams [114] and buffer-gas cooling [30]. The former relies on the Stark or Zeemaneffect to slow down molecules by gradients of time-dependent fields. It is partic-ularly well suited for molecules with large electric or magnetic dipole moments.The buffer-gas cooling technique relies on the thermalization of molecules in acold buffer gas [22]. The thermalization process is driven by momentum transferin elastic collisions of molecules with the buffer-gas atoms. While the buffer-gascooling technique appears to be general and applicable to any molecular system,the collisional thermalization at low temperatures may be plagued by various un-wanted processes, such as resonant scattering and the ensuing formation of atom -molecule clusters.The probability of cluster formation is clearly related to the complexity ofmolecules. The interaction of a slowly moving atom with an infinitely extendedsurface deposits the kinetic energy into phonons, trapping the atom on the surface[115]. On the other hand, the scattering of an atom by a diatomic molecule occurson a fast time scale. Falling between these two limits, the collision dynamics ofcomplex polyatomic molecules with low temperature atoms is characterized by amanifold of scattering resonances [116]. These resonances enhance the lifetimeof the collision complex, allowing three-body recombination and the formationof clusters [117–119]. The density of resonances, and consequently the collisionlifetime and the probability of clustering, are expected to increase as the collisionpartners become more complex [30]. The complexity of the collision partners canbe quantified by two parameters: the polarizability determining the inter-molecularor atom - molecule interaction strength and the density of accessible states. It isgenerally expected that the collision lifetimes must increase with the polarizabilityand the density of internal states of the collision partners. For example, Patter-son and coworkers [30] recently reported the results of the cooling dynamics fornaphthalene (C10H8) in the He buffer gas. They observed no clustering of the Heatoms on molecules and attributed the observed loss of molecules to the formationof naphthalene - naphthalene clusters, assuming that collisions of molecules with59molecules are more complex.In the present work, we examine the role of the polarizability of the collisionpartners and the complexity of the molecules on the collision lifetimes at tem-peratures between 5 and 10 Kelvin that correspond to the buffer-gas cooling ex-periments of Patterson and coworkers [30]. In order to do this, we carry out theclassical trajectory calculations for collisions of benzene (C6H6) molecules withrare gas atoms He – Xe and for collisions of benzene molecules with benzenemolecules. Surprisingly, we find that, despite the much stronger interaction en-ergy, the lifetimes for benzene - benzene collisions are comparable to those forbenzene - Ne collisions and are significantly shorter than the lifetimes of benzene- Xe collisions. Our analysis shows that the rigid geometry of benzene reduces thelifetimes of benzene - benzene collisions by constraining the molecules to a partof the configuration space. This finding indicates that the loss of molecules due tomolecule - molecule collisions in buffer-gas cooling experiments with rigid poly-atomic molecules may be less significant than previously thought. This result isalso important for the prospect of evaporative cooling of complex molecules, a pro-cess in which the Maxwell-Boltzmann distribution of a molecular gas is disturbedby removing the high-energy tail and the molecules return to thermal equilibriumby means of molecule - molecule collisions. This work complements a recent studyby Li et al. [47, 48] elucidating the effects of the molecular size and geometry onthe formation of long-lived collision complexes with He atoms.5.2 MethodologyWe employ the numerical method for the classical trajectory calculations describedby Li and Heller [47]. This approach uses the formulation of the classical equationsof motion in terms of the quaternions and is particularly well suited for collisionsof rigid polyatomic molecules at low temperatures. We consider binary collisionsof benzene molecules with rare gas (Rg) atoms He, Ne, Ar, Kr and Xe as well asbenzene - benzene collisions. In all calculations, the benzene molecule is treatedas a rigid rotor.605.2.1 Molecular dynamics simulation of rigid molecule-atomcollisionsThe collisional dynamics of a rigid molecule and a structureless atom like C6H6-Rgsystem can be divided into the relative translational motion between the collisioncomplex and the rotational motion of the rigid molecule[120]. The former is de-scribed in a space-fixed frame, where the origin O coincides with the center of massof the collision complex, whereas the latter is described using a body-fixed framewith axes corresponding to the molecular principle axes and the origin O′ locatedat the center of mass of the molecule. The relative translational motion depends onthe interaction force between collision complex FSF, which is given by the gradientof the interaction potential V (R),FSF =−∇V (R) (5.1)The construction of the interaction potential surface for the collision complexeswill be described in detail in Section 5.2.3. The relative translational motion at anytime t is specified by the position vector of atom relative to the center of mass ofmolecule R and linear momentum P obtained from the coupled ordinary differen-tial equations (ODE)dRdt= pµ (5.2)dPdt= FSF (5.3)where µ is the reduced mass of the collision complex.The rotational motion of a rigid molecule is governed by the net torque actingon its center of mass. The net torque in the space-fixed frame is given byNSF = R× (−FSF) (5.4)A full description of the rotational motion needs to specify the orientation ofthe molecule and the angular velocity ω . The components of the angular velocityωx,y,z can be expressed in terms of Euler angles as61ω =ωxωyωz =cosψ sinθ sinψ 0−sinψ sinθ cosψ 00 cosθ 1θ˙φ˙ψ˙= Ξ α˙ (5.5)where θ˙ , φ˙ and ψ˙ are the first time derivatives of the Euler angles.Since the simulation of the rotational motion at each time interval involves thenumerical integration of θ˙ , φ˙ and ψ˙ using the components of the angular velocityωx,y,z at the beginning of the time interval. Eq(5.5) and its inverse α˙ = Ξ−1ω willbe frequently used in the process of integration. When θ = 0 or pi , two of therotation axes will coincide with each other, and the matrix Ξ becomes singular,leading to numerical instability. To avoid the problem of singular matrices, weuse an alternative representation of the orientation and transform Euler angles intoquaternions parameters byq0 = cos(θ/2)cos((φ +ψ)/2) (5.6)q1 = sin(θ/2)cos((φ −ψ)/2) (5.7)q2 = sin(θ/2)sin((φ −ψ)/2) (5.8)q3 = cos(θ/2)sin((φ +ψ)/2) (5.9)The quaternions satisfy3∑m=0q2m = 1 (5.10)Correspondingly, the components of the angular velocity ωx,y,z can be replaced bylinear combinations of first time derivatives of quaternions q˙0,1,2,3 from the relationωxωyωz0= 2Wq˙0q˙1q˙2q˙3(5.11)62whereW =−q1 q0 q3 −q2−q2 −q3 q0 q1−q3 q2 −q1 q0q0 q1 q2 q3(5.12)is an orthogonal matrix.Based on the Euler’s rotation equation, the first time derivative of each compo-nent of the angular velocity ω˙x,y,z has the general form asω˙x =1Ix[NBFx +(Iy− Iz) ·ωx ·ωz] (5.13)ω˙y =1Iy[NBFy +(Iz− Ix) ·ωz ·ωx] (5.14)ω˙z =1Iz[NBFz +(Ix− Iy) ·ωx ·ωy] (5.15)where Ix,y,z are the principal moments of inertia of the molecule and NBFx,y,z are com-ponents of the torque in the body-fixed frame. The transformation of the torqueNSFx,y,z from a space-fixed frame to a body-fixed frame is achieved byNBF = RNSF (5.16)where R is the rotation matrix defined in terms of the quaternions as follows:R =12 −q22−q23 q1q2 +q0q3 q1q3−q0q2q1q2−q0q3 12 −q21−q23 q2q3 +q0q1q1q3 +q0q2 q2q3−q0q1 12 −q21−q22 (5.17)Since we employ the quaternions q0,1,2,3 and their first time derivatives q˙0,1,2,3in the description of the rotational motion of the molecule, and the change of therotational state over time is determined by the inverse of Eq (5.11) and a set of63second order differential equations:q¨0q¨1q¨2q¨3= 12−q1 −q2 −q3 q0q0 −q3 q2 q1q3 q0 −q1 q2−q2 q1 q0 q3ω˙xω˙yω˙y−2∑3m=0 q˙2m(5.18)where q¨0,1,2,3 are the second time derivatives of q0,1,2,3.Since ω˙x,y,z can be written in terms of quaternions and their first time deriva-tives, the Euler angles and the components of angular velocity will not appear inthe simulation.Finally, we denote the combined dynamic state of the collision complex in thetwo coordinates at time t by a vector consisting of 14 variables Y(t), namely, thecomponents of relative distance RX ,Y,Z(t) and linear momenta PX ,Y,Z(t) along theaxes in the space-fixed frame, quaternions q0,1,2,3(t) and their first time derivativesq˙0,1,2,3(t),Y(t) =(RX(t),RY (t),RZ(t),PX(t),PY (t),PZ(t),q0(t),q1(t),q2(t),q3(t),q˙0(t), q˙1(t), q˙2(t), q˙3(t))>(5.19)The time derivative of the dynamic state Y˙(t) ≡ dY(t)/dt determines the rateof change of the dynamic state over time, and given byY˙(t) =(R˙X(t), R˙Y (t), R˙Z(t), P˙X(t), P˙Y (t), P˙Z(t), q˙0(t), q˙1(t), q˙2(t), q˙3(t),q¨0(t), q¨1(t), q¨2(t), q¨3(t))>=(PX(t)µ ,PY (t)µ ,PZ(t)µ ,FSFX ,FSFY ,FSFZ , q˙0(t), q˙1(t), q˙2(t), q˙3(t),q¨0(t), q¨1(t), q¨2(t), q¨3(t))>(5.20)where FSFX ,Y,Z are the projections of the interaction force along the coordinate axesof the space-fixed frame.64During the simulation, we use the fourth-order Runge-Kutta method to solvethe 14 ODEs of motion defining Y˙(t), and the numerical integration at each timeinterval ∆t is based on the dynamic state at the beginning of the time interval Y(t).5.2.2 Initialization of the dynamic state of a collision complexThe initial configuration of Y(t = 0) is chosen as follows:We first initialize the relative position R0 and linear momentum p0 in a space-fixed frame. Assuming the virtual particle with collision energy Ecol and mass µequal to the reduced mass of the the collision complex is directed along a space-fixed Z-axis, the relative distance R0 is set to be large enough so that the interactionpotential is negligible. The impact parameter bl is chosen randomly from the inter-val [0,bmaxl ]. The lifetime of the quasi-bound collision complex is always zero forbl > bmaxl .RX(0) = 0;RY (0) = bl;RZ(0) = −√R20−b2l ;pX(0) = 0;pY (0) = 0;pZ(0) =√2µEcol;The dynamic state of the molecule in the body-fixed frame is initialized interms of Euler angles and the angular velocity, which are then transformed toquaternions and their first time derivatives according to Eqs (5.6–5.9) and the in-verse of Eq (5.11). The Euler angles θ ,φ ,ψ and the magnitudes of ωx,y,z are sam-pled randomly with the constraint Ixω2x + Iyω2y + Iyω2y = 2Erot. The initial rotationalenergy Erot is chosen randomly according to the Boltzmann distribution and has acutoff at 10 cm−1Pr(Erol) ∝ exp(−ErolkTrol) (5.21)65where k is the Boltzmann constant, and Trot is the rotational temperature.5.2.3 Potential energy surfaces for C6H6-Rg and C6H6-C6H6The potential energy surfaces for the interaction of benzene with rare gas atomswere constructed using the semi-empirical approach of Pirani and coworkers [121].This method treats a polyatomic molecule like benzene as an ensemble of diatoms,coinciding with the bonds, and represents the potential surfaces as the sums ofpairwise Rg - CH bond and Rg - CC aromatic bond interaction energies, optimizedbased on accurate ab initio calculations and measurements of bond polarizabilities.Each Rg-CH and Rg-CC aromatic bond interaction is a simple analytical functionof the pair radial and angular coordinates as well as the corresponding eight poten-tial parameters, and the dependencies are expressed asVab(r,α) = ε(α)[(32+2x2)·(1x)10+4x2−(5+2x22+2x2)·(1x)6](5.22)whereε(α) = ε⊥ sin2(α)+ ε‖ cos2(α) (5.23)rm(α) = rm⊥ sin2(α)+ rm‖ cos2(α) (5.24)r is the distance between the atom and the center of the bond, and x is the reduceddistance x = rrm(α) , where rm(α) is the position of the potential well and ε(α) isits depth. The parameters ε⊥,ε‖,rm⊥ and rm‖ represent the well depth and locationfor the parallel and perpendicular approaches of the probe atom to the diatom. α isthe angle between the bond axis and the vector connecting the atom to the centerof the bond axis. The above eight potential parameters for each Rg-CH and Rg-CCaromatic bond are given in Table 5.1. Table 5.2 lists the value of the equilibriumdistance Re and the well depth De for some limiting geometries of the Rg-C6H6collision complexes calculated by this method.The geometries represented in Table 5.2 are illustrated in Figure 5.1 (b).Theradial dependence of the potential energy for three limiting configurations illus-trated in Figure 5.1 (b) is shown in Figure 5.2. In general, the interaction strength66System Rm⊥ Rm‖ ε⊥ ε‖CCar-He 3.583 4.005 0.860 0.881CH-He 3.234 3.480 1.364 1.016CCar-Ne 3.629 4.015 1.706 1.862CH-Ne 3.316 3.549 2.544 1.960CCar-Ar 3.879 4.189 3.895 4.910CH-Ar 3.641 3.851 4.814 3.981CCar-Kr 3.997 4.284 4.824 6.365CH-Kr 3.782 3.987 5.667 4.781CCar-Xe 4.163 4.425 5.654 7.849CH-Xe 3.976 4.175 6.233 5.384Table 5.1: Potential parameters for the Rg-CH and Rg-CC aromatic bondpair. The perpendicular and parallel component of Rm and ε are in A˚and in meV, respectivelySystemout-of-plane vertex-in plane side-in-planeRe De Re De Re DeC6H6-He 3.25 -78.95 4.80 -42.59 5.23 -25.52C6H6-Ne 3.30 -157.2 4.87 -82.69 5.29 -50.37C6H6-Ar 3.58 -355.0 5.16 -179.0 5.53 -117.1C6H6-Kr 3.70 -439.5 5.28 -221.0 5.65 -148.2C6H6-Xe 3.88 -513.3 5.46 -258.5 5.80 -179.4Table 5.2: The equilibrium distance Re in A˚ and well depth De in cm−1 forthree relevant cross sections of the C6H6–Rg potential energy surface.The out-of-plane, vertex-in-plane and side-in-plane geometries are illus-trated in Figure 5.1.increases from He to Xe. The global minimum of the potential energy occurs in theout-of-plane configuration with the Rg atom placed on the 6-fold symmetry axis ofC6H6, for all Rg atoms.For the benzene dimer, we borrow the interaction potential from Ref. [122].The method proposed in this work represents each molecule as a combination oftwo concentric rings of atoms, an inner carbon ring and an outer hydrogen ring (seeFigure 5.3 (a)), leading to a computationally efficient model of the interaction po-tential surface as the sum of four ring-ring interactions: a carbon ring-carbon ringinteraction (EC−C), a carbon ring-hydrogen ring interaction (EC−H), a hydrogen67xyzRθφ(a)out-of-planevertex-in-plane side-in-plane(b)θ = 0◦θ = 90◦φ = 0◦θ = 90◦φ = 30◦Figure 5.1: (a) The internal coordinates used to describe the Rg – C6H6 colli-sion complex: R is the relative distance between centers of mass, θ (thepolar angle) is the angle between R and the 6-fold symmetry axis of ben-zene, and φ (the azimuthal angle) is the angle between the projection ofR onto the molecular plane and the x axis, which is perpendicular to aC-C bond of benzene. (b) The relative geometries of the Rg–C6H6 col-lision complex for the out-of-plane, vertex-in-plane, and side-in-planeconfigurations.68−600−400−2000Interaction Potential(cm−1 ) Out−of−PlaneHeNeArKrXe C6H6−600−400−2000Interaction Potential (cm−1 ) Vertex−in−Plane3 4 5 6 7 8 9 10−600−400−2000R  (A° )Interaction Potential (cm−1 ) Side−in−PlaneFigure 5.2: Radial dependence of the interaction energy of benzene with He(black solid line), Ne (red dashed line), Ar (blue dotted line), Kr (greendot-dashed line), and Xe (purple long-dashed line) compared with thecorresponding benzene-benzene interaction averaged over the normalangles of the molecules (magenta dot-dashed line). Upper panel: out-of-plane configuration; central panel: vertex-in-plane configuration; lowerpanel: side-in-plane configuration. The geometries are illustrated inFigure 5.1.69ring-carbon ring interaction (EH−C) and a hydrogen ring-hydrogen ring interaction(EH−H). Thus, the total interaction potential energy isEtotal = EC1−C2 +EC1−H2 +EC2−H1 +EH1−H2 (5.25)In a spherical polar coordinate system, we fix the first benzene ring on the x-yplane, and make its center coincide with the origin O of the coordinate system,so the geometry of the benzene dimer is determined by the relative position of thesecond ring to the first ring, which can be described by the vector ~R≡ (R,θ ,φ) con-necting the center of mass of the two rings, and the angle γ between the principalz-axes of the two rings.(see Figure 5.3 (b)).The interaction energy between a given pair of rings is calculated byEr1−r2 = r1r2η1η2(−AJ3 +BJ6) (5.26)where A and B are constants (Table 5.3); r1 and r2 are radii of the first ring andsecond rings, respectively; η1 and η2 are the atomic line densities of the first ringand second rings, respectively. (Table 5.4).A B(kcal/mol×A˚6) (kcal/mol×A˚12)Carbon-Carbon 560.44 1121755.66Carbon-Hydrogen 129.67 91727.95Hydrogen-Hydrogen 24.10 4841.34Table 5.3: Attractive and repulsive constants.Radius r(A˚) Atomic density η (A˚−1)Carbon ring 1.40 0.682Hydrogen ring 2.48 0.385Table 5.4: Radii and atomic densities of carbon ring and hydrogen ring.70RRR1R2Sandwich (S)T-shaped (T) Parallel-Displaced (PD)x(c)yzRθφhydrogen ringcarbon ring(a) (b)γFigure 5.3: (a) Benzene as a combination of two rings: inner carbon ringand outer hydrogen ring. (b) Model for the interaction between the tworings: R is the relative distance between the centers of mass, θ (the polarangle) is the angle between R and the 6-fold symmetry axis of benzene,φ (the azimuthal angle) is the angle between the projection of R ontothe molecular plane and the x axis, which is perpendicular to a certainC-C bond of benzene, and γ (the normal angle) is the angle between the6-fold symmetry axes of the two molecules. (c) The relative geometriesof the benzene dimer for the sandwich, T-shaped and parallel-displaced(PD) configurations.71Jn (n = 3,6) can be evaluated analytically as followsJn = 4pi∞∑m=0(4r21)m(n2)m(n+12 )mm!×m∑i=02i∑j=02m−2i∑k=0(2i)!(2m−2i)!h2m− j−kr j+k2 cos2i− j φ sin2m−2i−k φ cosk γi!(m− i)!(2i− j)!(2m−2i− k)!×j∑l=0k∑p=0k+l−p∑q=0(−1) j−l+q2 j−l+p+q(k+ l− p)!cosl+pψ0 sin j+k−l−pψ0l!p!q!( j− l)!(k− p)!(k+ l− p−q)!× Γ(β )Γ(γ−β )Γ(γ)F(α,β ;γ;µ)(r21 + r22 +h2 +Z2 +K5)n+2m(5.27)whereZ = Rcosθ (5.28)h = Rsinθ (5.29)K5 = 2r2[(hcosφ)2 +(hcosγ sinφ +Z sinγ)2]1/2 (5.30)sinψ0 = 2r2(hcosγ sinφ +Z sinγ)/K5cosψ0 = 2r2hcosφ/K5α = n+2mβ = ( j− l + p+2q+1)/2γ = j− l + p+q+1µ = 2K5r21 + r22 +h2 +Z2 +K5and Γ(z) is the gamma function and (a)n = a(a+1) · · ·(a+n−1) is the Pochham-mer symbol, and F(α,β ;γ;µ) = ∑∞i=0 (α)i(β )i(γ)ii! µi is the hypergeometric function.This model of the benzene - benzene interaction potential reproduces the mainfeatures of the potential energy surface, as shown by Tables 5.5 and 5.6 summariz-ing the dimer geometries and the well depths at the minima of the interaction en-ergies. The optimized equilibrium structure of the benzene dimer roughly matchesthe theoretical results calculated by various quantum chemistry methods except forthe parallel-displaced configuration due to effects of the electrostatic interactionand induction in a pi−pi system.72Method S TPDR1 R2MP2 [123] 3.8 5.0 3.4 1.6CCSD(T) [124] 4.0 5.1 3.6 1.8DF-LMP2 [125] 3.72 4.89 3.37 1.57DF-SCS-LMP2 [125] 3.89 5.04 3.55 1.63SAPT [126] 3.816 4.97 3.48 1.84Current PES 3.58 5.18 3.53 1.0Table 5.5: The optimized equilibrium structure (A˚) of the benzene dimer.The sandwich, T-shaped and parallel-displaced (PD) geometries are il-lustrated in Figure 5.3.Method S TPDR1 R2MP2 [123] -989.8 -1049.3 -1441.0CCSD(T) [124] -465.2 -783.5 -776.5DF-LMP2 [125] -1114.7 -1139.2 -1548.4DF-SCS-LMP2 [125] -622.9 -781.4 -910.8SAPT [126] -619.8 -948.5 -939.4Current PES -1186.7 -583.4 -1162.2Table 5.6: The interaction energies (cm−1) at the minima of the potential en-ergy surface for the benzene dimer. The sandwich, T-shaped and parallel-displaced (PD) geometries are illustrated in Figure 5.3.735.2.4 Details of numerical calculationsThe classical equations of motion were propagated using the fourth-order Runge-Kutta algorithm [47]. For the Rg – benzene collisions, we used the fixed integra-tion step size: 10−6 ns for He, 10−7 ns for Ne and Ar, and 10−8 ns for Kr and Xe.For benzene-benzene collisions, we decreased the stepsize for the computation ofeach trajectory until the rigidity of the molecules and the total energy remainedunperturbed throughout the trajectory. Longer and more complicated trajectoriesrequired smaller integration stepsize. The trajectories are terminated once the sepa-ration between the colliding particles reaches 50 A˚. The collision lifetime is definedas the time between the first and the last turning points in the relative translationalenergy. The initial conditions were chosen to cover the entire phase space and allpossible orientations of the benzene molecules. The trajectory computations werecarried out for fixed, randomly sampled, impact parameters b, with the largest im-pact parameter bmax determined so that the collision lifetimes from the calculationswith b > bmax are all zero (i.e. the trajectories with b > bmax have no or only oneturning point). The results presented are averages over 105 trajectories for eachcomplex at a given collision energy, yielding convergence of the mean lifetime tobetter than 0.1 %. The weighted mean lifetime at each collision energy is aver-aged over the Boltzmann distribution of the initial rotational energy of moleculesand the collision energy dependence of the collision lifetimes thus obtained is in-tegrated with the Maxwell-Boltzmann distribution, to yield the thermally averagedobservables at fixed temperature.5.3 ResultsWhen immersed in a cold buffer gas of noble gas atoms, polyatomic molecules as-sume the temperature of the buffer gas by loosing or gaining energy via momentum-transfer collisions. Some collisions may lead to scattering resonances, resulting inthe formation of long-lived collision complexes. If the lifetime of the collisioncomplexes is longer than the mean time between collisions (determined by the gasdensity), they may encounter another atom and form bound states by the processof three-body recombination. This leads to clustering and impedes the thermaliza-tion. The average lifetime of the collision complexes is thus a critical parameter74that determines the feasibility of a buffer-gas cooling experiment. The mean timebetween collisions can be estimated as τ = 1/N , whereN = nσ√8kBTpiµ (5.31)is the number of collisions per second. In this equation, n is the density of thebuffer gas, σ is the collision cross section and µ is the reduced mass of the collisioncomplex. Assuming the collision cross section σ = 10−13 cm2 [30], the reducedmass µ equal to the mass of He and a temperature of 6.5 Kelvin, the time betweencollisions is about 1300 ps for the density n = 4× 1017 atoms/cm3. This densitywas used in the buffer gas cooling experiments of Patterson et al [30]. Note that thedensity in buffer-gas cooling experiments is typically lower than the density usedin this particular experiment so the value 1300 ps can be used as an estimate of thelower threshold for clustering.Figures 5.4 and 5.5 illustrate that the lifetimes of Rg - benzene collisions in-crease monotonically down the column of the Rg atoms from He to Xe. Figure 5.4shows the average lifetime of the Rg - benzene collision complexes as a functionof the collision energy. The data points at each collision energy are averaged overthe Boltzmann distribution of the rotational states of benzene corresponding to thetemperature 6 Kelvin. For each collision system, the average lifetime is monoton-ically decreasing with the collision energy. The inset of Figure 5.4 compares thefraction of short trajectories (τ < 10 ps) and the fraction of long trajectories (τ >100 ps) for each Rg at a low collision energy of 2 cm−1. The fraction of long tra-jectories increases from 3 % for He to 43 % for Xe, while the contribution of shorttrajectories decreases from 78 % for He to 47 % for Xe.Figure 5.5 shows the collision lifetimes as functions of temperature. The dataof Figure 5.5 were obtained by first computing the energy dependence, similar tothe one displayed in Figure 5.4 but averaged over the Boltzmann distribution withthe corresponding temperatures, and then integrating the energy dependence withthe Maxwell-Boltzmann distribution. When the temperature is reduced from 10 to5 Kelvin, the fully averaged collision lifetimes increase by a factor of ≈ 1.8 for Heand ≈ 1.5 for Xe.75l l l l l l l l l l l l l l l l l l l l l l l l l0 5 10 15 20 250100200300400Collision Energy (cm−1) Lifetime  (ps)l HeNeArKrXe< 10 ps > 100 ps0. 5.4: The mean lifetimes of the quasi-bound collision complexes ofbenzene with He (black circles), Ne (red squares), Ar (blue diamonds),Kr (green up-triangles), Xe (purple down-triangles) as functions of thecollision energy. The inset shows the distributions of short trajectoriesand long trajectories for each rare gas at a collision energy of 2 cm−1.For each collision energy, the results are averaged over the rotationalenergy of benzene with the Boltzmann distribution corresponding to thetemperature 6 K.76l l l l l l l l l l l5 6 7 8 9 10050100150200250Temperature (K)Lifetime (ps)l HeNeArKrXeFigure 5.5: The mean lifetimes of the quasi-bound collision complexes ofrare gas atoms with benzene as functions of temperature.77There are two factors that may contribute to the change of the collision life-times from He to Xe: the increase of the reduced mass and the increase of theRg atom polarizability, and hence of the Rg – benzene interaction strength. Thelarger binding energy is expected to result in more energy conversion between thetranslational and internal degrees of freedom of the molecule, thereby enhancingthe collision lifetimes. The effect of the mass is less clear. In order to disentanglethe effects of the binding energy and the mass change, we compute the collisionlifetimes for an artificial system with the reduced mass and all other collision pa-rameters set to be identical to those in the Xe - benzene calculations, but with theinteraction potential set to be the same as for the He - benzene system. The com-parison of the results displayed in the left panel of Figure 5.6 with those for Xe -benzene collisions and He - benzene collisions shows that (i) the increase of themass and the polarizability both enhance the collision lifetimes; (ii) the increaseof the polarizability has a much more significant effect on the collision lifetimes,particularly at collision energies > 10 cm−1.In order to obtain a more detailed information on the relative effects of massand polarizability, we computed the collision lifetimes for a series of collision sys-tems with varying reduced mass and interaction potentials. Figure 5.7 presents theresults for a collision energy of 5 cm−1. Each curve shown in Figure 5.7 is com-puted with a fixed interaction potential. Each data point corresponds to a differentreduced mass. These calculations confirm that the increase of mass and polariz-ability both enhance the collision lifetimes. They also illustrate that the collisionlifetimes are more sensitive to the reduced mass when the polarizability is larger.The data of Figure 5.7 can be well approximated by the following expression:τc = a|De|+bµ+ c|De|µ, (5.32)where τc is the mean collision lifetime in ps, a = 0.0617 ps/cm−1, b = 0.658 psmol/g and c= 0.00573 ps mol / g cm−1, De is the energy of the interaction potentialfor the Rg - benzene systems at the global minimum (column 1 of Table I) and µis the reduced mass of the collision system. The fit (5.32) is a multiple linearregression with the adjusted R2 coefficient R2 = 0.97, illustrating that the meancollision lifetime is a nearly linear function of both the interaction strength and the78reduced mass.As can be seen in Tables 5.1, 5.2, 5.5, 5.6 and Figure 5.2, the potential forthe benzene - benzene interaction is significantly stronger than the Rg - benzeneinteractions. In particular, the binding energy of the benzene dimer is almost twotimes larger than the binding energy of the Xe – benzene complex. Extendingthe density-of-states and interaction-strength arguments to molecule - moleculesystems, one should expect that the mean lifetime of benzene - benzene collisionsmust be much larger than that for Xe - benzene collisions. However, contraryto this expectation, the results presented in Figure 5.6 show that the lifetime ofbenzene – benzene collisions is much shorter than that of Xe – benzene collisions.In addition, we observe that the lifetime of the benzene – benzene collisions is aslowly varying function of the collision energy, quickly converging to a constantvalue at the collision energy 5 cm−1. The energy dependence of the mean lifetimeof the benzene – benzene collisions is similar to that of the He – benzene collisionsand is much weaker than that of the Xe – benzene collisions.In order to elucidate the dynamical details of the Rg - benzene and benzene- benzene collisions, we compute the time dependence of the separation distanceand the relative configurations of the collision partners as well as the translational,rotational and potential energy of the molecule illustrating energy exchange. Fig-ure 5.8 illustrates a typical long trajectory of the Xe – benzene collision pair withlifetime longer than 600 ps. The trajectory is characterized by the formation of thequasi-bound collision complex that lasts from 65 ps to 702 ps. During this time,there are two episodes of strong interaction between Xe and benzene, each longerthan 200 ps. The quasi-bound complex is characterized by strong energy exchangebetween the translational and rotational motions of the molecule. The episodesof strong interaction are separated by an excursion to large atom - molecule sep-arations, where the interaction energy of the Xe - benzene complex decreases to−0.22 cm−1 at R = 16.22 A˚. Such excursions to long-range parts of the interactionpotential are common for low-temperature collisions, when most (but not all) ofthe kinetic energy returns from the internal energy of the molecule to the relativeatom - molecule motion. It is interesting to observe a quasi-periodic structure ofatom - molecule separations, probably originating from the dynamical features ofthe rotational motion of benzene within the complex.79l l l l l l l l l l l l l l l l l l l l l l l l l0 5 10 15 20 250100200300Collision Energy (cm−1) Lifetime  (ps)lll l l l l l l l l l l l l l l l l l l l l l lllHeBenzeneArtificialXel l l l l l l l l l l5 6 7 8 9 100100200300Temperature (K)Lifetime (ps)l HeBenzeneNeXeFigure 5.6: Left panel: The mean lifetimes of the quasi-bound collision com-plexes of benzene with He (black circles), an artificial atom describedin the text (green crossed circles), benzene (blue stars), and Xe (pur-ple down-triangles) as functions of collision energy. Right panel: Themean lifetimes of the quasi-bound collision complexes of benzene withHe (black circles), benzene (blue stars), Ne (red squares), and Xe (pur-ple down-triangles) as functions of temperature.The interaction potential of the Xe – benzene complex is characterized by threesignificant minima with the energies −513 cm−1, −179 cm−1 and −258 cm−1,corresponding to the out-of-plane, vertex-in-plane and side-in-plane configurations(see Table 5.1 and Figure 5.2). The atom is allowed to sample the entire topologyof the potential surface hopping from one minimum to another, in a process remi-niscent of an isomerization reaction. In order to illustrate this, we present in Figure5.8 the internal configuration of the Xe – benzene collision complex by plotting theangles defined in Figure 5.1 as functions of time. We find that both the polar angleand the azimuthal angle change periodically from 0 degree to 180 degrees through-out the collision process, and vary faster when there is strong energy exchange dueto the substantial gain in the rotational energy from the potential energy.Figure 5.9 shows a short segment of the trajectory from 80 ps to 120 ps. It canbe seen that the polar angle (θ ) is correlated with the intermolecular distance to80He Ne Ar Kr Xe050100150200Reduced Mass Lifetime  (ps)He PotentialNe PotentialAr PotentialKr PotentialXe PotentialFigure 5.7: The mean lifetimes of the quasi-bound collision complexes ofbenzene with atoms of varying mass (indicated on the horizontal axis)computed with the potential surfaces for the interaction of benzene withHe (black circles), Ne (red squares), Ar (blue diamonds), Kr (green up-triangles), Xe (purple down-triangles) at a collision energy of 5 cm−1.81051015202530R  (A° )0 200 400 600 800−600−3000200500Time (ps)Energy  (cm−1 )Translational EnergyRotational EnergyPotential EnergyFigure 5.8: The time dependence of the relative distance R, the translationalenergy (black solid line), the rotational energy of the molecule (reddashed line), and the interaction potential (blue dotted line) for a col-lision of Xe with benzene at a collision energy of 2 cm−1.82345678R  (A° )060120180Polar Angle (degree)80 90 100 110 120060120180Azimuthal Angle (degree)Time (ps)Figure 5.9: The time dependence of the relative distance R, the polar angleand the azimuthal angle as defined in Figure 5.1 for a collision of Xewith benzene at a collision energy of 2 cm−1.83approach 0 or 180 degrees when R is smaller than 5 A˚. This illustrates a prefer-ence for an out-of-plane configuration at small atom - molecule separations. Theazimuthal angle tends to be uncorrelated from R and vary uniformly in the intervalbetween zero and 180 degrees, indicating that the molecule rotates freely aroundits 6-fold symmetry axis (Cν ) during the collision event. When R > 5 A˚, the out-of-plane configuration is no longer favorable, and Xe exhibits a preference forinteracting with benzene in the plane of the molecule, either in the vertex-in-plane(φ = 30,90,150) or side-in-plane (φ = 0,60,120) configuration. This is consistentwith the observations in molecular beam scattering experiments. [127].Benzene - benzene collisions are very different from collisions of Xe with ben-zene. Figure 5.10 presents the time dependence of the relative separation as well asthe translational, potential and rotational energy during a long trajectory for a ben-zene – benzene collision. As in the case of Xe - benzene collisions, the trajectory ischaracterized by multiple episodes of strong energy exchange. However, unlike inthe Xe - benzene case, these episodes are very short-lived, less than 10 ps each. Inaddition, the collision complex spends half of the collision time rambling outsidethe region of strong interaction.Figure 5.11 shows the time dependence of the relative configuration anglesfor the benzene - benzene complex. It can be seen that the polar and azimuthalangles (defined in Figure 5.3) are allowed to sample the entire range of valuesbetween 0 and 180 degrees. However, the change of the normal angle is verysmall even when there is large deposition of energy into the rotational motion ofthe molecules. The change of the normal angle describes the relative rotation ofthe molecules around their C2 symmetry axes (see Figure 5.3). By examining thefull topography of the benzene-benzene potential energy surface, we find that thechange of the normal angle at small intermolecular separations is always impededby large energy barriers, leading to hindered rotation around the in-plane C2 axes.Recent theoretical and experimental studies [128] show that this hindered rotationcauses the splitting of the rotational spectral lines of the benzene dimer.In order to show that the hindered rotation is common in benzene-benzene col-lisions, we sample a set of 400 trajectories with different initial Euler angles forthe rotational angular momentum of the molecules J = 0 and a collision energyof 2 cm−1. Figure 5.11 illustrates the dependence of the collision lifetime on the84051015202530R  (A° )0 50 100 150 200 250−100001000Time (ps)Energy  (cm−1 )Translational EnergyRotational Energy ARotational Energy BPotential EnergyFigure 5.10: The time dependence of the distance R, the translational energy(black solid line) of the relative motion, the rotational energy of thefirst molecule (red dashed line), the rotational energy of the secondmolecule (green dash-dotted line), and the interaction potential (bluedotted line) for a benzene - benzene collision with a collision energyof 2 cm−1. 85llllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllll0102030R  (A° )060120180Polar Angle (degree)060120180Azimuthal Angle (degree)0 50 100 150 200 25060708090Time (ps)Normal Angle (degree)Figure 5.11: The time dependence of the distance R, polar angle θ , azimuthalangle φ and normal angle γ as defined in Figure 3 for a benzene -benzene collision with a collision energy of 2 cm−1.86difference between the maximum normal angle and the minimum normal angleduring the collision process γmax− γmin. We find that the trajectories with largechanges in the normal angle are more likely to have longer lifetime. Frequent fail-ure to overcome the local potential barriers makes the excursion from one energyminimum to another less likely, and requires detours outside the strong interactionregion, breaking up the collision complex. Our calculations thus show that benzenemolecules, when colliding at low energies, are restricted to a part of the configura-tion space. This reduces the average lifetime of the benzene - benzene collisions.5.4 ConclusionIn this chapter, we use the classical trajectory calculations to illustrate the effect ofthe interaction strength and the geometry of the collision partners on the formationof long-lived collision complexes at low collision energies. We first compare theresults of the calculations for collisions of benzene molecules with the rare gasatoms He, Ne, Ar, Kr and Xe. The comparison illustrates that the mean lifetimeof the collision complexes increases monotonically with the strength of the Rg -benzene interactions. As the strength of the atom - molecule potential increases,the probability of energy exchange between the translational and rotational motionsof the molecule increases, thus increasing the mean lifetime of the collision events.We find that low-energy collisions of Rg atoms with benzene are generally ergodic,i.e. the atom is allowed to sample the entire potential energy surface.We next compare the results of the Rg – benzene calculations with those forbenzene - benzene collisions. Replacing the Rg atom with the benzene moleculedoubles the density of rotational states accessible at a given energy. In addition,the large polarizability of benzene leads to the interaction strength of the benzene- benzene potential about two times larger than the energy at the minimum of theXe – benzene potential. A large density of accessible states should lead to moreenergy exchange during collisions and enhanced collision lifetimes. However, weobserve that the mean lifetimes of benzene – benzene collisions are much shorterthan those of Xe – benzene collisions and are similar to those for Ne – benzenecollisions.87lllllllll llllllllllllll lll lllllll lllllllllll lllllllllllllllllllllll llllllllllll lllll llllllllllll llll lll llll llllllllllllllllllllllllllll lllllllllllllllllllllllllllllllllllll llllllll llllllllllllllllllllllllllllllllllllllllllllllllllllllllll lllllllllllll llllllllllllllllllllllllllllllll lllll llll ll lllllllll llllllllllllllllllll ll10−2 10−1 100 101 102100200300400γmax − γminLifetime (ps)Figure 5.12: The lifetime of 400 randomly sampled trajectories for benzene- benzene collisions as a function of the difference between the maxi-mum and minimum normal angles γmax− γmin sampled during the col-lision.88In order to explain this observation, we present a detailed analysis of the en-ergy exchange and configuration change during the molecule - molecule collisions.We find that the benzene - benzene collisions at low energies are significantly non-ergodic, i.e. the molecules are allowed to sample only a part of the accessibleconfiguration space. The change from one part of the configuration space to an-other must proceed through the separation of molecules to large intermoleculardistances, a process that most collision complexes cannot survive. Our results thusillustrate that the geometry of molecules plays a crucial role in the formation ofquasi-bound collision complexes at low temperatures and that the intermolecularinteraction strength cannot be used as an isolated parameter to predict the colli-sion lifetimes. Most importantly, our results indicate that the loss of moleculesdue to molecule - molecule collisions in buffer-gas cooling experiments with rigidpolyatomic molecules may be less significant than previously thought [129]. Ourcalculations are also important for the prospect of evaporative cooling of complexmolecules. The evaporative cooling experiments rely on energy transfer in elas-tic molecule - molecule collisions and are impeded by the formation of long-livedcollision complexes. Our results show that the detrimental effect of the long-livedcollision complexes may not be as prohibitive as often assumed [129].89Chapter 6Statistical background ofGaussian process models6.1 IntroductionCurrently, theoretical studies of complex physical phenomena are generally car-ried out through the following steps: first developing a mathematical model forthe process of interest, then transforming the mathematical model to a numericalmodel by implementing the former in form of a computer program and translat-ing experimental conditions and other context-specific parameters into computerinputs, and then running numerical models to produce outputs corresponding toany given input parameters. As alternatives to expensive laboratory experiments,numerical models have been successfully used to study mechanisms of complexphysical processes and for predicting the feasibility of untried experiments withmolecular gases at cold and ultracold temperatures. For example, as discussedin previous chapters, rigorous quantum calculations can be used to study colli-sional decoherence [130] or explore possible candidates for evaporative cooling inexternal field traps [131], and classical trajectory calculations [132] can be em-ployed for studying the effects of the interaction strength and the geometry ofrigid polyatomic molecules on the formation of long-lived collision complexes atlow collision energies. As the systems of interest become more complex, fromsimple diatomic molecule-structureless atom colliding pairs to complicated poly-90atomic molecule collision complexes, those numerical calculations become ex-tremely time-consuming, even on high performance computing clusters. For com-plex molecules with many degrees of freedom, it is often impossible to compute thethermally averaged cross sections, which are relevant for real experiments, sincethe evaluation of each data point is very expensive. Another problem with currentmolecular dynamics calculations is the inaccuracy of the potential energy surfaces(PES). Unfortunately, even the most sophisticated quantum chemistry calculationsproduce PES with uncertainties that lead to significant (and often unknown) errorsin the dynamical calculations. This sensitivity to PES inaccuracies is especiallydetrimental for low temperature applications (cold molecules, ultracold chemistry,astrophysics and pressure standards) [133, 134]. One possible solution is to incor-porate the inaccuracy of the PES into calculations by sampling multiple potentialsurfaces. However, it brings us back to the first problem: the computational cost.In order to solve all problems in a single step, we need to find a fast surrogatemodel for those slow numerical models. Here, we introduce a flexible and suf-ficient non-parametric model – a Gaussian Process (GP) model for this purpose.[135–139].6.2 Relevant statistical theories6.2.1 Definition and properties of multivariate normal distributionTable 6.1: Related Statistical NotationTerms Emphasis Examplesrandom variable upper case italic X ,Y,Z(non-random) number lower case italic x,y,zrandom vector boldface upper case italic X,Y,Z(non-random) vector boldface upper case italic x,y,zrandom function upper case italic F(·)1(non-random) function lower case italic f (·)matrix boldface capital A,Bestimator or estimate hat symbol θ̂ , β̂91Table 6.2: Notations for special real-valued functionsTerms Examplesprobability density function (pdf) f (X) or f (X|θ )likelihood function L (θ ) or L (θ |X)prior distribution function pi(θ )posterior distribution function p(θ |x)Let’s assume that a p-dimensional random vector Yp = (Y1,Y2, . . . ,Yp)> hasa multivariate normal distribution with mean vector µ p = (µ1,µ2, . . . ,µp)> anda positive definite covariance matrix Σ, which is denoted as Yp ∼ MVN(µ p,Σ).This means that each element Yi of a random vector Yp is normally distributed withmean µi and variance Σi,i, and the covariance of two elements Yi and Yj is Σi, j,for i = 1,2, . . . , p,and j = 1,2, . . . , p [140]. The joint density function of Yi’s isexpressed asf (Yp) = det(Σ)−1/2(2pi)−p/2exp(−12(Yp−µ p)>Σ−1(Yp−µ p))(6.1)If we partition the p-dimensional random vector Yp into two subsets as Yp =(Yq1,Yp−q2 ) with q < p, then Eq (6.1) can be rewritten asYp =(Yq1Yp−q2)∼ MVN((µ q1µ p−q2),(Σ11 Σ12Σ21 Σ22)). (6.2)whereΣi j = E{(Yi−µ i)(Y j−µ j)>}, for i = 1,2 and j = 1,2 (6.3)The marginal distribution of subset Yi is multivariate normal with mean µ i andcovariance matrix Σii:Yi ∼ MVN(µ i,Σii), (6.4)1except for cumulative distribution functions F(x) and quantile distribution functions Q(x),which by convention are non-random functions in this thesis.92and the conditional distribution of Yi given Y j still belongs to the same distributionfamily but with a shifted mean and an updated covariance matrixYi | Y j ∼ MVN(µ i +Σi jΣ−1j j (Y j−µ j),Σii−Σi jΣ−1j j Σ ji) (6.5)Properties of conditional distributions of a multivariate normal distribution will befrequently used in this and the next chapter.The simplest case of a multivariate normal distribution is a bivariate normaldistribution where a random vector Y has two elements Y1 and Y2. Figure 6.1 shows1000 Y’s randomly drawn from a bivariate normal distribution with mean(00)and covariance matrix(10 88 10)and histograms of their elements Y1 (blue) andY2 (red). In Figure 6.1, we find that most elements of Y1 and Y2 are spread from-10 and 10 with mean around 0, and their histograms are roughly consistent withnormal distribution N(0,10). We also recognize that at any vertical cross section, across section where either Y1 or Y2 is fixed, the spread of the corresponding elementYi is decreased and no longer centred around 0, which is in agreement with theproperties of conditional distributions we discussed before.6.2.2 Random function (“Process”)Since statisticians usually use process as a synonym for random function (or ran-dom field, stochastic process), we will use these terms alternately throughout thischapter. From now on, we only consider real functions whose values are real, un-less otherwise stated. A random function, or a process is a function F(·) whichcan take on a set of different possible functions f (·) with associated probabilities[141]. Or, a random function F(·) can be regarded as a mapping from outcomes ωiin a sample space Ω to a given set of functions fi(·), which is denoted byF(ωi) = fi(·). (6.6)A random function can be partially or fully described by a mean function µ(·)93●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●● ●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●● ●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●● ●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●● ●●●●●●●●●●●●●●●●●●●●●● ●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●● ●●●●●●●●●●●●●●● ●●●●●●●●●●●●●●●●●●●●●●●●●● ●●●●●●●●●●● ●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●● ●●●●●●●●●●● ●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●● ●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●● ●●●●●●●●●●●●● ●● ●●●●●●●●●●●●●●● ●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●● ●●●●● ●●●●●●●●●● ●●●●● ●●●●●●●●●●●●●●●●●●●●●●●●●●●●● ●●●●●●●● ●●●●●●●●●●●●●●●●●●●●●●●●●●● ●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●● ●●●●●●●●●●●●●●●●● ●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●● ●●●●● ●●●●●●●●●●●●●●●−10 −5 0 5 10−10−50510yY1Y2HistogramofY2Histogram of Y1Figure 6.1: Scatter plot of 1000 samples from a bivariate normal distribution.Each circle represents a sampled vector y with elements y1 and y2 plot-ted on the X-axis and Y-axis, respectively.and a covariance function (covariance kernel or autocovariance) C(·, ·). Let usconsider a simple example of a random function in the form ofF(x) = β0 +β1x+β2x2 +β3x3where the coefficients βi for i = 0,1,2,3 are independent normal random variableswith mean 1 and variance 22. Since the sum of independent normal random vari-ables is still normal, with the mean being the sum of the individual means andthe variance being the sum of the individual variances, we can easily show thatgiven x, F(x) is normally distributed with the mean function µ(x) = 1+x+x2 +x3and the variance σ2(x) = 4x6 + 4x4 + 4x2 + 4 > 0, and the covariance functionC(xi,x j) = 4(x2i x2j + 1)(xix j + 1). Figure 6.3 presents five random draws (or so-94!1!2!i............⌦F (·)f(·)f1(·)f2(·)fi(·)Figure 6.2: A random function F(·) mapping from a set Ω of outcomes ω ina sample space to deterministic functions f (·).called sample paths) from Eq (6.7), where every draw is a cubic equation. Wefind that the farther x is away from 0, the more widely the possible values of F(x)spread out, and when xi and x j have the same sign, they tend to increase or decreasetogether.The most important random functions used in nonparametric modelling aresecond-order stationary random functions. A second-order stationary random func-tion Z(·) has a constant mean µ and a constant variance σ2, and the covariancebetween two variables Z(xi) and Z(x j) is given byCov{Z(xi),Z(x j)}=C(|xi−x j|) = σ2R(|xi−x j|) (6.7)where C(·) and R(·) are the covariance and correlation functions, and they solelydepend on |xi−x j|.95−5 −3 −1 1 3 5−200−1000100200xF(x)f1(x)f2(x) f3(x)f4(x) f5(x)Figure 6.3: Five draws from a random function6.2.3 Gaussian random functionGaussian random function (GRF), or Gaussian Process (GP) is, by definition, arandom function F(·) for any choice of inputs x1, x2, . . . xn from its domain Xwhich is a subset of q-dimensional Euclidean space having positive q-dimensionalvolume, such that the outputs F(x1), F(x2) . . . F(xn) jointly follow a multivari-ate normal distribution fully described by a mean function µ(·) and a covariancefunction C(·, ·).We will frequently use second-order stationary Gaussian random functions Z(·)with zero mean in this thesis. The covariance function of a second-order stationaryGRF, as a product of constant variance σ2 and a correlation function R(·, ·) [141,142], only depends on the absolute difference between input vectors |xi− x j|. Inaddition, it must be positive definite, since the variance of any linear combination96of outputs F(xi) must be non negative.Var(n∑i=1λiF(xi))=n∑i=1n∑j=1Cov(F(xi),F(x j))≥ 0 (6.8)There are a wide variety of correlation functions for response surfaces withdifferent levels of smoothness [143–146]. One of the widely-used families of cor-relation functions are power exponential correlation functions defined byR(x,x′) = exp{−q∑i=1ωi|xi− x′i|pi}(6.9)where q is the number of input parameters x(x′), and ωi and pi represent thesmoothness parameter and the scale parameter for each input parameter.Eq (6.9) is often rewritten in terms of so-called characteristic length-scalesθi, which is physically interpretable in the same unit as the corresponding inputparameter xi,R(x,x′) = exp{−q∑i=1( |xi− x′i|θi)pi}. (6.10)A special member of the family of power exponential correlation functions isthe Gaussian exponential correlation functions, where the scale parameter p equals2 for all input parameters:R(x,x′) = exp{−q∑i=1ωi(xi− x′i)2}(6.11)Gaussian exponential correlation functions make the processes differentiable toall orders, so they lead to smooth responses. With all pi = 1 corresponding toexponential correlation functions, the process is only continuous, so it is useful fornon-differentiable response surfaces.Another popular correlation function family is the Mate´rn correlation functions97defined asR(x,x′) =q∏i=11Γ(ν)2ν−1ανi Kν(αi) (6.12)where αi = 2ωi√ν |xi− x′i| and Kν(·) is a modified Bessel function of order ν .When ν = 12 + k where k is a non-negative integer, the Mate´rn correlation func-tions have analytical expressions. We usually fix the value of ν to ν = 5/2, whichreduces Eq (6.12) toR(x,x′) ={q∏i=1(1+√5|xi− x′i|θi+ 5(xi− x′i)23θ 2i)exp(−√5|xi− x′i|θi)}. (6.13)With the Mate´rn correlation function with parameter ν , the associated processes aredifferentiable to order k (k < ν). Thus, with ν = 5/2, the GP is twice differentiable.We later find that the Mate´rn correlation function yields more accurate predictionsfor the present applications in molecular collisions than the Gaussian correlationfunctions we used in chapter Latin hypercube samplingBefore we start any numerical calculations, we should plan it carefully. Some-times, a good computation design can save us a tremendous amount of time anddrudgery. Numerical simulations often involve sampling issues: selecting inputsites at which to compute outputs. Simple random (SR) sampling is a widelyused scheme in numerical simulations. This approach uses random or pseudo-random numbers to sample from a probability distribution of each input parameter.For a common low-dimensional sampling, a SR method is a good choice sincemany computer languages have their own library functions for generating randomor pseudo-random numbers. However, for simulations involving multiple inputparameters, SR sampling is not quite efficient, asking for a large sample size toobtain converged results, since some sampled points might be very closed to eachother, failing to spread over the entire input space evenly. For example, for theBenzene-Rg collision model, the mean lifetime of collision complexes depends on11 input variables, including temperature, collision energy, mass and the 8 poten-98tial parameters. The calculation of a single mean lifetime at each 11-tuple inputconfiguration requires an average over 105 trajectories, which takes 100-300 hoursof runtime even on a high performance computing cluster. In order to accelerate thecalculation for a single mean value and fit surrogate models with a small numberof mean values, we apply a more efficient sampling scheme: the Latin Hypercubesampling (LHS) [147–149].1 2 4 3 2 3 1 4 3 4 2 1 4 1 3 2 1 3 4 2 2 4 3 1 4 2 1 3 3 1 2 4 !"#$%&'(")*% +%,%+%&(-./(%Figure 6.4: (a) a 4× 4 Latin square (b) a completed 4×4 Sudoku gridWe first introduce LHS for a two-dimensional sampling case. Suppose we wantto sample n points from a square.Step 1: Divide each axis into n equally spaced intervals, leading to n2 cells of equalsize.Step 2: Fill these cells with integers 1,2, . . . ,n to form a latin square, each integeroccurring exactly once in each row and exactly once in each column.Step 3 : Choose an integer from 1 to n at random, and then randomly select one99point from each cell filled with the selected integer, leading to a typical LHS ofsize n.If you are familiar with Sudoku puzzles, you will find that any completed Su-doku grid is a special Latin square with additional restraints on sub-squares. Figure6.4 gives an example of LHS of size 4 and a 4×4 Sudoku grid. The cells filled withthe same integer are also highlighted by the same colour. A Latin square guaranteesthat no matter which integer you pick, the centres of the chosen cells are alwaysspread evenly over each input dimension. In contrast, using SR sampling, we mightend up with some points clustered closely due to complete randomness.Now, we describe a general approach for getting a LHS sample of size n froma p-dimensional sample space X p, where each sample is a random vector withp components like Ypi = (Yi1,Yi2, . . . ,Yip)> for i = 1,2, . . . ,n. We first generate ann× p matrix, and fill p columns with p different randomly selected permutationsof {0,1, . . . ,n-1}. Next, we convert the (i, j)th matrix element, which is filledwith integers Zi j, to the jth component of ith random vector Yi j by following theoperation:Yi j = Q j(1n(Zi j + εi j)), for i = 1,2, . . . ,n and j = 1,2, . . . , p (6.14)where Q(·) j is the quantile function of Yj, the jth component of Y, and εi j is inde-pendent and identically distributed N(0,12).The above generic LHS will lead to numerous candidate samples with requiredsize due to randomly selected permutations of {0,1, . . . ,n-1} and the addition ofε . In this chapter, we use maximin LHS from R package lhs [150], which is fur-ther optimized by maximizing the minimum distance between design points fromcandidate samples [151].6.2.5 Maximum likelihood estimationMaximum likelihood estimation (MLE) is a frequentist-based method for estimat-ing the parameter(s) of a statistical model (distribution) [152]. The essence ofthe MLE method is finding the value(s) of the parameter(s) that makes the ob-served data most likely to occur. Given a probability model, its likelihood func-100tion takes the same form as its joint probability function (probability density func-tion for continuous random variables and mass probability function for discreterandom variables), but the probability function is a function of observed data asf (Y1,Y2, . . . ,Yn | θ 0), whereas the likelihood function is a function of model pa-rameter variable(s) θ conditioned on observed dataL (θ |Y1,Y2, . . . ,Yn). The MLEestimator θ̂ for the true value(s) of the parameter(s) θ 0 maximizes the likelihoodfunctionθ̂ ≡ argmaxθ∈ΘL (θ | Y1,Y2, . . . ,Yn) (6.15)In practice, it is easier to work with log-likelihood function - the natural log ofthe likelihood function logL (θ | Y1,Y2, . . . ,Yn). We will drop the additive con-stant terms in log-likelihood functions in the following analysis, unless indicatedotherwise. If the log-likelihood function is a concave function of θ , which hasat most one maximum, we take the first derivative with respect to each parameterθ1,θ2, . . . ,θn, respectively. If the value of θi leads to ∂ logL∂θi = 0, it is the estima-tor θ̂i of the true value of θi. For some complex models, log-likelihood functionscannot be maximized analytically, and we resort to numerical methods.Here, we need to emphasize that the observed data Y1,Y2, . . . ,Yn are randomvariables, which change from sample to sample, so the estimator θ̂ representsrandom variable(s) over repeated samples. For a specific sample of observationsy1,y2, . . . ,yn, θ̂ is real-valued, and called an estimate. Unfortunately, statisticiansuse the same notation θ̂ for both estimators and estimates. Since it is an almost uni-versal convention in statistical community, and we can distinguish between thembased on the context, we will keep the same notation throughout the thesis.6.3 Fitting a Gaussian process modelAll statistical inferences in this section are made within frequentist framework.The other widely-used statistical paradigm, Bayesian inference will be reviewedand used in section 6.5. Studies of statistics generally use simulators to represent anumerical model (computer experiment) and refer to a statistical surrogate modelas emulators or metamodels. The output of a typical numerical model generallycontains multiple values. To simplify our analysis, we restrict the output to be a101scalar y, while the input x can be either a scalar or a vector.6.3.1 Why use a Gaussian process model?To facilitate the studies of collision dynamics of complex molecules, it is desirableto compute a limited number of outputs at sampled input sites, and then use a surro-gate model to rapidly predict the outputs at any other input sites . The deterministicnature of most numerical models makes interpolation preferable to regression forthis prediction purpose.There are many choices of interpolation methods, including polynomial in-terpolation, spline and statistical interpolation. For low-dimensional cases, if thetraining sites are densely and uniformly distributed throughout the input domain,we will get satisfactory predictions regardless of interpolation methods. However,it is difficult to extend common interpolators, such as polynomial and spline, tohigh-dimensional nonparametric numerical models. To overcome this problem,we introduce Gaussian process models (analogous to kriging in geostatistics andspatial analysis [144, 145] due to some advantages:First, Gaussian process models are flexible, that is, they are adaptable to anyshapes of functions with up to 100 input parameters [136, 153, 154].Second, Gaussian process models are efficient. A decent fitting usually onlyrequires 10 times the number of input parameters [155].Last but not least, Gaussian process models provide prediction errors alongwith predictions for values of unknown functions [136].A Gaussian process model, as a kind of statistical interpolation, treats the un-known function f (·) representing a numerical model as a realization (sample path)of a random function F(·). The random function interpretation does not mean theunknown function f (·) is random, but rather it provides the prediction f̂ (x0) atany input site x0 with statistical meaning. The uncertainty associated with the pre-diction does not result from an actual random function but from the incompleteinformation about the underlying deterministic function. Given a limited numberof sampled data, there exit numerous functions fi(·) passing though all sampleddata Yn. At an arbitrary input site x0, the spread of all possible fi(x|Yn) counts forthe estimation variance (right panel in Figure 6.5).102XYl ll lllllllXFigure 6.5: left panel: 20 functions fi(x) randomly drawn from a randomfunction F(x); right panel: 20 functions fi(x|yn) passing through allsampled data yn (red circle) randomly drawn from the same randomfunction.To illustrate the reason for choosing a Gaussian process instead of other kindof random functions, we first drop the Gaussian assumption, and construct a non-parametric model F(x) consisting of a polynomial regression component (so-calledtrend function in geostatistics) h(x)>β and a second-order stationary random func-tion Z(x) with zero mean:F(x) =k∑j=1hj(x)βj +Z(x) = h(x)>β +Z(x) (6.16)where x=(x1,x2, . . . ,xq)> is the column vector of q input parameters, h(x) consistsof k known regression functions, called regressors, h1(x),h2(x), . . . ,hk(x) (q and kare not necessarily equal) and the location parameter β = (β1,β2, . . . ,βk)> is acolumn vector of unknown regression coefficients.The model in Eq (6.16) is equivalent to universal Kriging (UK) in geostatistics.A special case of UK is ordinary Kriging (OK) where the regression componenth(x)>β reduces to a constant β .Suppose we have already evaluated the original numerical model at input sites103x1, x2, . . . , xn, and obtained the corresponding outputs yn =( f (x1), f (x2), . . . , f (xn))>,which are realizations of random variables Yn = (F(x1),F(x2), . . . ,F(xn))>. Thenfor a new input site x0, we would like to predict a random variable Y0 = F(x0)based on data Yn. Given the properties of the non-parametric model F(x) definedin Eq (6.16), the joint distribution of Y0 and Yn is(Y0Yn)∼((h(x0)>H)β ,σ2z(1 A>0A0 A))F(6.17)where F represents a given family of distributions (not necessarily multivariateGaussian), and H is a n×k design matrix with ith row filled with the h(xi)>, σ2z isthe constant variance of Z(·), A is a n×n matrix defined asA =1 R(|x1−x2|) . . . R(|x1−xn|)R(|x2−x1|) 1....... . .R(|xn−x1|) . . . 1, (6.18)and A0 = (R(|x0 − x1|),R(|x0 − x2|), . . . ,R(|x0 − xn)|)> is also specified by thecorrelation function R(·) of Z(·).Since Y0 and Yn are dependent variables, intuitively, Yn must contain someinformation about Y0. We apply the basic idea of Kriging to predict the value of Y0by using a weighted average of the values of Yn:Ŷ0 =n∑i=1λiYi = λ>Yn. (6.19)Ŷ0 is a linear predictor of Y0. An optimal or so-called the best linear unbiasedpredictor (BLUP) of Ŷ0 should satisfy the following two criteria:First, the predictor Ŷ0 should be unbiased. That means that the expectation of Ŷ0is the same as the expectation of Y0 with respect to a given joint distribution of104(Y0,Yn):E(Ŷ0)F = E(Y0)F . (6.20)Second, the predictor Ŷ0 minimizes the mean squared prediction error (MSPE) de-fined byMSPE = E[(Y0− Ŷ0)2]F. (6.21)In order to obtain the BLUP of Y0, we need to find the minimizer λ ∗ of Eq(6.21) subject to the constraint in Eq (6.20). Based on the assumption that thecorrelation function R(·) is known, if the correlation matrix A of Yn is invertible,λ ∗ exists and is unique, and the BLUP is written asŶ0 = h(x0)>β̂ +A>0 A−1(Yn−Hβ̂), (6.22)whereβ̂ = (H>A−1H)−1H>A−1Yn (6.23)is the generalized least squares estimator of β .We notice that the BLUP coincides with the conditional mean of Y0 given Ynwhen Z(x) in Eq (6.16) is Gaussian. That is why we use a Gaussian process under-lying a non-parametric model.6.3.2 Prediction based on Gaussian process modelPrediction [135–137] is an important topic in the statistical analysis of computerexperiments and the basis of all other applications, such as the sensitivity anal-ysis [156] or the model calibration [138, 157]. The essence of prediction is theinterpolation of training data.In a Gaussian process model, Z(x) in Eq (6.16) is a Gaussian second-order sta-tionary process with zero mean and constant variance σ2z . The correlation func-tion R(·|ψ ) is parametrized by ψ , which contains all the information definingthe correlation function, and the correlation between any Z(xi) and Z(x j) solely105depends on |xi − x j| (called lag in Kriging). The joint distribution of the out-put at an untried input site Y0 = F(x0) and the outputs at training sites Yn =(F(x1),F(x2), . . . ,F(xn))>in Eq (6.17) becomes multivariate normal(Y0Yn)∼ MVN((h(x0)>H)β ,σ2z(1 A>0A0 A))(6.24)The choice of the regression component or the trend function h(x)>β reflectssome prior knowledge about the shape of the function F(·) which relates inputsto outputs. Unlike parametric models, such as the linear regression model, wheremissing or spurious regressors lead to poor fitting, a Gaussian process model, asa non-parameteric model, can suit any shape of the function F(·). Therefore, Eq(6.16) can be simplified asF(x) = β +Z(x) (6.25)where β is the overall constant mean of F(x).Since we usually use the exponential power correlation functions with thesame power parameter p for all input dimensions or the Mate´rn correlation func-tions with a fixed order of the modified Bessel function ν throughout the the-sis, the correlation function is only specified by the smoothness parameters Ω =(ω1,ω2, . . . ,ωq)> (or the characteristic length-scales θ = (θ1,θ2, . . . ,θq)>) and asingle power parameter p or an order parameter ν , leading to ψ ≡ (Ω(or θ )).According to the properties of multivariate normal distribution, the conditionaldistribution of Y0 = F(x0) given Yn follows a normal distributionY0|Yn,β ,σ2z ,Ω ∼ N(m(x0)∗,σ∗2z (x0)) (6.26)wherem(x0)∗ = h(x0)>β +A>0 A−1(Yn−Hβ ) (6.27)106σ∗2z (x0) = σ2z (1−A>0 A−1A0) (6.28)Given β and ψ , Ŷ0 = m(x0)∗ can be proved to be the BLUP [136]. It can beinterpreted as the marginal mean h(x0)>β plus a correction term, which is a linearcombination of the residuals of previous data (Yn−Hβ ).In the frequentist framework, the model’s hyperparameters β , σ2z and ψ areunknown but fixed constants. There are several methods for estimating the model’sparameters (β ,σ2z ,ψ ), such as the maximum likelihood estimation(MLE) [135],Higdon’s Bayesian method [157, 158] and Treed Gaussian process hierachical ap-proach (TGP) [159, 160]. We will use the frequentist-based-MLE.6.3.3 Maximum likelihood estimation of model parametersGiven β ,σ2z ,ψ , the probability density function of a Gaussian process model iswritten asf (Yn|β ,σ2z ,Ω) = det(A)−1/2(2piσ2z )−n/2exp[−(Yn−Hβ )>A−12σ2z(Yn−Hβ )](6.29)Correspondingly, its log likelihood function (up to an additive constant) islogL (β ,σ2z ,Ω|Yn) =−12[nlogσ2z + log(det(A))]− 12[(Yn−Hβ )>A−1σ2z(Yn−Hβ )](6.30)Given Ω, the maximum likelihood estimators of β and σ2z have closed-form solu-tions, and are expressed asβ̂ (Ω) = (H>A−1H)−1H>A−1Yn (6.31)σ̂2z (Ω) =1n(Yn−Hβ )>A−1(Yn−Hβ ) (6.32)Since both β̂ (Ω) and σ̂2z (Ω) depend on σ2z , we plug them back into Eq. (6.30),107leading tologL (Ω|Yn) =−12[nlogσ̂2z + log(det(A))+n](6.33)which is a function of Ω only. Next, we need to maximize the rewritten log-likelihood function logL (Yn|Ω) numerically, which requires repeated computa-tion of the determinant of the correlation matrix |A| and its inverse A−1. If anypair of training sites are very close together, the computation of |A| and A−1 willbe unstable due to near-singularity. To avoid near-singularity, a nugget term ε isadded to the GP model as F(·)+ε , which will be discussed in detail in section Two-step Gaussian process modelsSuppose we have a fast numerical model and a slow numerical model, both ofwhich describe the same physical process at different levels of complexity, and ob-tain a large quantity of fast model outputs y1,y2, . . . ,yN at input sites x∗1,x∗2, . . . ,x∗Nand scanty slow model output q1,q2, . . . ,qn (n N) at input sites x1,x2, . . . ,xn.Take our present scattering problems for example, the former is a fully classicalmodel, while the latter gives a rigorous quantum treatment to the same collisionevents. Generally speaking, the quantum model provides a more faithful represen-tation of the true physical process , therefore, it is desirable to construct a Gaussianprocess model for the quantum model for the studies of relationship between in-puts and outputs; that is, we treat the quantum model e(·) as a realization of aGaussian process E(·), and interpolate its outputs directly. However, we usuallycannot afford enough data points to fit a Gaussian process model with desired pre-dictive accuracy especially for high-dimensional input cases. Moreover, quantummodels are more sensitive to input variables than their classical counterparts, sothey require more training sites to build good overall models. On the other hand,the classical model f (·) can generate sufficient but less accurate data points withina short time. Since the classical model and the quantum model simulate the samephysical process, their output should be correlated to some degree.We heuristically model E(·) as a linear combination of two Gaussian processesF(·) and G(·) plus a normal random variable ε with zero mean and unknown vari-108ance λ asE(·) = ρF(·)+G(·)+ ε. (6.34)whereF(·) ∼ GP(h>1 (·)β 1,C1(·, ·|ψ 1))(6.35)G(·) ∼ GP(h>2 (·)β 2,C2(·, ·|ψ 2))(6.36)ε ∼ N(0,λ ) (6.37)The Gaussian processes F(·) and G(·) models the classical model f (·) and thediscrepancy between the quantum model and the classical model g(·), respec-tively. ρ is an unknown regression coefficient. F(·), G(·) and ε are independentof each other. We denote all hyperparameters related to the covariance functionsof F(·) and G(·) as ψ 1 and ψ 2 , respectively, and represent (β 1,β 2), (ψ 1,ψ 2) and(ρ,λ ,ψ ) collectively as β , ψ and φ .(σ21 ,Ω∗X )≡ ψ 1 (6.38)(σ22 ,ΩX)≡ ψ 2 (6.39)(ψ 1,ψ 2)≡ ψ (6.40)(ρ,λ ,ψ )≡ φ (6.41)There are two purposes for adding a normal random variable ε to Eq (6.34):First, the Gaussian process modelling is made based on the assumption of con-tinuity in the outputs of a computer model. However, due to the existence of reso-nances, especially in quantum dynamics calculations, a slight variation in the inputsite x may lead to a dramatic change in the output e(x). Those discontinuities canlead to numerical instability due to near-singularity in the covariance matrix. Thenear-singularity problem is usually solved by using a nugget term ε: adding anidentity matrix up to a multiplicative constant λ to the covariance matrix Σ asΣ′ = Σ+λ I (6.42)109However, a large ε can lead to unnecessary over-smoothing, affecting predictiveaccuracy. A linear combination of two Gaussian processes F(·) and G(·) is help-ful to maintain the balance between numerical stability and minimizing the over-smoothing of the model predictions, since two processes are expected to contributemore to the total variance of observed quantum data e(x) than a single one.Second, ε can also serve as a homogeneous/heterogeneous noise term. Al-though computer models, different from experimental observations, is technicallydeterministic, however, we sometimes cannot obtain the real value of f (x) butf (x)+ε due to the partial convergence. For example, numerical results often showsome small kinks in the data curves. The addition of a small noise term to the GPmodel is proved to improve its predictive accuracy [161]. Here, a major differencebetween the nugget effect and the noise effect is that a GP model with a nuggetterm exactly interpolates all training sites. Outside the training sites, the predic-tions based on a GP models with a nugget term and the one with a homogeneousnoise term are the same.2The combined outputs DN+n from YN = (F(x∗1),F(x∗2), . . . ,F(x∗N))> and Qn =(E(x1),E(x2), . . . ,E(xn))> jointly follow a multivariate normal distribution writtenasDN+n =(YNQn)∼ MVN(N+n)((µ 1µ 2),Σ)= MVN(N+n)(H(β 1β 2),Σ)(6.43)The design matrix for the full data set isH =(H1(D1) 0ρH1(D2) H2(D2))(6.44)2Some papers do not distinguish between the nugget term and the noise term, and call both ofthem the nugget.110whereH1(D1) =h11(x∗1) h12(x∗1) . . . h1q1(x∗1)h11(x∗2) h12(x∗2) . . . h1q1(x∗2)............h11(x∗N) h12(x∗N) . . . h1q1(x∗N)=h1(x∗1)>h1(x∗2)>...h1(x∗N)>N×q1,(6.45)H2(D2) =h21(x1) h22(x1) . . . h2q2(x1)h21(x2) h22(x2) . . . h2q2(x2)............h21(xn) h22(xn) . . . h2q2(xn)=h2(x1)>h2(x2)>...h2(xn)>n×q2,(6.46)H1(D2) =h1(x1)>h1(x2)>...h1(xn)>n×q1(6.47)The covariance matrix for the full data set is written asΣ =(Σ11(D1) Σ12(D1,D2)Σ12(D1,D2)> Σ22(D2))(6.48)whereΣ11(D1)i j =C1(x∗i ,x∗j |ψ 1) (6.49)Σ12(D1,D2)i j = ρC1(x∗i ,x j|ψ 1) (6.50)Σ22(D2)i j = ρ2C1(xi,x j|ψ 1)+C2(xi,x j|ψ 2)+δi jλ (6.51)Since fitting E(·) is carried out through two steps, we call it two-step model inthe following discussion.Step 1: We first fit F(·) completely based on classical data yN by using MLEmethod or Bayesian inference introduced in section 6.5, and then substitute the111corresponding estimates β̂ 1 and ψ̂ 1 into Eq (6.43).Step 2: Given β , φ and YN , the conditional distribution of Qn still follows amultivariate distribution with shifted mean vector µ ′2(D2) and updated covariancematrix Σ′22(D2)µ ′2(D2) = ρH1(D2)β 1 +H2(D2)β 2 +Σ12(D1,D2)>Σ11(D1)−1(YN−H1(D1)β 1),(6.52)Σ′22(D2) = Σ22(D2)−Σ12(D1,D2)>Σ11(D1)−1Σ12(D1,D2). (6.53)The likelihood function of unknown hyperparameters for β 2, ψ 2, ρ and λbecomesL(β 2,ψ 2,λ |YN ,Qn, β̂ 1, ψ̂ 1)∝ det(Σ′22)− 12exp[− 12(Qn−µ ′2)>Σ′−122(Qn−µ ′2)]= det(Σ′22)− 12exp[− 12(Q′n−H2(D2)β 2)>Σ′−122(Q′n−H2(D2)β 2)](6.54)whereQ′n = Qn−ρH1(D2)β 1−Σ12(D1,D2)>Σ11(D1)−1(YN−H1(D1)β 1)(6.55), and the corresponding log-likelihood function is written as (up to an additiveconstant)logL(β 2,ψ 2,λ |YN ,Qn, β̂ 1, ψ̂ 1)=−12log(∣∣Σ′22∣∣)−12(Q′n−H2(D2)β 2)>Σ′−122(Q′n−H2(D2)β 2)(6.56)Given ψ 2, ρ and λ , the maximum likelihood estimators of β 2 have the closed-formsolution asβ̂ 2(ψ 2,ρ,λ ) =(H2(D2)>Σ′−122 H2(D2))−1H2(D2)>Σ′−122 Q′n (6.57)112By substituting β̂ 2(ψ 2,ρ,λ ) into Eq (6.56), we can numerically minimize the fol-lowing function with respect of ψ 2, ρ and λlog(∣∣Σ′22∣∣)+(Q′n−H2(D2)β̂ 2)>Σ′−122(Q′n−H2(D2)β̂ 2). (6.58)The minimizers ψ̂ 2, ρ̂ and λ̂ are their corresponding MLE estimators.Once we have the fitted two-step model E(·), we can interpolate the limiteddata from the quantum model to predict the output of e(x0) at an untried input sitex0. The unknown E(x0) together with combined outputs DN+n jointly follow a newmultivariate normal distribution(E(x0)DN+n)∼ MVN(N+n+1)((h0(x0)>H)β ,Σ′)(6.59)whereh0(x0) =(ρh1(x0) h2(x0))>(6.60)β̂ =(β̂ 1 β̂ 2)>(6.61)Σ′ =(σ20 t0(x0)>t0(x0) Σ̂)(6.62)The (N +n)×1 covariance vector of E(x0) with DN+n is defined ast0(x0) =(t01(x0) t02(x0))>(6.63)wheret01(x0)i = ρC1(x0,x∗i ) i = 1,2, . . . ,N (6.64)t02(x0) j =ρ2C1(x0,x j)+C2(x0,x j) if ε is a noise termρ2C1(x0,x j)+C2(x0,x j)+λδ|x0−x j|(0) if ε is a nugget termj = 1,2, . . . ,n(6.65)113, and δ|xi−x j|(0) is a Dirac measure defined byδ|xi−x j|(0) =0 |xi−x j| 6= 01 |xi−x j|= 0(6.66)From Eq (6.65), we find that the model with a nugget term still interpolate alltraining sites xi, since λ appears in the covariance vector t0(x0).Given DN+n, the conditional mean of E(x0) is written asm(E(x0)|DN+n)= h0(x0)>β̂ + t0(x0)>(Σ′)−1(DN+n−Hβ̂), (6.67)and it serves as the two-step model predictor of for the unknown e(x0).6.5 Bayesian calibration of numerical modelsIn this section, we discuss the use of bayesian calibration to make inference aboutunknown model parameters and provide accurate predictions based on dynamicalcalculations plus a few physical observations. The traditional parameter estimationis usually an ad hoc search for the best-fitting values, which not only requires agreat many calculations but also takes no account of possible uncertainties in com-puter models, such as model inadequacy and observation errors. Due to the pres-ence of those uncertainties, the numerical model with the best-fitting values mightstill fail to make satisfactory predictions in the future. In contrast, the Bayesianapproach allows for all sources of uncertainties such as measure error, model inad-equacy and scaling factor.6.5.1 Introduction to Bayesian inferenceThere are two main approaches to statistical inference: frequentist and Bayesian,each of which offers some advantages over the other. The methods we have dis-cussed so far are all frequentist (particularly the approach based on the MaximumLikelihood) since they are easy to use and both theorists and experimentalists arefamiliar with some frequentist concepts, like the confidence intervals and ofteninclude them in data analyses, while their Bayesian analogue, the credibility inter-114vals, are rarely used outside the statistics and biostatistics community. In frequen-tist theory, we consider that the parameter θ of a given model (distribution) is afixed constant, albeit unknown, and the data to be random. That is why we shouldnever interpret our computed 95% confidence interval as the interval encompassingthe true value of parameter θ with probability 95%. On the contrary, the Bayesianapproach can address the uncertainty of the model parameter(s) more directly andassign a probability distribution to θ based on the observed data since it treatsparameters to be random, but the data fixed. Bayesian inference is a process ofupdating our initial belief about something with objective new information to get anew and improved belief. We will intensively use Bayesian estimation method inmodel calibration since some model parameters have a meaning only over a limitedrange, and we can incorporate this information by confining parameters to physicalregions in order to ensure better estimations.The basis for Bayesian inference is Bayes’ theorem, often called Bayes’ Rule.Bayes’ theorem relates the conditional probability Pr(A|B) of event A given anotherevent B to the reversely conditional probability Pr(B|A) by the expressionPr(A|B) = Pr(B|A)Pr(A)Pr(B)= Pr(A∩B)Pr(A∩B)+Pr(A¯∩B) (6.68)where Pr(A∩B) is the probability that event A and event B both occur, and A¯ is A’scomplement.Take a spam filter for example. An incoming email containing certain suspi-cious strings might be a spam. Suppose we receive a new email containing “Con-grats”, we want to know the probability that the email is a spam by using Bayes’theoremPr(Spam|Congrats) = Pr(Congrats|Spam)Pr(Spam)Pr(Congrats)= Pr(Congrats|Spam)Pr(Spam)Pr(Congrats|Spam)Pr(Spam)+Pr(Congrats|non-Spam)Pr(non-Spam)The Bayesian treats probability as beliefs, not frequencies. To apply Bayes’115theorem to Bayesian inference, we replace event B with observed data y, event Awith model parameter θ , and rewrite Eq. (6.68) asp(θ |y) = pi(θ) f (y|θ)f (y) (6.69)wheref (y) =∑θpi(θ) f (y|θ) (6.70)if θ is discrete, orf (y) =∫θpi(θ) f (y|θ)dθ (6.71)if θ has continuous values.pi(θ) is the prior (before data collection) distribution for θ , representing ourinitial uncertainty, f (y|θ) is the likelihood of y given θ , and p(θ |y) is the pos-terior (after data collection) distribution representing our final uncertainty aboutθ after taking both initial belief and observations into account. When θ is high-dimensional and continuous, the explicit integration in f (y) = ∫θ pi(θ) f (y|θ)dθis often mathematically intractable, however, it is just a normalization constant.Since the posterior p(θ |y) is a probability density function for θ , we only need toconsider about terms involving θ . We can rewrite the posterior asp(θ |y) ∝ pi(θ) f (y|θ) (6.72)If the posterior has a standard distribution whose normalization constant is al-ready known, we can conduct the posterior inference directly and compute poste-rior quantities of interest, such as posterior mean, posterior variance and Bayesiancredibility intervals. On the other hand, when the posterior models are too com-plicated to work with directly, we use the Markov Chain Monte Carlo (MCMC)method to estimate the posterior distributions by iterative sampling.Here, we will apply Bayesian techniques for fitting and calibrating Gaussianprocess models because they are able to quantify the estimation uncertainty.1166.5.2 Metropolis-Hastings MCMC algorithmMarkov Chain Monte Carlo (MCMC) Method is nowadays widely used for numer-ical approximations for high-dimensional probability distributions. The main ideaof MCMC is to construct a Markov chain whose stationary distribution convergesto the target distribution. There are several MCMC algorithms, such as the Gibbssampler, Metropolis algorithm and Metropolis-Hastings (M-H) algorithm. We willfocus on the M-H algorithm since it is simple and applicable to any complicatedtarget distributions. I outline the details of a generic M-H algorithm belowStep1 : Initialize θ0 at some arbitrary value.Step2 : Let θi−1 be the current value of θ , and make a new draw θnew from aproposal distribution q(θi|θi−1) centred on θi−1,θnew ∼ q(θi|θi−1) (6.73)Step3 : Compute the M-H acceptance ratioα = min{1, pi(θnew) f (y|θnew)q(θi−1|θnew)pi(θi−1) f (y|θi−1)q(θnew|θi−1)}(6.74)If the proposal distribution is symmetric, satisfying q(θi−1|θnew) = q(θnew|θi−1),the acceptance ratio reduces toα = min{1, pi(θnew) f (y|θnew)pi(θi−1) f (y|θi−1)}(6.75)Step4: Simulate a uniform random number ui from [0,1]: ui ∼ Uniform(0,1) Setθi = θnew, if α > ui, otherwise reject the new draw and set θi = θi−1Step5: Iterate steps 2-4.In this way, the chain starts from an arbitrary value and moves towards thetarget distribution. To eliminate the effect of the starting point, we delete the earlyiterations of the simulation runs. This process is known as the burn-in.6.5.3 Bayesian modelling with Gaussian processesIn the Bayesian framework, the model’s hyperparameters β , σ2z and Ω are not fixedanymore, and each of them is assigned with a proper prior distribution pi(·). The117power parameter p is fixed at 2. There are several choices of prior distributions forβ , σ2z and Ω. Here, I only list one of them [157, 158]pi(σ2z ) ∝1σ2z(6.76)pi(βi) ∝ 1 for i = 1,2 . . . ,k (6.77)pi(ωi) ∝1ωifor i = 1,2 . . . ,q (6.78)where k is the number of regressors, and q is the input dimension, as mentionedbefore.Their priors are independent of each other.pi(β ,σ2z ,ψ ) = pi(β )pi(σ2z )pi(Ω) (6.79)We are unable to deal with the uncertainty in the smoothness parameter Ωanalytically, thererfore, we first integrate β and σ2z out from the joint posteriorp(β ,σ2z ,Ω|yn), and then use the posterior mode of p(Ω|yn) as point estimates forΩ.We first fix Ω, and write the likelihood function of β and σ2z asf (yn|β ,σ2z ) = det(A)−1/2(2piσ2z )−n/2exp[− 12σ2z(yn−Hβ )>A−1(yn−Hβ )](6.80)Using Bayes’ theorem, the joint distribution of β and σ2z conditional on yn has anormal inverse gamma distribution asp(β ,σ2z |yn) ∝ pi(β )pi(σ2z ) f (yn|β ,σ2z )∝ σ−(n+2)z exp{− 1σ2z(β − β̂ )>H>A−1H(β − β̂ )+(n− k−2)σ̂2}(6.81)118whereβ̂ = (H>A−1H)−1H>A−1yn (6.82)σ̂2 =(yn)>{A−1−A−1H(H>A−1H)−1H>A−1}ynn− k−2 (6.83)Since the correlation matrix A is specified by Ω, the joint distribution of β and σ2zis also conditional on Ω.According to the properties of a normal inverse gamma distribution, we canderive thatβ |σ2z ,Ω,yn ∼ MVN(β̂ ,σ2z (H>A−1H)−1)(6.84)σ2z |Ω,yn ∼ (n− k−2)σ̂2χ−2n−k (6.85)To estimate Ω, we write the full joint posterior distribution using Bayes’ theo-remp(β ,σ2z ,Ω|yn) ∝ pi(β )pi(σ 2z )pi(Ω) f (yn|β ,σ2z ,Ω) (6.86)and then analytically integrate out β and σ2z to finally get the marginal posteriordistribution of Ω and use its mode as the point estimate of Ωp(Ω|yn) ∝ pi(Ω)(σ̂2z )−n−k2 det(A)−1/2det(H>A−1H)−1/2, (6.87)whereσ̂2z =1n− k (yn−Hβ̂ )>A−1(yn−Hβ̂ ). (6.88)In addition, we can numerically simulate the posterior distribution of Ω to approx-imate the posterior mean, the posterior variance and credibility intervals.Given outputs at training sites and model hyperparameters β , σ2z and Ω, thepredictive distribution of the output Y0 at the untried site x0 follows a normal dis-tribution. By integrating out β and σ2z , we obtain that the predictive distribution of119Y0|yn is a non-central t distribution with n− k degrees of freedom.Y0|yn ∼ tn−k(m(x0)∗∗,σ∗∗2z (x0)) (6.89)wherem(x0)∗∗ = h(x0)>β̂ +A>0 A−1(yn−Hβ̂ ) (6.90)σ∗∗2z (x0) = σ̂2z (1−A>0 A−1A0)+ σ̂2z{(h(x0)−H>A−1A0)>(H>A−1H)−1(h(x0)−H>A−1A0)}(6.91)Thus, the Bayesian estimator of Y0 isŶ0(Ω̂) = h(x0)>β̂ (Ω̂)+A>0 A−1(yn−Hβ̂ (Ω̂))(6.92)where Ω̂ is the posterior mode of p(Ω|yn)6.5.4 Bayesian model calibrationIn our calibration model, the inputs are divided into two groups : control vari-ables and calibration variables. Control variables x = (x1,x2, . . . ,xq1)> are speci-fied in both lab experiments and computer simulations, and calibration variablest = (t1, t2, . . . , tq2)> are substitutes for the unknown true calibration parametersθ = (θ1,θ2, . . . ,θq2)>, which are fixed in the real physical process, and will belearned by fitting computer outputs to experimental observations.We denote the true value of the physical process as ζ (x), which depends oncontrol input x alone since the true calibration parameters are fixed. The experi-mental observations Z(x) are generally subject to measurement errors εZ(x) = ζ (x)+ ε (6.93)where ε’s are identical and independently distributed as N(0,λ ).Suppose we have n field observations performed under n different controlled120conditions, we denote them as z = (z1,z2, . . . ,zn), where zi is for ith observation atknown variable input xi. In addition, we have N (much larger than n) runs of thecomputer code, and denote their outputs as η(xi, ti). The relationship between theexperimental observations and code outputs in the general cases can be expressedasZ(x) = ζ (x)+ ε = ρη(x,θ )+δ (x)+ ε (6.94)where ρ is an unknown scaling parameter and δ (·) is a model inadequacy function.We treat δ (·) also as a Gaussian random function, which is independent ofthe code output η(·, ·). As Gaussian random functions, δ (·) and η(·, ·) are bothcharacterized by their own mean functions and covariance functions:η(·, ·)∼ GP ( m1(·, ·),c1((·, ·),(·, ·)) ) (6.95)δ (·)∼ GP ( m2(·),c2(·, ·) ) (6.96)wherem1(x, t) = h1(x, t)>β 1 (6.97)m2(x) = h2(x)>β 2 (6.98)We use Gaussian exponential correlation functions for both η(·, ·) and δ (·). Sincethe power parameters are fixed at 2, the covariance functions are only specified bythe variance σ2 and smoothness parameters Ω = (ω1,ω2, . . . ,ωq), and expressedasc(x,x′) = σ2exp(−q∑i=1ωi(xi− x′i)2) (6.99)We denote all hyperparameters related to the covariance functions of η(·, ·) andδ (·) as ψ 1 and ψ 2 , respectively, and represent (ψ 1,ψ 2) and (ρ,λ ,ψ ) collectivelyas ψ and φ .121(σ21 ,Ω∗X ,Ωt)≡ ψ 1 (6.100)(σ22 ,ΩX)≡ ψ 2 (6.101)(ψ 1,ψ 2)≡ ψ (6.102)(ρ,λ ,ψ )≡ φ (6.103)where ΩX (or Ω∗X ) represents hyperparameters for control variables x (or x∗), andΩt represents hyperparameters for calibration variables t.The set of points at which the computer code run is denoted as design D1 ={(x∗1, t1), . . . ,(x∗N , tN)}, while the experimental measurement is carried out at pointsD2(θ ) = {(x1,θ ),(x2,θ ), . . . ,(xn,θ )}The full set of calibration data, including both code outputs and the experimen-tal observations, should jointly follow a multivariate normal distribution asDN+n =(YNZn)∼ MVNN+n (H(θ )β ,C(θ )) (6.104)The design matrix for the full data set isH(θ ) =(H1(D1) 0ρH1{D2(θ )} H2(D2))(6.105)whereH1(D1) =h1(x∗1, t1)>h1(x∗2, t2)>...h1(x∗N , tN)>N×q1,H2(D2) =h2(x1)>h2(x2)>...h2(xn)>n×q2,(6.106)H1{D2(θ )} =h1(x1,θ )>h1(x2,θ )>...h1(xn,θ )>n×q1(6.107)122The covariance matrix for the full data set is written asC(θ ) =(V1(D1) ρC1{D1,D2(θ )}>ρC1{D1,D2(θ )} λ In +ρ2V1{D2(θ )}+V2(D2))(6.108)whereV1(D1)i j = c1{(x∗i , ti),(x∗j , t j)} (6.109)V2(D2)i j = c2(xi,x j) (6.110)C1{D1,D2(θ )}i j = c1{(xi,θ ),(x∗j , t j)} (6.111)V1{D2(θ )}i j = c1{(xi,θ ),(x j,θ )} (6.112)The priors for the location parameters β and the true calibration parameter θand φ are independent of each other.pi (θ ,β ,φ ) ∝ pi(θ )pi(β )pi(φ ) (6.113)The purpose of our calibration is to make inference about the calibration pa-rameters θ based on its marginal posterior distribution p(θ |d). That is to updateour prior belief of values of true calibration parameters by combined informationfrom computer outputs and field observations.Based on Bayes’ theorem, the joint posterior p(θ ,β ,φ |d) is proportional to theproduct of priors and the likelihood.p(θ ,β ,φ |d) ∝ pi(θ )pi(β )pi(φ ) f (d|θ ,β ,φ ) (6.114)∝ pi(θ )pi(β )pi(φ )|V(θ )|−1/2exp{−12[d−H(θ )β ]>V(θ )−1[d−H(θ )β ]}(6.115)123We assume an uninformative prior on the location parameters β aspi(β 1,β 2) ∝ 1 (6.116)which is a conjugate prior for the likelihood f (d|θ ,β ,φ ), and get the posterior ofβ conditional on θ ,φ ,dβ |θ ,φ ,d ∼ N(β̂ (θ ),W(θ )) (6.117)whereβ̂ (θ ) =W(θ )H(θ )>V(θ )−1d, (6.118)W(θ ) =(H(θ )>V(θ )−1H(θ ))−1(6.119)By using Eq. (6.117), we can analytically integrate β out from the full joint poste-rior distribution to get p(θ ,φ |d)p(θ ,φ |d) =∫βp(θ ,β ,φ |d)dβ ∝ pi(θ )pi(φ )|V(θ )|−1/2W(θ )|1/2×exp{−12[d−H(θ )β̂ (θ )]>V(θ )−1[d−H(θ )β̂ (θ )]}(6.120)Ideally, the marginal posterior p(θ |d) should be obtained by integrating φ outfrom the joint distribution p(θ ,φ |d). However, the integration is mathematicallyintractable. Instead, we derive the posterior modes φ̂ of hyperparameters φ ,andplug them into p(θ ,φ |d), leading to the posterior given both d and the estimatedvalues of φ̂ . Finally, we make inference about the true calibration parameters basedon p(θ |d, φ̂ ).The model calibration is carried out through three stages:Stage 1: The complete set of φ includes ψ 1, ψ 2, ρ and λ . Since code outputsyN only depend on ψ 1, we can first estimate hyperparameters ψ 1 for the correlationfunctions c1((·, ·),(·, ·)) of computer model η(x, t) by using maximum likelihoodestimation or Bayesian methods.124Stage 2: Given estimates for ψ̂ 1 obtained from stage 1 plus the experimentalobservations and code outputs, we derive the joint posterior for the rest of theunknown hyperparameters. Based on Bayes’ theorem, we can writep(β 2,ρ,λ ,ψ 2|zn,yN , ψ̂ 1) ∝ pi(β 2,ρ,λ ,ψ 2) f (zn,yN , ψ̂ 1|β 2,ρ,λ ,ψ 2) (6.121)Since yN and ψ̂ 1 are independent of β 2, ρ , λ and ψ 2, we can rewrite the jointposterior asp(β 2,ρ,λ ,ψ 2|zn,yN , ψ̂ 1) ∝ pi(β 2,ρ,λ ,ψ 2) f (zn|β 2,ρ,λ ,ψ 2,yN , ψ̂ 1) (6.122)Instead of directly using the distribution of zn|β 2,ρ,λ ,ψ 2,yN , ψ̂ 1, which is math-ematically intractable, we assume zn|β 2,ρ,λ ,ψ 2,yN , ψ̂ 1 to have a multivariatenormal distribution as zn|y,β 2,φ ,θ . The mean vector and the elements in thecovariance matrix are approximated by integrating θ out from conditional meanE(Z|y,β 2,φ ,θ ) and conditional covariance matrix V′2(Z|y,β 2,φ ,θ ). Finally, weintegrate out β 2 to obtain the approximationp(ρ,λ ,ψ 2|d, ψ̂ 1) ∝ pi(ρ,λ ,ψ 2)|V2′ |−1/2|W2|1/2×exp{12(z−H2(D2)β̂ 2−ρη̂((D2))>|V2′ |−1(z−H2(D2)β̂ 2−ρη̂((D2))}(6.123)whereβ̂ 2 = W2H2(D2)>V2′−1(z−ρη̂((D2)), (6.124)W2 = (H2(D2)>V2′−1H2(D2))−1 (6.125)η̂((D2)i = E(η(xi,θ |y)) (6.126)Stage 3: CalibrationConditioned on the estimated hyperparameters φ̂ , we regard the posterior dis-125tribution of the calibration parameters to bep(θ |φ ,d) ∝ |V(θ )|−1/2|W(θ )|1/2exp{−12(d−H(θ )β̂ )>V(θ )−1(d−H(θ )β̂ )}pi(θ )(6.127)Then we use the MCMC numerical simulations to make inference about θ , suchas MCMC methods.The purpose of calibration is to use the calibrated model for predicting the realprocess in the future. The problem of predicting the true process Z(x) at somespecified variable inputs x can be seen as interpolating the function Z(·).The posterior distribution of Z(·) conditional on the estimated hyperparametersφ̂ and the calibration parameters θ is a Gaussian process. Its mean function is givenbyE(z(x0|θ ,φ ,d)) = h(x0,θ )>β̂ (θ )+ t(x0,θ )>V(θ )−1(d−H(θ )β̂ ), (6.128)whereh(x0,θ ) =(ρh1(x0,θ )h2(x0))(6.129)t(x0,θ ) =(ρV1((x0,θ ),D1)ρ2V1((x0,θ ),D2(θ ))+V2(x0,D2))(6.130)Its covariance function is given bycov(z(x0),z(x′0)|θ ,φ ,d) = ρ2c1((x0,θ ),(x′0,θ ))+ c2(x0,x′0) (6.131)− t(x0,θ )>V(θ )−1t(x′0,θ )+ (h(x0,θ )−H(θ )>V(θ )−1t(x0,θ ))>W(θ )(h(x′0,θ )−H(θ )>V(θ )−1t(x′0,θ ))Integrating out θ with respect to its posterior distribution, we can make inferencesabout Z(x) by using numerical methods.126Chapter 7Gaussian Process Model forCollision Dynamics of ComplexMoleculesIn this chapter, we show that a Gaussian Process model can be combined with asmall number of scattering calculations to provide an accurate multi-dimensionaldependence of scattering observables on the experimentally controllable param-eters (such as the collision energy, temperature or external fields) as well as thepotential energy surface parameters. This can be used for solving the inverse scat-tering problem, the prediction of collision properties of a specific molecular sys-tem based on the information for another molecule, the efficient calculation ofthermally averaged observables and for reducing the error of the molecular dy-namics calculations by averaging over the potential energy surface variations. Weshow that, trained by a combination of classical and quantum dynamics calcula-tions, the model provides an accurate description of the scattering cross sections,even near scattering resonances. In this case, the classical calculations stabilize themodel against uncertainties arising from wildly varying correlations of resonantlyenhanced quantum results.1277.1 IntroductionThe outcome of a molecular collision in a thermal gas is a complicated functionof many parameters: the collision energy, the internal energy of molecules, thedetails of the intermolecular potential energy surface (PES) and the relative orien-tation of the colliding species. Providing a mapping between these parameters andthe scattering observables is a very complex task, paramount to solving the inversescattering problem [162, 163] and critical for the interpretation as well as the phys-ical interpolation or extrapolation of experimental data. Knowing the dependenceof observables on the PES is also critical for assessing the accuracy of the scatter-ing calculations and averaging the dynamical results over PES variations to reducethe uncertainties [133, 134]. For simple systems, this mapping can be providedby a series of classical [164] or quantum scattering [165] calculations on a grid ofinternal energy, collision energy, relative angular momentum and/or the PES pa-rameters [133, 134]. This approach quickly becomes impossible as the complexityof molecules increases. Even for simple molecules, the role of the potential is of-ten examined by scaling the PES by a single multiplicative factor, which does notchange the topology of the PES [90]. For polyatomic molecules and to explore theeffect of the PES topology, it is necessary to examine the observables as functionsof all (usually many) parameters determining the PES. For example, a realistic PESfor the interaction of benzene with structureless atoms must be parametrized by≥ 8constants [48, 121] and it is clearly impossible to carry out dynamical calculationson an 8-dimensional grid of the PES parameters.An alternative approach could be to develop an analytical function of the scat-tering observables on the underlying parameters by fitting the calculated values,for example, with spline functions or regressions [166]. However, any such fitwould require many data points as a function of every parameter to infer the properanalytical dependence. This becomes impossible as the complexity of moleculesincreases and when the collision dynamics is affected by resonances leading todramatic variations of observables, extremely sensitive to PES [70, 167]. Here, wepropose an approach combining a small number of scattering calculations with aGaussian Process (GP) model [135, 137]. Widely used in engineering technologies[157, 158], the GP model learns from correlations between data points and pro-128vides a non-parametric dependence of observables on the underlying parameters.We show that, trained by the classical trajectory (CT) and/or quantum scatteringcalculations, the GP model provides an efficient and accurate model of molecu-lar collisions, giving the simultaneous dependence of the scattering observableson all of the underlying parameters, without the need to fit any data by analyticalfunctions. A combination of the CT and quantum calculations provides an accu-rate GP model of the scattering cross sections, even near scattering resonances. Inthis case, the CT calculations stabilize the model against uncertainties arising fromwildly varying correlations of resonantly enhanced quantum results.To illustrate the accuracy and efficiency of the model we consider the scat-tering cross sections and the formation of long-lived collision complexes of ben-zene molecules placed in a cold environment of rare gas (Rg = He – Xe) atoms.Collisions of complex organic molecules with Rg atoms at low temperatures haverecently received much attention due to the buffer gas cooling experiments [30,31, 168]. In these experiments, molecules are thermalized by momentum-transfercollisions with the buffer gas atoms, while the formation of long-lived collisioncomplexes leads to clustering and impedes the thermalization. We show that theGP model can be used to obtain the dependence of the collision lifetimes on therotational temperature and the collision energy of the molecule, the mass of the Rgatom as well as on the details of the PES for atom - molecule interactions.7.2 GP model of scattering observablesWe consider a scattering observable O as a function of q parameters described byvector x. The components of the vector x = (x1,x2, · · · ,xq)> can be the collisionenergy, the internal energy and/or the parameters representing the PES. We assumethat O is known from a classical or quantum dynamics computation at a smallnumber of x values. Our first goal is to construct an efficient model that, given afinite set of O(x), produces a global dependence of the scattering observable on x.If the same observable is measured as a function of some parameters xi – e.g., thecollision energy – we show that the model can be adjusted to produce the globaldependence of O on x that reproduces the experimental data, even if the dynamicalcalculation method is inaccurate.129We assume that the unknown function underlying scattering observable f (·)in this multi-dimensional parameter space is a realization of a Gaussian processF(·), fully characterized by a mean function µ(·), a constant variance σ2 and acorrelation function R(·, ·). Individually, the output of F(·) at an arbitrary x is nor-mally distributed with mean µ(x) and variance σ2. Jointly, the outputs at differentpoints follow a multivariate normal distribution completely defined by µ(·), σ2,and R(·, ·) [141, 142]. We assume the following form for the correlation function[143–146]:R(x,x′) = exp{−q∑i=1ωi|xi− x′i|p}0 < p≤ 2 (7.1)where p reflect the smoothness of realizations of a GP process: larger p givessmoother responses.We write the model asF(x) =k∑j=1h j(x)β j +Z(x) = h(x)>β +Z(x), (7.2)where h(x)= (h1(x), ...,hk(x))> is a vector of k regression functions, β =(β1,β2, · · · ,βk)>is a vector of unknown coefficients, and Z(·) is a second-order stationary Gaussianrandom function with zero mean. The problem is thus reduced to finding β , p andΩ = (ω1,ω2, · · · ,ωq)>.We spread n input vectors x1, ...,xn evenly throughout a region of interest andcompute the desired observable O at each xi with a classical or quantum dynamicsmethod. The outputs of a GP at these points Y n =(F(x1),F(x2), · · · ,F(xn))>follow a multivariate normal (MVN) distributionY n ∼ MVN(Hβ ,σ2A) (7.3)with the mean vector Hβ and the covariance matrix σ2A. Here, H is a n×k design130matrix with ith row filled with h(xi)> , and A is a n×n matrix defined asA =1 R(x1,x2) · · · R(x1,xn)R(x2,x1) 1....... . .R(xn,x1) · · · 1(7.4)Given Ω, the maximum likelihood estimators (MLE) of β and σ2 have closed-form solutions [135]:β̂ (Ω) = (H>A−1H)−1H>A−1Yn (7.5)σ̂2(Ω) = 1n(Yn−Hβ )>A−1(Yn−Hβ ) (7.6)To find the MLE of Ω, we fix p and maximize the log-likelihood function (up to aadditive constant)logL (Ω|Y n) =−12[nlogσ̂2 + log(det(A))+n](7.7)numerically by an iterative computation of the determinant |A| and the matrix in-verse A−1.The goal is to make a prediction of the scattering observable at an arbitraryx = x0. Because the values Y0 = F(x0) at x0 and the outputs at training sites Y n =(F(x1),F(x2), · · · ,F(xn))>are jointly distributed, as(Y0Y n)∼ MVN((h(x0)>H)β ,σ2(1 A>0A0 A))(7.8)where A0 = (R(x0,x1),R(x0,x2), · · · ,R(x0,xn))>.The conditional distribution of Y0 = F(x0) given the values Y n and hyperpa-rameters β , σ2 and Ω is a normal distributionY0|Y n,β ,σ2,Ω ∼ N(m(x0)∗,σ∗2z (x0)) (7.9)131with the conditional mean and variance given bym(x0)∗ = h(x0)>β +A>0 A−1(Y n−Hβ ) (7.10)σ∗2(x0) = σ2(1−A>0 A−1A0), (7.11)By substituting the MLE estimators β̂ , σ̂2 and Ω̂ into Eq. (7.10), we can getthe GP model prediction for the value of the scattering observable at x0.7.3 Gaussian process for interpolation of scattering dataTo illustrate the applicability and accuracy of a GP model, we first compute thecollision lifetimes of benzene molecules with Rg atoms [48, 132]. This is a criticalparameter that determines the probability of benzene molecules to collect atomsinto molecule-centered clusters by three-body recombination. We use the classi-cal trajectory (CT) method described in Ref. [48]. As shown in Ref. [121], theC6H6 - Rg interaction potentials can be expressed as a sum over terms describingthe interaction of Rg with the C-C and C-H bond fragments, characterized by 8parameters.7.3.1 Mean lifetime dependence analysis for the C6H6 – Ar systemWe first fix the PES parameters to describe the C6H6 - Ar system and focus on thedependence of the lifetimes on two parameters: the collision energy E from 2 to10 cm−1 and the rotational temperature Tr from 5 to 10K. Loeppky et al. [155]suggest that 10 times the input dimension (n = 10×d) is a good rule of thumb fordetermining the size of a space-filling design for building a good GP model whenthe response surface is smooth. To fit the GP model, we choose a 20-point and a50-point maximin Latin hypercube designs implemented in the R package lhs, andcompute the mean lifetimes at those training sites. Figure 7.1 shows the results ofthe CT calculations for the 50-point design. Of those data, only collision energyappears to have a strong negative relationship with the mean lifetime, that is, themean lifetime decreases with increasing collision energy, which is consistent withour previous work (See Chapter 5). On the other hand, the plot of the mean lifetimeagainst rotational temperature exhibits a fairly random pattern, and it is hard to tell132whether they are negatively correlated or not, and we still lack the knowledge oftheir joint effects on the lifetime. In addition, the vastly different gradients of the Trand E dependence may make the conclusions based on calculations at fixed valuesof one of the parameters misleading.lll lllll llllll lllllllllllllll ll llll lllllllll llllll5 6 7 8 9 10100150200Rotational Temperature (K)Lifetime  (ps)(a)lllll lllllllll lllllllll llllllllllll lllllllllllllll2 3 4 5 6 7 8 9 10Collision Energy (cm−1)(b)Figure 7.1: (a) and (b): The lifetime dependence on the rotational tempera-ture and collision energy for C6H6 – Ar collisions.To demonstrate the accuracy and convenience of the model, we use a 20-pointdata set and a 50-point data set to fit the GP model, separately. The regression partof the model is set to be a constant as F(Tr,E) = β +Z(Tr,E). The log-likelihoodfunctions are optimized by a multi-gradient based search algorithm implementedin R package GPfit. To verify the goodness-of-fit of the GP models, we predict themean lifetime on a regular 100-point grid consisting 10×10 equally spaced pointsover [5K,10K]× [2cm−1,10cm−1].Figures 7.2 and 7.3 show the predicted mean lifetime based on 20-point designand 50-point design versus true mean lifetime of the equispaced grid of 100 pointswhen power parameter is set to be 1.95 and 2, respectively.All the points are roughly scattered around y = x, indicating that the fitted GPmodels work well. For p = 1.95, the GP model based on the 50-point designseems to make better predictions than the one based on the 20-point design, whilethe effect of sample size on the predictive performance is insignificant for p = 2.133lllllllllllllllllllllllllllllllllllllllllllllllllllllllllll lllllllllllllllllllllllllllllllllllllll70 110 150 19070110150190True Lifetime (ps)Predicted Lifetime  (ps)l 20 runs50 runsp = 1.95Figure 7.2: The predicted mean lifetime based on a 20-point design and a50-point design versus true mean lifetime of the equispaced grid of 100points when power parameter is set to be 1.95To quantitatively compare the prediction accuracy of the two designs, we calculatethe empirical root mean squared error (ERMSE) and the scaled root mean squarederror (SRMSE) given by134llllllllll llllllllll llllllllll llllllllllllllllllllllllllllll llllllllllllllllllll llllllllll llllllllll70 110 150 19070110150190True Lifetime (ps)Predicted Lifetime  (ps)l 20 runs50 runsp = 2Figure 7.3: The predicted mean lifetime based on a 20-point design and a50-point design versus true mean lifetime of the equispaced grid of 100points when power parameter is set to be 2.ERMSE =√1nn∑i=1(yi− yˆi)2 (7.12)SRMSE = 1ymax− ymin√1nn∑i=1(yi− yˆi)2 (7.13)135Sample Size ERMSE (p=1.95) SRMSE (p=1.95) ERMSE (p=2) SRMSE (p=2)n = 20 9.36 7.93% 8.09 6.85%n = 50 5.17 4.38% 8.37 7.09%Table 7.1: The summary of ERMSE and SRMSE values of n = 20 and n = 50The ERMSE results are similar to those evaluated by another popular measurefor prediction accuracy – the mean absolute percentage error (MAPE) defined byMAPE = 1nn∑i=1|yi− yˆi|yi(7.14)The data are summarized in Table 7.1. The prediction error decreases as samplesize n increases for p = 1.95 (SRMSE: 7.93% vs. 4.38%), but on the whole a20-point design is enough for building a good surrogate model for our classicaltrajectory model. In addition, due to randomness, different samples give differentprediction accuracy even for the sample size. This explains the slight increase inthe prediction error when the sample size increase from 20 to 50 for p = 2.For a system with only two input parameters, it is desirable to have a global sur-face to explore the relationship between input parameters and outputs. However, itgenerally requires running numerical models at a very fine grid, so it is impracticalfor those complicated systems. Instead, we can construct a predicted surface basedon a GP model fitted by a few data points. Figure 7.4 shows the global surface ofthe lifetime as a function of E and Tr obtained from the GP model based on 20-point design and p set to 1.95. Since the SRMSE of the GP model is smaller than8%, we assume the predicted surface could be an acceptable surrogate for the truesurface. In Figure 7.4, we notice that the predicted mean lifetime is monotonicallydecreasing with both rotational temperature and collision energy. The effect of therotational temperature is much weaker especially when E > 5 cm−1 and there isno strong joint interaction between Tr and E. The GP model surface can be used toevaluate thermally averaged collision lifetimes by integrating the E-dependence atgiven Tr.136Rotational Temperature Collision EnergyLifetimellllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllFigure 7.4: The surface produced by the GP model. The lines connect thevalues (circles) computed from the classical trajectories with the valuespredicted by the GP model.7.3.2 Anaysis of Interaction between Mass and PolarizabilityThe GP model can be extended to multiple collision systems for the predictionsof the collision properties of a specific molecule based on the known collisionproperties of another molecule. To illustrate this, we consider the lifetimes of thelong-lived complexes formed by benzene in collisions with Rg atoms He – Xe.As the collision system is changed from C6H6 - He to C6H6 - Xe, there are two137varying factors that determine the change of the collision dynamics: the reducedmass and the PES. These factors must be correlated but the correlations cannot beclearly determined from the scattering calculations for different systems becausethe dependence of the scattering observable on these two parameters is very differ-ent.As before, we use the GP model F(x) = β + Z(x), with x now representingthe atomic mass µA and the interaction strength De at the global minimum of theatom - molecule PES obtained by scaling the Ar - C6H6 PES. To build GP modelsto approximate the effects of mass and polarizability efficiently, we fix Tr = 4 Kand the collision energy E = 4 cm−1, and compute the collision lifetimes at two40-point maximin LH design in the interval of µA and De [4g/mol,130g/mol]×[80cm−1,520cm−1], which covers all of the Rg – C6H6 systems, and then use theGP models fitted by one design to predict the outputs of the other design. Thepredicted mean lifetime surface and contour based on one 40-points design areplotted in Figure 7.5 and Figure 7.6, which are in accordance with our previousconclusion based on the multiple linear regression model.Although in our previous studies, a multiple linear regression model with in-teraction effects x1x2 as y = β0 +β1x1 +β2x2 +β12x1x2 fits the training data withthe adjusted R2 value close to 1, it does not prove that the model is adequate forprediction. In addition, the previous training sites used for fitting the linear regres-sion model is a Cartesian product of the true atomic masses and the true potentialdepths of five rare atoms, which might not be a good design for modelling thewhole area of interest. To quantitatively compare the predictive performance of thelinear regression models and the GP models, we measure the ERMSE and SRMSEin those models. For linear models, we do the same as in the GP models, usingthe model fitted by one design to predict the outputs of the other design. The datasummarized in Table 7.2 show that the RMSE in the linear models are two timeslarger than that in the GP models, the GP models thus outperformance the linearmodels in prediction accuracy.138MassPotential depthLifetimeFigure 7.5: The 3D global response surface produced by the GP model forC6H6 – Rg collision lifetimes with the atomic mass and the PES depthplotted on the x- and y-axis at Tr = 4 K and E = 4 cm− Fitting a High-Dimensional ModelThe GP model can be applied to explore the role of the individual PES parameterson the scattering observables. To illustrate this application, we now consider thatthe mean lifetime of collision complexes as a function of E, Tr and 8 input variablesgiving the analytical form of the Rg - C6H6 PES [121]. We first create an artificialRg - C6H6 system with each potential parameter in the range between He and Xe13910509013017021025020 40 60 80 100100200300400Mass (g mol)Potential Depth (cm−1 )Figure 7.6: The contour plot produced by the GP model for C6H6 – Rg colli-sion lifetimes with the atomic mass and the PES depth plotted on the x-and y-axis at Tr = 4 K and E = 4 cm−1.Model SRMSE ERMSEGaussian Process 1 11.43 5.09%Gaussian Process 2 13.68 5.71%Linear Regression 1 27.25 12.13%Linear Regression 2 26.92 11.25%Table 7.2: The summary of ERMSE and SRMSE values of the GP modelsand linear regression modelsat rotational temperature from 4 to 10K and collision energy from 2 to 10 cm−1 butkeep the atomic mass fixed at the value of Ar, since artificial system will be reusedin the Bayesian calibration in the next section. According to the 10 times the input140dimension rule [155], we use 200 training data sampled by LH scheme to fit a GPmodel in the form as F(x) = β +Z(x). To illustrate the prediction accuracy of theGP model thus produced, we plot in Figure 7.7 the predicted mean lifetimes foranother set of 70 randomly selected data points against their computed values. InFigure 7.7, we find that most points are scattered very close to the line of y = x,and the values of ERMSE and SRMSE are 6.08 and 4%, respectively.llllllllllllllllllllllllllllllllllllllllllllllllllll llllllllllll50 100 150 20050100150200Computed Lifetime (ps)Predicted Lifetime  (ps)Figure 7.7: Accuracy of the GP model with variable PES parameters for theprediction of the collision lifetimes. The scatter plot compares the pre-dicted values with the computed values. The error of the GP model isthe deviation of the points from the diagonal line.The 10-parameter GP model contains a wealth of information on the depen-dence O(x). With the help of the GP model, we perform sensitivity analysis byusing the Functional Analysis of Variance (FANOVA) decomposition [169–171],to determine which of the input parameters have the strongest impact on the observ-141Tr E t1 t2 t3 t4 t5 t6 t7 t8contribution to total variance10−410−310−210−1100Figure 7.8: Relative effect of the variation of Tr, E and the PES parameters onthe collision lifetimes. The blue area of the bars shows the uncorrelatedcontribution of the corresponding variable and the green area – the effectthat depends on one or more other variables.able (see Figure 7.8). Large percentage of the total functional variance indicatesstrong influence on the mean lifetime. The bar plot in Figure 7.8 illustrates that theinitial collision energy plays a dominant role (78%) in the change of lifetime, andamong eight potential parameters, the location of the potential well of C-C bondfor the parallel approach of Rg to Benzene ε‖ is the most important factor (12%)in the mean lifetime. The model can also use be used to compute the uncertaintiesdue to global variation of the PES. Figure 7.9 shows the interval of the lifetimesobtained by the simultaneous ±3 % variation of all 8 PES parameters.142l l l l l l l l l l l l l l l l l l l l l l l l l l l l l l l l l l l l l l l l l2 3 4 5 6 7 8 9 1050100150200250Collision energy (cm−1)Lifetime  (ps)Figure 7.9: Energy dependence of the collision lifetime for Ar - C6H6 withthe error interval obtained by varying all the individual PES parametersby ±3 %.7.4 Bayesian calibration of dynamical calculations7.4.1 Fitting experimental dataThe GP model can be extended to fit the experimental observations by varyingthe PES, with the fit absorbing both the measurement error and the error of thedynamical calculations. To illustrate this, consider an ensemble of randomly gen-erated pseudo-experimental data and assume that the molecular dynamics calcula-tion is affected by unknown errors. Suppose experimental observations are madeover two input parameters, x and y, representing rotational temperature and ini-tial collision energy, and the observable of interest is the mean lifetime of colli-143sion complexes. The theoretical model requires eight more parameters collectivelydescribing the interaction potential surface: θ = (θ1,θ2, · · · ,θ8). The values ofθ are constant but unknown. To fabricate pseudo experimental observations, wefirst sample 50 training sites in the ranges of temperature and collision energy[4K,8K]× [4cm−1,8cm−1]. As a fact of matter, there is no need to make the ex-perimental measurements spread out evenly. In reality, we take whatever exper-imentalist could provide, and then make our numerical calculations to cover theregion of their observations. We take into account of all possible uncertainties, andexpress the relationship between experimental observations E(·) and dynamicalcalculations F(·, ·) as [138]E(x) = ρF(x,θ )+G(x)+ ε (7.15)where F(·) and G(·) are independent Gaussian random functions, with G(·) char-acterizing the inaccuracy of the scattering calculations, ρ is a scaling parameter,and ε is a normal random variable characterizing the measurement errors in exper-iments.The 50 pseudo experimental observations are generated based on 50 dynamicalcalculation results, and all the parameters about model inadequacy and measure-ment error and scaling parameter are set arbitrarily, and treated as unknowns duringthe calibration process. In Figure 7.10, the off-set of the data from the dashed di-agonal line represents the inaccuracy of the molecular dynamics calculations.We choose a 500-point Maximin LH code design over the ranges of rotationaltemperature and collision energy of the pseudo experimental observations, and therange of potential parameters [3cm−1,6cm−1]4× [3A˚, 6A˚]4. The calibration iscarried out in three steps. First, the scattering calculations are used to train the GPmodel F(·) of the scattering calculation data. In the second step, the experimentaland calculated data are used together to train the model G(·) in Eq. (7.15), usingthe parameters of F(·) from the first step averaged over possible values of t andtreating ρ and ε as variable parameters. This fixes the models F(·) and G(·) aswell as ρ and ε . Given these, Eq. (7.15) can be used to find the best-fitting PESparameters by using Markov chain Monte Carlo methods [172]. Figure 7.11 showsthe data extracted from MCMC simulations.144llllllllllllllllllll lllllllllllllllllllllllllllllllllllllllllllllllllllllllllll100 120 140 160100300500700900Dynamical calculation (ps)Observation  (ps)llllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllll700 800 900700800900Prediction (ps)Observation  (ps)Figure 7.10: Randomly generated pseudo-observation data (red circles) vsthe results obtained from the CT calculations. The off-set of the datapoints from the dashed line represents the inherent inaccuracy of thescattering calculations. Inset: Accuracy of the GP model.Generally, the purpose of calibration is to use the calibrated model for pre-dicting the real process. The problem of predicting the true process E(x) at somespecified input site x can be seen as interpolating the function E(·). The model(7.15) is also very useful when the number of the experimental data points is smallso that no meaningful dependence of the experimental data on the experimental pa-rameters can be determined. In this case, the model (7.15) provides a global surfaceof the scattering observable as a function of t with the functional dependence onthe experimental parameters determined by the molecular dynamics calculations.The inset of Figure 7.10 illustrates the accuracy of the GP model (7.15) withoutfurther optimization of the parameters t. The scatter plot of the inset in Figure 7.101450 2000 6000 10000 1400044.14.24.3Number of MCMC SimulationsSimulated Values0 1000 2000 3000 4000 50003.544.55Number of MCMC SimulationsFigure 7.11: The total 15000 MCMC simulations for θ2(left panel), last10000 simulations (out of 15000) of all eight potential parameters. Theburn-in period is 5000.has the error εS = 6.3%.7.4.2 Fitting quantum calculationsEq. (7.15) can also be used to model time-consuming quantum scattering calcu-lations with the help of much more efficient classical dynamics calculations. Toillustrate this, we consider the cross sections for rotationally inelastic C6H6 - Hecollisions [48] shown in Figure 7.12. Sixty randomly chosen quantum results areused to train the model (7.15) with ε set to zero and G(·) modelling the inaccuracyof the CT calculations. Figure 7.4 shows that the model provides an accurate en-ergy dependence of the cross sections, even near scattering resonances. The scat-tering resonances make the direct model (7.2) unstable, as shown in Figure 7.12(inset). The CT calculations in a two-function model (7.15) stabilize the model,removing the errors arising from the resonant variation of the quantum results.7.5 ConclusionIn summary, we have shown that a Gaussian Process model combined with asmall number of scattering calculations can be used to obtain an accurate multi-146lllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllll1 3 5 7 90100200300400500600Collision Energy (cm−1)Cross section  (a.u.) llllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllll1 3 5 7 90200400600Collision Energy (cm−1)Cross section  (a.u.)Figure 7.12: GP model (solid curve) of quantum scattering cross sections(squares and circles) for C6H6 - He collisions with CT calculations(not shown) serving as an intermediate model. The 60 squares arethe quantum results used to train the models. The inset shows the GPmodel (solid curve) trained by the quantum results without CT calcu-lations.dimensional dependence of the scattering observables on the experimentally con-trollable parameters and the PES parameters. In the absence of the experimentaldata, the GP model can be trained by a limited set of scattering calculations to pro-vide a multi-parameter fit of the calculation results. This can be used to integrate147the results over some of the parameters in order to evaluate efficiently thermally av-eraged quantities or to average the scattering calculations for complex moleculesover variations of the PES parameters, thereby minimizing the uncertainties due toPES inaccuracies. If experimental data are available, they can be reproduced by acombination of two independent GP models: one trained to model the scatteringcalculations and one to model the uncertainty of the scattering calculation method.The same model can be used to connect classical and quantum dynamics calcula-tions, allowing one to use classical calculations to interpolate quantum results. Themodel described here is expected to find a wide range of applications, from fittingthe interaction potentials by solving the inverse scattering problem, to analyzingthe dependence of scattering observations on external parameters, to calibratingthe accuracy of the scattering calculation methods.148Chapter 8OutlookIn this thesis, we answer several important questions pertaining to collision dy-namics of cold molecules. We first demonstrate that the collisional decoherence ofatoms or molecules can be tuned by a magnetic field, and propose a new methodfor measuring the scattering lengths of ultracold particles in excited states. In or-der to provide insights into the possibility of evaporative cooling of molecules, wetake into account the uncertainty in the potential energy surface, and perform rig-orous quantum calculations of the cross sections for diatomic molecule–diatomicmolecule collisions with multiple potential surfaces. We conclude that the major-ity of 2Σ molecules can be trapped and evaporatively cooled at millikelvin tem-peratures. We then extend our research to more complicated collision systems:atom–polyatomic molecule and polyatomic molecule–polyatomic molecule. Weuse classical dynamics to reduce the complexity of numerical calculations, andcarry out molecular dynamics simulations for rigid polyatomic molecules such asbenzene in a cold rare gas environment. We illustrate that the mean lifetimes ofthe collision complexes increase monotonically with the strength of the atom –molecule interaction, while the mean lifetimes of polyatomic molecule-moleculecollision complexes are significantly reduced due to the restricted configurationspace.During the process of solving the above problems, we realized that the ma-jor challenge of scattering calculations is the computational cost. More often thannot, the scattering calculations for complicated collisions systems are very time-149consuming. They may require several weeks of computer time for a single datapoint. Thus, within a reasonable time period, we cannot obtain enough data merelyfrom scattering calculations to obtain a good understanding of collision dynam-ics, especially for complicated systems involving multiple parameters. To over-come this challenge, inspired by statisticians dealing with computer experiments,we model an implicit function relating the collision outcomes to input factors as asample path of a stochastic process. Given a small number of scattering calcula-tions obtained from a space-filling sampling design, we built a Gaussian Processmodel to statistically approximate collision outcomes, dramatically reducing thecomputational time. To illustrate the feasibility of Gaussian Process models as anefficient surrogate for scattering calculations, we use the Benzene-Rare gas colli-sion model outlined in Chapter 6. We show that Gaussian Process models trainedby a small number of scattering calculations can provide good predictions of col-lision outcomes for a complicated system involving multiple parameters. We fur-ther extend the applications of Gaussian Process models for the collision dynamicsto the model calibration, uncertainty analysis of the potential energy surface andthe sensitivity analysis of input parameters. For example, given a few experimen-tal observations, we can use Gaussian Process models to calibrate the scatteringcalculations for improved predictions of untried experimental observations. Us-ing a Gaussian Process model as an efficient integrand, we can easily determinethe probability distribution of collision outcomes by integrating out potential pa-rameter variables V ∼ pi(·) by using simple Monte Carlo methods or numericalquadratures in smooth problems.In addition to predicting collision observables, another purpose of our theoret-ical studies of collision dynamics is to interpret the input-output relationships incollision systems. It is desirable to identify the important parameters and quantifytheir corresponding effects: the main effects of individual parameters or the jointeffects of multiple parameters on the collision outcomes. However, it is not aneasy task since those relationships are sometimes not very obvious in complicatedcollision systems. Despite their efficiency and flexibility in the interpolation ofscattering calculations, Gaussian Process models, as a class of non-parametric re-gression models, cannot provide an easily interpretable input-output relationship.In Chapter 6, we perform sensitivity analysis by using the Functional Analysis150of Variance (FANOVA) decomposition. Our future work will concentrate on theidentification of important collision parameters and the quantitative measure oftheir effects on collision outcomes. To achieve this purpose, we plan to decom-pose the complex input-output relationships y(x) ≡ y(x1,x2, · · · ,xq) in scatteringproblems into contributions from the main effects of individual parameters µi(xi)and joint interaction effect of two different parameters µi j(xi,x j) and all the otherhigher-order interaction effects involving multiple parameters [173]y(x) = µ0 +q∑i=1µi(xi)+q−1∑i=1q∑j=i+1µi j(xi,x j)+ · · ·+µ1···q(x) (8.1)whereµ0 =∫Xy(x)ω(x)dx (8.2)is an overall average of a collision outcome, andω(x) =q∏iωi(xi) (8.3)is a product of weight functions of individual parameters. The (corrected) maineffect of a given parameter xi is expressed asµi(xi) = y¯i(xi)−µ0 (8.4)wherey¯i(xi) =∫y(x)∏j 6=iω j(x j)dx j (8.5)is the marginal effect of parameter xi computed by integrating out all the otherparameters x j ′s. The (corrected) joint interaction of any two parameters xi and x jis defined byµi j(xi,x j) = y¯i j(xi,x j)− y¯i(xi)− y¯ j(x j) (8.6)151where the corresponding marginal effect of (xi,x j) isy¯i j(xi,x j) =∫y(x) ∏k 6=i, jωk(xk)dxk (8.7)Consider, for example, the artificial Rg - C6H6 system previously discussed inSection 7.3.3. This system has varying potential parameters ranging from He to Xe,and fixed atomic mass corresponding to the mass of Ar. The collision outcomes ofthe artificial system y(x), such as the mean lifetime of the collision complexes andthe scattering cross sections, are functions of ten parameters, including the collisionenergy, the rotational temperature and eight potential parameters. Using Eq. (8.1),we can write y(x) as a sum of 210 terms, including one total overall average µ0,ten (corrected) main effects of each parameters µi(xi) for i = 1,2, · · · ,10, forty-five joint interaction effects µi j(xi,x j), and a large number of higher-order effects.However, it is completely unnecessary to compute all the effects. The sensitivityanalysis done in Section 7.3.3 has already screened out the important effects, withthe main effects being the collision energy µ(Ecol) and the location of the potentialwell of C–C bond for the parallel approach of the atom to benzene µ(ε‖). UsingEqs. (8.2-8.4), we can not only quantify those effects but also visualize them todetermine the types of those effects, linear or non-linear.At this point, y(·) is still an unknown function. Given the estimator of anunknown output y(x0) based on a Gaussian Process model in Eq. (7.10), we cancompute the estimated marginal effect by carrying out the required integration inEq. (7.10) for the important effects and the corresponding corrected effects bysubtracting all lower-order effects from the estimated marginal effects.In conclusion, the introduction of Gaussian Process models to our current the-oretical studies in collision dynamics of cold molecules offers us a wealth of in-formation on collision systems of interest with only a small number of classical orquantum dynamics calculations.152Bibliography[1] S. Chu, L. Hollberg, J. E. Bjorkholm, A. Cable, and A. Ashkin.Three-dimensional viscous confinement and cooling of atoms by resonanceradiation pressure. Phys. Rev. Lett., 55(1):48, 1985. → pages 1[2] M. H. Anderson, J. R. Ensher, M. R. Matthews, C. E. Wieman, and E. A.Cornell. Observation of bose-einstein condensation in a dilute atomicvapor. Science, 269(5221):198, 1995. → pages 1[3] K. B. Davis, M. O. Mewes, M. R. Andrews, N. J. Van Druten, D. S. Durfee,D. M. Kurn, and W. Ketterle. Bose-einstein condensation in a gas ofsodium atoms. Phys. Rev. Lett., 75(22):3969, 1995. → pages 1[4] B. DeMarco and D. S. Jin. Onset of fermi degeneracy in a trapped atomicgas. Science, 285(5434):1703, 1999. → pages 1[5] C. A. Regal, M. Greiner, and D. S. Jin. Observation of resonancecondensation of fermionic atom pairs. Phys. Rev. Lett., 92(4):040403,2004. → pages 1[6] M. R. Andrews, C. G. Townsend, H.-J. Miesner, D. S. Durfee, D. M. Kurn,and W. Ketterle. Observation of interference between two bosecondensates. Science, 275(5300):637, 1997. → pages 1[7] A. G. Truscott, K. E. Strecker, W. I. McAlexander, G. B. Partridge, andR. G. Hulet. Observation of fermi pressure in a gas of trapped atoms.Science, 291(5513):2570, 2001. → pages 1[8] M. Cadoret, E. de Mirandes, P. Clade´, S. Guellati-Khe´lifa, C. Schwob,F. Nez, L. Julien, and F. Biraben. Combination of bloch oscillations with aramsey-borde´ interferometer: New determination of the fine structureconstant. Phys. Rev. Lett., 101(23):230801, 2008. → pages 1153[9] D. DeMille. Quantum computation with trapped polar molecules. Phys.Rev. Lett., 88(6):067901, 2002. → pages 1[10] S. Ospelkaus, K.-K. Ni, D. Wang, M. H. G. De Miranda, B. Neyenhuis,G. Que´me´ner, P. S. Julienne, J. L. Bohn, D. S. Jin, and J. Ye.Quantum-state controlled chemical reactions of ultracoldpotassium-rubidium molecules. Science, 327(5967):853, 2010. → pages 1[11] E. S. Shuman, J. F. Barry, and D. DeMille. Laser cooling of a diatomicmolecule. Nature, 467(7317):820, 2010. → pages 1[12] J. Doyle, B. Friedrich, R. V. Krems, and F. Masnou-Seeuws. Editorial: QuoVadis, cold molecules? Eur. Phys. J. D, 31(2):149, 2004. → pages 2[13] L. Carr, D. DeMille, R. V. Krems, and J. Ye. Cold and ultracold molecules:science, technology and applications. New J. Phys., 11(5):055049, 2009.→ pages 2, 22, 31, 39, 40, 58[14] J. T. Bahns, P. L. Gould, and W. C. Stwalley. Formation of cold (T < 1K)molecules. Adv. At. Mol. Opt. Phys., 42:171, 2000. → pages 2[15] S. Inouye, M. R. Andrews, J. Stenger, H.-J. Miesner, D. M. Stamper-Kurn,and W. Ketterle. Observation of feshbach resonances in a bose–einsteincondensate. Nature, 392(6672):151, 1998. → pages 2[16] J. Werner, A. Griesmaier, S. Hensler, J. Stuhler, T. Pfau, A. Simoni, andE. Tiesinga. Observation of feshbach resonances in an ultracold gas of52Cr. Phys. Rev. Lett., 94(18):183201, 2005. → pages 2[17] K. Xu, T. Mukaiyama, J. R. Abo-Shaeer, J. K. Chin, D. E. Miller, andW. Ketterle. Formation of quantum-degenerate sodium molecules. Phys.Rev. Lett., 91(21):210402, 2003. → pages 2[18] A. Fioretti, D. Comparat, A. Crubellier, O. Dulieu, F. Masnou-Seeuws, andP. Pillet. Formation of cold Cs2 molecules through photoassociation. Phys.Rev. Lett., 80(20):4402, 1998. → pages 2[19] H. R. Thorsheim, J. Weiner, and P. S. Julienne. Laser-inducedphotoassociation of ultracold sodium atoms. Phys. Rev. Lett., 58(23):2420,1987. → pages 2[20] F. Masnou-Seeuws and P. Pillet. Formation of ultracold molecules (T≤200µK) via photoassociation in a gas of laser-cooled atoms. Adv. At. Mol. Opt.Phys., 47:53, 2001. → pages 2, 40154[21] K.-K. Ni, S. Ospelkaus, D. Wang, G. Que´me´ner, B. Neyenhuis, M. H. G.De Miranda, J. L. Bohn, J. Ye, and D. S. Jin. Dipolar collisions of polarmolecules in the quantum regime. Nature, 464(7293):1324, 2010. → pages2[22] J. Weinstein, R. deCarvalho, T. Guillet, B. Friedrich, and J. M. Doyle.Magnetic trapping of calcium monohydride molecules at millikelvintemperatures. Nature, 395:148, 1998. → pages 2, 39, 40, 42, 59[23] P. Bowe, L. Hornekær, C. Brodersen, M. Drewsen, J. S. Hangst, and J. P.Schiffer. Sympathetic crystallization of trapped ions. Phys. Rev. Lett.,82(10):2071, 1999. → pages 2[24] W. Ketterle and N. J. Van Druten. Evaporative cooling of trapped atoms.Adv. At. Mol. Opt. Phys., 37:181, 1996. → pages 2, 40[25] K. Maussang, D. Egorov, J. S. Helton, S. V. Nguyen, and J. M. Doyle.Zeeman relaxation of CaF in low-temperature collisions with Helium.Phys. Rev. Lett., 94(12):123002, 2005. → pages 2, 40, 52[26] W. C. Campbell, E. Tsikata, H. Lu, L. D. van Buuren, and J. M. Doyle.Magnetic trapping and zeeman relaxation of NH (X 3Σ−). Phys. Rev. Lett.,98(21):213001, 2007. → pages 2, 31, 40[27] E. Tsikata, W. C. Campbell, M. T. Hummon, H.-I. Lu, and J. M. Doyle.Magnetic trapping of NH molecules with 20 s lifetimes. New J. Phys.,12(6):065028, 2010. → pages 2, 40[28] D. Egorov, J. D. Weinstein, D. Patterson, B. Friedrich, and J. M. Doyle.Spectroscopy of laser-ablated buffer-gas-cooled PbO at 4 K and theprospects for measuring the electric dipole moment of the electron. Phys.Rev. A, 63(3):030501, 2001. → pages 2[29] M. Stoll, J. M. Bakker, T. C. Steimle, G. Meijer, and A. Peters. Cryogenicbuffer-gas loading and magnetic trapping of CrH and MnH molecules.Phys. Rev. A, 78(3):032707, 2008. → pages 2[30] D. Patterson, E. Tsikata, and J. M. Doyle. Cooling and collisions of largegas phase molecules. Phys. Chem. Chem. Phys., 12(33):9736, 2010. →pages 2, 59, 60, 75, 129[31] D. Patterson and J. M. Doyle. A slow, continuous beam of coldbenzonitrile. Phys. Chem. Chem. Phys., 17(7):5372, 2015. → pages 2, 129155[32] T. Hasegawa and T. Shimizu. Resonant oscillation modes ofsympathetically cooled ions in a radio-frequency trap. Phys. Rev. A,66(6):063404, 2002. → pages 2[33] K. Mølhave and M. Drewsen. Formation of translationally cold MgH+ andMgD+ molecules in an ion trap. Phys. Rev. A, 62(1):011401, 2000. →pages 2[34] T. Baba and I. Waki. Chemical reaction of sympathetically laser-cooledmolecular ions. J. Chem. Phys., 116(5):1858, 2002. → pages 2[35] Matthew T. Hummon, Wesley C. Campbell, Hsin-I Lu, Edem Tsikata,Yihua Wang, and John M. Doyle. Magnetic trapping of atomic nitrogen(14N) and cotrapping of NH (X 3Σ−). Phys. Rev. A, 78(5):050702, 2008. →pages 2[36] M. Lara, J. L. Bohn, D. Potter, P. Solda´n, and J. M. Hutson. UltracoldRb-OH collisions and prospects for sympathetic cooling. Phys. Rev. Lett.,97(18):183201, 2006. → pages 2, 41[37] A. O. G. Wallis and J. M. Hutson. Production of ultracold NH moleculesby sympathetic cooling with Mg. Phys. Rev. Lett., 103(18):183201, 2009.→ pages 2, 6, 31, 41[38] P. Solda´n, P. S. Z˙uchowski, and J. M. Hutson. Prospects for sympatheticcooling of polar molecules: NH with alkali - metal and alkaline - earthatoms-a new hope. Faraday Disc., 142:191, 2009. → pages 2, 31[39] M. L. Gonza´lez-Martı´nez and J. M. Hutson. Effect of hyperfineinteractions on ultracold molecular collisions: NH (3Σ−) with Mg (1S) inmagnetic fields. Phys. Rev. A, 84(5):052706, 2011. → pages 2[40] B. K. Stuhl, M. T. Hummon, M. Yeo, G. Que´me´ner, J. L. Bohn, and J. Ye.Evaporative cooling of the dipolar hydroxyl radical. Nature,492(7429):396, 2012. → pages 3, 41[41] R. V. Krems. Molecules near absolute zero and external field control ofatomic and molecular dynamics. Int. Rev. Phys. Chem., 24(1):99, 2005. →pages 3[42] R. V. Krems and A. Dalgarno. Quantum-mechanical theory ofatom-molecule and molecular collisions in a magnetic field: Spindepolarization. J. Chem. Phys., 120(5):2296, 2004. → pages 4, 5, 18, 29,32, 40, 41, 42, 43, 45, 47, 50, 52, 56156[43] R. V. Krems, A. Dalgarno, N. Balakrishnan, and G. C. Groenenboom.Spin-flipping transitions in 2Σ molecules induced by collisions withstructureless atoms. Phys. Rev. A, 67(6):060703, 2003. → pages 4, 40, 50[44] H. Cybulski, R. V. Krems, H. R. Sadeghpour, A. Dalgarno, J. Kłos, G. C.Groenenboom, A. Van Der Avoird, D. Zgid, and G. Chałasin´ski.Interaction of NH (X 3Σ−) with He: Potential energy surface, bound states,and collisional Zeeman relaxation. J. Chem. Phys., 122(9):094307, 2005.→ pages 4, 6, 40[45] J. Ku¨pper, F. Filsinger, and G. Meijer. Manipulating the motion of largeneutral molecules. Faraday Disc., 142:155, 2009. → pages 4, 59[46] Ch. Daussy, T. Marrel, A. Amy-Klein, C. T. Nguyen, Ch. J. Borde´, and Ch.Chardonnet. Limit on the parity nonconserving energy difference betweenthe enantiomers of a chiral molecule by laser spectroscopy. Phys. Rev.Lett., 83(8):1554, 1999. → pages 4, 59[47] Z. Li and E. J. Heller. Cold collisions of complex polyatomic molecules. J.Chem. Phys., 136(5):054306, 2012. → pages 4, 60, 74[48] Z. Li, R. V. Krems, and E. J. Heller. Collision dynamics of polyatomicmolecules containing carbon rings at low temperatures. J. Chem. Phys.,141(10):104317, 2014. → pages 4, 60, 128, 132, 146[49] L. M. C. Janssen, A. van der Avoird, and G. C. Groenenboom. On the roleof the magnetic dipolar interaction in cold and ultracold collisions:numerical and analytical results for NH (3Σ−) + NH (3Σ−). Eur. Phys. J. D,65(1):177, 2011. → pages 5, 6, 41[50] L. M. C. Janssen, P. S. Z˙uchowski, A. van der Avoird, G. C. Groenenboom,and J. M. Hutson. Cold and ultracold NH-NH collisions in magnetic fields.Phys. Rev. A, 83(2):022713, 2011. → pages 5, 6, 41, 45[51] L. M. C. Janssen, P. S. Zuchowski, A. van der Avoird, J. M. Hutson, andG. C. Groenenboom. Cold and ultracold NH-NH collisions: the field-freecase. J. Chem. Phys., 134(12):124309, 2011. → pages 5, 6, 45[52] J. R. Taylor. Scattering theory: the quantum theory of nonrelativisticcollisions. Courier Corporation, 2012. → pages 10, 25[53] A. M. Arthurs and A. Dalgarno. The theory of scattering by a rigid rotator.Proc. R. Soc. London A, 256(1287):540, 1960. → pages 10, 43157[54] R. V. Krems, W. C. Stwalley, and B Friedrich. Cold molecules: theory,experiment, applications. CRC press, 2009. → pages 10, 16, 22, 39[55] D. J. Griffiths. Introduction to quantum mechanics. Addison-Wesley, 2005.→ pages 11, 13, 14[56] W. A. Lester. Dynamics of molecular collisions. Plenum, New York, 1976.→ pages 12[57] B. R. Johnson. The multichannel log derivative method for scatteringcalculations. J. Comput. Phys., 13(3):445, 1973. → pages 19, 45[58] D. E. Manolopoulos, M. J. Jamieson, and A. D. Pradhan. Johnson’s logderivative algorithm rederived. J. Comput. Phys., 105(1):169, 1993. →pages 19, 45[59] M. Bu¨chner G. Tre´nec A. Miffre, M. Jacquey and J. Vigue´. Atominterferometry. Phys. Scr., 74(2):C15, 2006. → pages 21[60] M. A. Nielsen and I. L. Chuang. Quantum computation and quantuminformation. Cambridge university press, 2010. → pages 21[61] G. K. Brennen A. Micheli and P. Z˙oller. A toolbox for lattice-spin modelswith polar molecules. Nature Physics, 2(5):341, 2006. → pages 21, 40[62] S. Mukamel. Principles of nonlinear optical spectroscopy. OxfordUniversity, New York, 1995. → pages 22[63] M. Shapiro and P. Brumer. Quantum Control of Molecular Processes.Wiley, 2012. → pages 22[64] D. Budker, W. Gawlik, D. F. Kimball, S. M. Rochester, V. V. Yashchuk, andA. Weis. Resonant nonlinear magneto-optical effects in atoms. Rev. Mod.Phys., 74(4):1153, 2002. → pages 22[65] H. Stapelfeldt and T. Seideman. Colloquium: Aligning molecules withstrong laser pulses. Rev. Mod. Phys., 75(2):543, 2003. → pages 22[66] A. J. Orr-Ewing and R. N. Zare. Orientation and alignment of reactionproducts. Annu. Rev. Phys. Chem., 45(1):315, 1994. → pages 22[67] E. Joos, H. D. Zeh, C. Kiefer, D. Giulini, J. Kupsch, and I.-O. Stamatescu.Decoherence and the appearance of a classical world in quantum theory.Springer Science & Business Media, 2003. → pages 22158[68] S. Ramakrishna and Tamar Seideman. Intense laser alignment indissipative media as a route to solvent dynamics. Phys. Rev. Lett.,95(11):113001, 2005. → pages 22[69] C. Reinhold, J. Burgdo¨rfer, and F. Dunning. Collisional decoherence invery-high-n rydberg atoms. Methods Phys. Res. B, 233(1):48, 2005. →pages 22[70] C. Chin, R. Grimm, P. Julienne, and E. Tiesinga. Feshbach resonances inultracold gases. Rev. Mod. Phys., 82(2):1225, 2010. → pages 22, 23, 58,128[71] K. Blum. Density Matrix Theory and Applications. Springer Science &Business Media, 1996. → pages 23[72] C. J. Hemming and R. V. Krems. Collisional decoherence of internal-statesuperpositions in a trapped ultracold gas. Phys. Rev. A, 81(5):052701,2010. → pages 24, 25, 26, 32[73] N. Balakrishnan, V. Kharchenko, R. C. Forrey, and A. Dalgarno. Complexscattering lengths in multi-channel atom-molecule collisions. Chem. Phys.Lett., 280(1):5, 1997. → pages 27[74] R. N. Zare. Angular Momentum: Understanding spatial aspects inchemistry and physics. Wiley, New York, 1988. → pages 31, 44[75] A. J. Moerdijk, B. J. Verhaar, and A. Axelsson. Resonances in ultracoldcollisions of 6Li, 7Li, and 23Na. Phys. Rev. A, 51(6):4852, 1995. → pages36[76] M. L. Gonza´lez-Martı´nez and J. M. Hutson. Ultracold atom-moleculecollisions and bound states in magnetic fields: Tuning zero-energyFeshbach resonances in He-NH (X 3Σ−). Phys. Rev. A, 75(2):022702,2007. → pages 36, 40, 41[77] H. L. Bethlem, G. Berden, and G. Meijer. Decelerating neutral dipolarmolecules. Phys. Rev. Lett., 83(8):1558, 1999. → pages 39[78] K. M. Jones, E. Tiesinga, P. D. Lett, and P. S. Julienne. Ultracoldphotoassociation spectroscopy: Long-range molecules and atomicscattering. Rev. Mod. Phys., 78(2):483, 2006. → pages 40[79] J. J. Hudson, B. E. Sauer, M. R. Tarbutt, and E. A. Hinds. Measurement ofthe electron electric dipole moment using YbF molecules. Phys. Rev. Lett.,89(2):023003, 2002. → pages 40159[80] R. V. Krems. Cold controlled chemistry. Phys. Chem. Chem. Phys.,10(28):4079, 2008. → pages 40, 50, 58[81] M. Lemeshko, R. V. Krems, J. M. Doyle, and S. Kais. Manipulation ofmolecules with electromagnetic fields. Mol. Phys., 111(12-13):1648, 2013.→ pages 40[82] M. Stoll, J. M. Bakker, T. C. Steimle, G. Meijer, and A. Peters. Cryogenicbuffer-gas loading and magnetic trapping of CrH and MnH molecules.Phys. Rev. A, 78(3):032707, 2008. → pages 40[83] A. Volpi and J. L. Bohn. Magnetic-field effects in ultracold molecularcollisions. Phys. Rev. A, 65(5):052712, 2002. → pages 40, 41[84] W. C. Campbell, T. V. Tscherbul, H.-I. Lu, E. Tsikata, R. V. Krems, andJ. M. Doyle. Mechanism of collisional spin relaxation in 3Σ molecules.Phys. Rev. Lett., 102(1):013003, 2009. → pages 40[85] M. Tacconi, E. Bodo, and F. A. Gianturco. Sympathetic cooling of NH(3Σ−) molecules by Rb and Cs atoms at ultralow energies. Phys. Rev. A,75(1):012708, 2007. → pages 40, 41[86] X. T. Xu, X. Shao, C. H. Yu, C. Y. Sun, W. Huang, and E. Feng. Cold andultracold collisions of MgH (2Σ+) with Helium. Eur. Phys. J. D, 65(3):383,2011. → pages 40[87] E. Feng, X. Shao, C. Yu, C. Sun, and W. Huang. Low energy collisions ofCN (x 2Σ+) with He in magnetic fields. J. Chem. Phys., 136(5):054302,2012. → pages 40[88] R. V. Krems and A. Dalgarno. Disalignment transitions in cold collisionsof 3P atoms with structureless targets in a magnetic field. Phys. Rev. A,68(1):013406, 2003. → pages 41[89] T. V. Tscherbul, Y. V. Suleimanov, V. Aquilanti, and R. V. Krems.Magnetic field modification of ultracold molecule-molecule collisions.New J. Phys., 11(5):055021, 2009. → pages 41[90] Y. V. Suleimanov, T. V. Tscherbul, and R. V. Krems. Efficient method forquantum calculations of molecule-molecule scattering properties in amagnetic field. J. Chem. Phys., 137(2):024103, 2012. → pages 41, 42, 43,44, 128160[91] T. V. Tscherbul and A. Dalgarno. Quantum theory of molecular collisionsin a magnetic field: Efficient calculations based on the total angularmomentum representation. J. Chem. Phys., 133(18):184104, 2010. →pages 42, 43[92] T. V. Tscherbul, T. A. Grinev, H.-G. Yu, A. Dalgarno, J. Kłos, L. Ma, andM. H. Alexander. Cold collisions of polyatomic molecular radicals withs-state atoms in a magnetic field: An ab initio study of He + CH2 collisions(X˜). J. Chem. Phys., 137(10):104302, 2012. → pages 43[93] J. Pe´rez-Rı´os, J. Campos-Martı´nez, and M. I. Herna´ndez. Ultracold O2+O2 collisions in a magnetic field: On the role of the potential energysurface. J. Chem. Phys., 134(12):124310, 2011. → pages 43, 50[94] P. S. Z˙uchowski and J. M. Hutson. Cold collisions of N (4S) atoms and NH(3Σ−) molecules in magnetic fields. Phys. Chem. Chem. Phys., 13(9):3669,2011. → pages 43[95] A. O. G. Wallis, E. J. J. Longdon, P. S. Z˙uchowski, and J. M. Hutson. Theprospects of sympathetic cooling of NH molecules with Li atoms. Eur.Phys. J. D, 65(1):151, 2011. → pages 43[96] T. V. Tscherbul. Total-angular-momentum representation foratom-molecule collisions in electric fields. Phys. Rev. A, 85(5):052710,2012. → pages 43[97] T. V. Tscherbul, Gˆ. Barinovs, J. Kłos, and R. V. Krems. Formation of slowmolecules in chemical reactions in crossed molecular beams. Phys. Rev. A,78(2):022705, 2008. → pages 45[98] L. M. C. Janssen, A. van der Avoird, and G. C. Groenenboom. Quantumreactive scattering of ultracold NH (X 3Σ−) radicals in a magnetic trap.Phys. Rev. Lett., 110(6):063201, 2013. → pages 50[99] T. V. Tscherbul, J. Kł os, and A. A. Buchachenko. Ultracold spin-polarizedmixtures of 2Σ molecules with S-state atoms: Collisional stability andimplications for sympathetic cooling. Phys. Rev. A, 84(5):040701, 2011. →pages 52[100] D. M. Egorov. Buffer-gas cooling of diatomic molecules. PhD thesis,Harvard University Cambridge, Massachusetts, 2004. → pages ix, 52, 56,57161[101] T. V. Tscherbul and R. V. Krems. Controlling electronic spin relaxation ofcold molecules with electric fields. Phys. Rev. Lett., 97(8):083201, 2006.→ pages 52[102] M. Mizushima. The Theory of Rotating Diatomic Molecules. Wiley, NewYork, 1975. → pages 56[103] H. R. Sadeghpour, J. L. Bohn, M. J. Cavagnero, B. D. Esry, I. I. Fabrikant,J. H. Macek, and A. R. P. Rau. Collisions near threshold in atomic andmolecular physics. J. Phys. B: At. Mol. Opt. Phys., 33(5):R93, 2000. →pages 58[104] N Balakrishnan and A Dalgarno. Chemistry at ultracold temperatures.Chem. Phys. Lett., 341(5):652, 2001. → pages 58[105] M. T. Bell and T. P. Softley. Ultracold molecules and ultracold chemistry.Mol. Phys., 107(2):99, 2009. → pages 58[106] M. A. Baranov, M. Dalmonte, G. Pupillo, and P. Zoller. Condensed mattertheory of dipolar quantum gases. Chem. Rev., 112(9):5012, 2012. → pages58[107] D. J. Nesbitt. Toward state-to-state dynamics in ultracold collisions:Lessons from high-resolution spectroscopy of weakly bound molecularcomplexes. Chem. Rev., 112(9):5062, 2012. → pages 58[108] D. S. Jin and J. Ye. Polar molecules in the quantum regime. Phys. Today,64(5):27, 2011. → pages 58[109] K.-K. Ni, S. Ospelkaus, D. J. Nesbitt, J. Ye, and D. S. Jin. A dipolar gas ofultracold molecules. Phys. Chem. Chem. Phys., 11(42):9626, 2009. →pages 58[110] S. Ospelkaus, K.-K. Ni, M. H. G. de Miranda, B. Neyenhuis, D Wang,S. Kotochigova, P. S. Julienne, D. S. Jin, and Ye. J. Ultracold polarmolecules near quantum degeneracy. Faraday Disc., 142:351, 2009. →pages 58[111] K.-K. Ni, S. Ospelkaus, M. H. G. de Miranda, A. Pe’er, B. Neyenhuis, J. J.Zirbel, S Kotochigova, P. S. Julienne, D. S. Jin, and J. Ye. A highphase-space-density gas of polar molecules. Science, 322(5899):231, 2008.→ pages 58162[112] D. Egorov, W. C. Campbell, B. Friedrich, S. E. Maxwell, E. Tsikata, L. D.van Buuren, and J. M. Doyle. Buffer-gas cooling of NH via the beamloaded buffer-gas method. Eur. Phys. J. D, 31(2):307, 2004. → pages 59[113] P. Blythe, B. Roth, U. Fro¨hlich, H. Wenz, and S. Schiller. Production ofultracold trapped molecular hydrogen ions. Phys. Rev. Lett.,95(18):183002, 2005. → pages 59[114] K. Wohlfart, F. Gra¨tz, F. Filsinger, H. Haak, G. Meijer, and J. Ku¨pper.Alternating-gradient focusing and deceleration of large molecules. Phys.Rev. A, 77(3):031404, 2008. → pages 59[115] A. Mody, M. Haggerty, J. M. Doyle, and E. J. Heller. No-sticking effectand quantum reflection in ultracold collisions. Phys. Rev. B, 64(8):085418,2001. → pages 59[116] M. Mayle, B. P. Ruzic, and J. L. Bohn. Statistical aspects of ultracoldresonant scattering. Phys. Rev. A, 85(6):062712, 2012. → pages 59[117] J. L. Bohn, A. V. Avdeenkov, and M. P. Deskevich. Rotational Feshbachresonances in ultracold molecular collisions. Phys. Rev. Lett.,89(20):203202, 2002. → pages 59[118] M. Mayle, G. Quemener, B. P. Ruzic, and J. L. Bohn. Scattering ofultracold molecules in the highly resonant regime. Phys. Rev. A,87(1):012709, 2013. → pages 59[119] J. F. E. Croft and J. L. Bohn. Long-lived complexes and chaos in ultracoldmolecular collisions. Phys. Rev. A, 89(1):012714, 2014. → pages 59[120] D. C. Rapaport. The art of molecular dynamics simulation. Cambridgeuniversity press, 2004. → pages 61[121] F. Pirani, M. Albertı´, A. Castro, M. Moix Teixidor, and D. Cappelletti.Atom-bond pairwise additive representation for intermolecular potentialenergy surfaces. Chem. Phys. Lett., 394(1):37, 2004. → pages 66, 128,132, 139[122] T. Tran-Duc. Mathematical modelling for three problems innanotechnology. PhD thesis, University of Wollongong, 2011. → pages 67[123] M. O. Sinnokrot and C. D. Sherrill. Highly accurate coupled clusterpotential energy curves for the benzene dimer: sandwich, T-shaped, andparallel-displaced configurations. J. Phys. Chem. A, 108(46):10200, 2004.→ pages 73163[124] C. D. Sherrill, T. Takatani, and E. G. Hohenstein. An assessment oftheoretical methods for nonbonded interactions: Comparison to completebasis set limit coupled-cluster potential energy curves for the benzenedimer, the methane dimer, benzene- methane, and benzene-H2S. J. Phys.Chem. A, 113(38):10146, 2009. → pages 73[125] J. G. Hill, J. A. Platts, and H.-J. Werner. Calculation of intermolecularinteractions in the benzene dimer using coupled-cluster and local electroncorrelation methods. Phys. Chem. Chem. Phys., 8(35):4072, 2006. →pages 73[126] A. Van Der Avoird, R. Podeszwa, K. Szalewicz, C. Leforestier, R. vanHarrevelt, P. R. Bunker, M. Schnell, G. von Helden, and G. Meijer.Vibration-rotation-tunneling states of the benzene dimer: An ab initiostudy. Phys. Chem. Chem. Phys., 12(29):8219, 2010. → pages 73[127] D. Cappelletti, M. Bartolomei, F. Pirani, and V. Aquilanti. Molecular beamscattering experiments on benzene-rare gas systems: Probing the potentialenergy surfaces for the C6H6-He, -Ne, and -Ar dimers. J Phys. Chem. A,106(45):10764, 2002. → pages 84[128] M. Schnell, U. Erlekam, P. R. Bunker, G. von Helden, J.-U. Grabow,G. Meijer, and A. Van Der Avoird. Structure of the benzenedimer-governed by dynamics. Angew. Chem. Int. Ed., 52(19):5180, 2013.→ pages 84[129] D. Patterson. Buffer gas cooled beams and cold molecular collisions. PhDthesis, Harvard University, 2010. → pages 89[130] J. Cui and R. V. Krems. Collisional decoherence of superposition states inan ultracold gas near a Feshbach resonance. Phys. Rev. A, 86(2):022703,2012. → pages 90[131] J. Cui and R. V. Krems. Elastic and inelastic collisions of 2Σ molecules in amagnetic field. Phys. Rev. A, 88(4):042705, 2013. → pages 90[132] J. Cui, Z. Li, and R. V. Krems. Collision lifetimes of polyatomic moleculesat low temperatures: Benzene-benzene vs benzene-rare gas atom collisions.J. Chem. Phys., 141(16):164315, 2014. → pages 90, 132[133] M. D. Frye and J. M. Hutson. Collision cross sections for thethermalization of cold gases. Phys. Rev. A, 89(5):052705, 2014. → pages91, 128164[134] M. L. Gonza´lez-Martı´nez and J. M. Hutson. Ultracold hydrogen atoms: aversatile coolant to produce ultracold molecules. Phys. Rev. Lett.,111(20):203004, 2013. → pages 91, 128[135] J. Sacks, W. J. Welch, T. J. Mitchell, and H. P. Wynn. Design and analysisof computer experiments. STAT. SCI., 4(4):409, 1989. → pages 91, 105,107, 128, 131[136] T. J. Santner, B. J. Williams, and W. I. Notz. The design and analysis ofcomputer experiments. Springer Science & Business Media, 2003. →pages 91, 102, 105, 107[137] C. E. Rasmussen. Gaussian processes in machine learning. In Advancedlectures on machine learning, page 63. Springer, 2004. → pages 91, 105,128[138] M. Kennedy and A. O’Hagan. Bayesian calibration of computer models. J.Roy. Statist. Soc. Ser. B, 63(3):425, 2001. → pages 91, 105, 144[139] M. Kennedy and A. O’Hagan. Supplementary details on bayesiancalibration of computer models. Technical report, Internal Report. URLhttps://stat.duke.edu/courses/Spring14/sta961.01/ref/KennOHag2000b.pdf,2001. → pages 91[140] B. Flury. A first course in multivariate statistics. Springer Science &Business Media, 1997. → pages 92[141] R. J. Adler. The geometry of random fields, volume 62. SIAM, 1981. →pages 93, 96, 130[142] H. Crame´r and M. R. Leadbetter. Stationary and related stochasticprocesses: Sample function properties and their applications. CourierCorporation, 2013. → pages 96, 130[143] T. Mitchell, M. Morris, and D. Ylvisaker. Existence of smoothed stationaryprocesses on an interval. Stoch. Proc. Appl., 35(1):109, 1990. → pages 97,130[144] N. A. Cressie. Statistics for spatial data. John Wiley & Sons, New York,2015. → pages 97, 102, 130[145] M. L. Stein. Interpolation of spatial data: some theory for kriging.Springer Science & Business Media, 1999. → pages 97, 102, 130165[146] M. Abt. Estimating the prediction mean squared error in Gaussianstochastic processes with exponential correlation structure. Scand. J. Stat.,26(4):563, 1999. → pages 97, 130[147] M. D. McKay, R. J. Beckman, and W. J. Conover. Comparison of threemethods for selecting values of input variables in the analysis of outputfrom a computer code. Technometrics, 21(2):239, 1979. → pages 99[148] M. Stein. Large sample properties of simulations using Latin hypercubesampling. Technometrics, 29(2):143, 1987. → pages 99[149] R. L. Iman. Latin hypercube sampling. In Encyclopedia of QuantitativeRisk Analysis and Assessment. John Wiley & Sons, 2008. → pages 99[150] R Carnell. lhs: Latin hypercube samples. R package version 0.5, 2009. →pages 100[151] M. E. Johnson, L. M. Moore, and D. Ylvisaker. Minimax and maximindistance designs. J. Stat. Plan. Inference, 26(2):131, 1990. → pages 100[152] G. Casella and R. L. Berger. Statistical inference. Duxbury Press, 2001. →pages 100[153] T. W. Simpson, J. D. Poplinski, P. N. Koch, and J. K. Allen. Metamodelsfor computer-based engineering design: survey and recommendations.Eng. Comput., 17(2):129, 2001. → pages 102[154] V. Dubourg, B. Sudret, and F. Deheeger. Metamodel-based importancesampling for structural reliability analysis. Probabilist. Eng. Mech., 33:47,2013. → pages 102[155] J. L. Loeppky, J. Sacks, and W. J. Welch. Choosing the sample size of acomputer experiment: A practical guide. Technometrics, 51(4):366, 2009.→ pages 102, 132, 141[156] J. E. Oakley and A. O’Hagan. Probabilistic sensitivity analysis of complexmodels: a Bayesian approach. J. Roy. Statist. Soc. Ser. B, 66(3):751, 2004.→ pages 105[157] D. Higdon, M. Kennedy, J. C. Cavendish, J. A. Cafeo, and R. D. Ryne.Combining field data and computer simulations for calibration andprediction. SIAM J. Sci. Comput., 26(2):448, 2004. → pages 105, 107, 118,128166[158] D. Higdon, J. Gattiker, B. Williams, and M. Rightley. Computer modelcalibration using high-dimensional output. JASA, 103(482):570, 2008. →pages 107, 118, 128[159] R. B. Gramacy. Bayesian treed Gaussian process models. PhD thesis,UNIVERSITY OF CALIFORNIA SANTA CRUZ, 2005. → pages 107[160] R. B. Gramacy and H. K. Lee. Bayesian treed gaussian process modelswith an application to computer modeling. JASA, 103(483):1119. → pages107[161] R. B. Gramacy and H. K. H. Lee. Cases for the nugget in modelingcomputer experiments. Stat. Comput., 22(3):713, 2012. → pages 110[162] S. Shi and H. Rabitz. Sensitivity analysis in molecular dynamics andinverse scattering. Comp. Phys. Rep., 10(1):3, 1989. → pages 128[163] R. Baer and R. Kosloff. Inversion of ultrafast pump-probe spectroscopicdata. J. Phys. Chem., 99(9):2534, 1995. → pages 128[164] R. B. Bernstein. Atom-molecule collision theory: a guide for theexperimentalist. Springer Science & Business Media, 2013. → pages 128[165] N. F. Mott and H. S. W. Massey. The theory of atomic collisions,volume 35. Clarendon Press Oxford, 1965. → pages 128[166] J. O. Ramsay. Functional data analysis. In Encyclopedia of StatisticalSciences. John Wiley & Sons, 2004. → pages 128[167] K. Liu, R. T. Skodje, and D. E. Manolopoulos. Resonances in bimolecularreactions. Phys. Chem. Comm., 5(4):27, 2002. → pages 128[168] D. Patterson, M. Schnell, and J. M. Doyle. Enantiomer-specific detectionof chiral molecules via microwave spectroscopy. Nature, 497(7450):475,2013. → pages 129[169] G. Pujol. Sensitivity: Sensitivity analysis. r package version 1.4–0, 2008.→ pages 141[170] A. Saltelli, M. Ratto, T. Andres, F. Campolongo, J. Cariboni, D. Gatelli,M. Saisana, and S. Tarantola. Global sensitivity analysis: the primer. JohnWiley & Sons, 2008. → pages 141167[171] O. Roustant, D. Ginsbourger, and Y. Deville. Dicekriging, Diceoptim: TwoR packages for the analysis of computer experiments by kriging-basedmetamodelling and optimization. J. Stat. Softw., 51(1):1, 2012. → pages141[172] W. K. Hastings. Monte Carlo sampling methods using Markov chains andtheir applications. Biometrika, 57(1):97, 1970. → pages 144[173] M. Schonlau and W. J. Welch. Screening the input variables to a computermodel via analysis of variance and visualization. In Screening, page 308.Springer, 2006. → pages 151168


Citation Scheme:


Citations by CSL (citeproc-js)

Usage Statistics



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"
                            async >
IIIF logo Our image viewer uses the IIIF 2.0 standard. To load this item in other compatible viewers, use this url:


Related Items