A QUASICLASSICAL TRAJECTORY STUDY OF H2 + H2 ENERGY
TRANSFER: A SURVEY OF APPLICABILITY OF DETAILED
BALANCE

by

MOSAMMAT HASINA AKTER
B.Sc. in Chemistry, Jahangirnagor University, 2006
M.Sc. in Chemistry, Jahangirnagor University, 2007

THESIS SUBMITTED IN PARTIAL FULFILLMENT OF
THE REQUIREMENTS FOR THE DEGREE OF
MASTER OF SCIENCE
IN
CHEMISTRY

THE UNIVERSITY OF NORTHERN BRITISH COLUMBIA
August 2016
c Mosammat Akter, 2016

Abstract
State-to-state cross sections and a survey of detailed balance for transitions
among the vibrational and rotational (v, j) states of H2 1

P

g + below 37 mEh as the

result of collisions in the H2 + H2 system were determined by using the quasiclassical
trajectory (QCT) method. Study of this system is necessary for an improved understanding of the kinetics in the interstellar medium to model processes occurring in
the molecular clouds.
The potential energy surface of Boothroyd et al. (J. Chem. Phys. 116, 666,
2002) was used for trajectory calculations. The Discrete Variable Explicit RungeKutta (DVERK) method was used for numerical integration. State-to-state cross
sections were examined for agreement with microscopic reversibility. In the majority
of cases the agreement was poor. Consequently the conclusion is that QCT is not
a viable method for the evaluation of state-to-state rate coefficients in the H2 + H2
system.

ii

Contents
Abstract

ii

Table of Contents

iv

List of Tables

v

List of Figures

xii

Acknowledgements

xiii

1 Introduction
1.1

1

The H2 + H2 system represents a useful prototype to understand interaction of two molecules . . . . . . . . . . . . . . . . . . . . . . . .

1

1.2

H2 is an important species in interstellar medium . . . . . . . . . . .

1

1.3

Computational methods in the study of molecular collisions

. . . . .

3

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

3

1.3.1

Quantum mechanical

1.3.2

Semiclassical

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

4

1.3.3

Quasiclassical trajectory (QCT) . . . . . . . . . . . . . . . . .

4

1.3.4

Advantages and disadvantages of the quasiclassical trajectory
method . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

4

1.4

Motivation for this study . . . . . . . . . . . . . . . . . . . . . . . . .

5

1.5

Research objectives . . . . . . . . . . . . . . . . . . . . . . . . . . . .

6

1.6

Thesis outline . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

6

2 Survey of H2 + H2 literature

7

2.1

H2 + H2 experimental . . . . . . . . . . . . . . . . . . . . . . . . . .

8

2.2

H4 potential Energy Surface and Calculations . . . . . . . . . . . . .

9

iii

3 Methods of calculation
3.1

14

An overview of the quasiclassical trajectory method . . . . . . . . . .

14

3.1.1

Quasiclassical trajectory method . . . . . . . . . . . . . . . .

14

3.1.2

Numerical integration . . . . . . . . . . . . . . . . . . . . . . .

15

3.1.3

Analysis of the quasiclassical trajectory results . . . . . . . . .

16

An overview of microscopic reversibility . . . . . . . . . . . . . . . . .

20

3.2.1

Microscopic reversibility . . . . . . . . . . . . . . . . . . . . .

20

3.2.2

Microscopic reversibility and rate coefficients . . . . . . . . . .

20

3.2.3

Microscopic reversibility in H + H2 system . . . . . . . . . . .

21

3.2.4

Microscopic reversibility in H2 + H2 system . . . . . . . . . .

22

3.3

Enumeration of transitions considered . . . . . . . . . . . . . . . . . .

22

3.4

Assessment of detailed balance . . . . . . . . . . . . . . . . . . . . . .

23

3.4.1

Statistical analysis . . . . . . . . . . . . . . . . . . . . . . . .

23

3.4.2

Interpolation of cross sections and errors . . . . . . . . . . . .

23

3.4.3

Histograms of internal energy distribution . . . . . . . . . . .

25

3.4.4

χ2 –test for the histograms of energy distribution . . . . . . . .

27

4 Insights from examination of microscopic reversibility calculation

28

3.2

4.1

4.2

Comparison of direct and indirect cross sections by using Z-test . . .

28

4.1.1

Z-test results . . . . . . . . . . . . . . . . . . . . . . . . . . .

28

Survey of histograms of energy distribution results . . . . . . . . . . .

29

4.2.1

Analysis of internal energy distribution . . . . . . . . . . . . .

29

4.2.2

Application of χ2 –test to histograms of energy distribution . .

49

5 Conclusions and Future Directions

53

5.1

Microscopic reversibility results in H2 + H2 system . . . . . . . . . .

53

5.2

Future directions . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

54

5.3

Concluding statement . . . . . . . . . . . . . . . . . . . . . . . . . . .

55

Bibliography

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

56

Appendix . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

59

iv

List of Tables
1.1

Collisional outcomes of H2 + H2 . . . . . . . . . . . . . . . . . . . . .

3

2.1

States below 37 mEh (1eV) . . . . . . . . . . . . . . . . . . . . . . . .

13

4.1

Detailed balance results for transitions of direct and indirect cross sections with respect to Z-test . . . . . . . . . . . . . . . . . . . . . . . .

28

4.2

Results for the cases of the Strongly Agree Class with respect to χ2 –test 50

4.3

Results for the cases of the Somewhat Agree Class with respect to χ2 –test 51

4.4

Results for cases of the Disagree Class with respect to χ2 –test . . . .

v

52

List of Figures
1.1

(v, j) levels of the ground state hydrogen molecule, H2 (1 Σ+
g ). The
dashed line is the minimum energy required for dissociation. Energy
levels above this line are quasibound. The dotted line is the internal
energy boundary for the states that are studied in this research work.

3.1

2

Schematic of binning method. The axes delineate (v 0 , j 0 ) space. The ×
indicates the trajectory result (v 00 , j 00 ). The dotted lines represent the
boundaries of the bins common to all methods.

3.2

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

18

Selected cases for three classes of agreement of direct and indirect cross
sections. The vertical axis is the cross section in a0 2 and the horizontal
axis is the translational energy in Eh . The square indicates direct cross
section and the circle is for indirect cross section. First three cases from
Strongly Agree Class for the transitions H2 (0,1) + H2 (1,4) → H2 (1,1) +
H2 (0,6), H2 (0,3) + H2 (0,6) → H2 (0,7) + H2 (0,2), and H2 (0,4) + H2 (0,9)
→ H2 (0,8) + H2 (0,5) (Top row a1 − a3 ). Second three from Somewhat
Agree Class for the transitions H2 (0,2) + H2 (0,9) → H2 (0,6) + H2 (0,5),
H2 (0,2) + H2 (1,5) → H2 (0,8) + H2 (0,3), and H2 (0,4) + H2 (0,1) →
H2 (1,2) + H2 (1,3) (Middle row b1 − b3 ). Last three from Disagree Class
for the transitions H2 (0,1) + H2 (0,2) → H2 (0,7) + H2 (0,4), H2 (0,3)
+ H2 (1,5) → H2 (0,7) + H2 (0,9), and H2 (0,3) + H2 (0,7) → H2 (0,3) +
H2 (1,3) (Bottom row c1 − c3 ). . . . . . . . . . . . . . . . . . . . . . .

vi

26

4.1

Distribution of final internal energy for the transition: H2 (0,1) + H2 (1,4)
→ H2 (1,1) + H2 (0,6) (From Strongly Agree Class). The vertical axis
is the probability and the horizontal axis is the final internal energy in
Eh . Figures a1 and a2 are for impact parameter range 0.0 - 1.0 and 1.0
- 2.0 a0 at translational energy 0.200 Eh respectively. Similarly figures
b1 and b2 are for translational energy 0.300 Eh and figures c1 and c2 are
for translational energy 0.396 Eh . . . . . . . . . . . . . . . . . . . . .

4.2

31

Distribution of final internal energy for the transition: H2 (0,6) + H2 (1,1)
→ H2 (1,4) + H2 (0,1) (From Strongly Agree Class). The vertical axis
is the probability and the horizontal axis is the final internal energy in
Eh . Figures a1 and a2 are for impact parameter range 0.0 - 1.0 and 1.0
- 2.0 a0 at translational energy 0.194 Eh respectively. Similarly figures
b1 and b2 are for translational energy 0.294 Eh and figures c1 and c2 are
for translational energy 0.390 Eh . . . . . . . . . . . . . . . . . . . . .

4.3

32

Distribution of final internal energy for the transition: H2 (0,3) + H2 (0,6)
→ H2 (0,7) + H2 (0,2) (From Strongly Agree Class). The vertical axis
is the probability and the horizontal axis is the final internal energy
in Eh . Figures a1 and a2 are for impact parameter range 0.0 - 1.0 and
1.0 - 2.0 a0 at translational energy 0.100 Eh respectively. Figures b1
to b3 are for impact parameter range 0.0 - 1.0, 1.0 - 2.0, and 2.0 - 3.0
a0 at translational energy 0.200 Eh respectively. Similarly figures c1
to c3 are for translational energy 0.300 Eh and figures d1 to d3 are for
translational energy 0.398 Eh . . . . . . . . . . . . . . . . . . . . . . .

vii

33

4.4

Distribution of final internal energy for the transition: H2 (0,2) + H2 (0,7)
→ H2 (0,6) + H2 (0,3) (From Strongly Agree Class). The vertical axis
is the probability and the horizontal axis is the final internal energy
in Eh . Figures a1 and a2 are for impact parameter range 0.0 - 1.0 and
1.0 - 2.0 a0 at translational energy 0.098 Eh respectively. Figures b1
to b3 are for impact parameter range 0.0 - 1.0, 1.0 - 2.0, and 2.0 - 3.0
a0 at translational energy 0.198 Eh respectively. Similarly figures c1
to c3 are for translational energy 0.298 Eh and figures d1 to d3 are for
translational energy 0.396 Eh . . . . . . . . . . . . . . . . . . . . . . .

4.5

34

Distribution of final internal energy for the transition: H2 (0,4) + H2 (0,9)
→ H2 (0,8) + H2 (0,5) (From Strongly Agree Class). The vertical axis
is the probability and the horizontal axis is the final internal energy
in Eh . Figures a1 and a2 are for impact parameter range 0.0 - 1.0 and
1.0 - 2.0 a0 at translational energy 0.100 Eh respectively. Figures b1
to b3 are for impact parameter range 0.0 - 1.0, 1.0 - 2.0, and 2.0 - 3.0
a0 at translational energy 0.200 Eh respectively. Similarly figures c1
to c3 are for translational energy 0.300 Eh and figures d1 to d3 are for
translational energy 0.400 Eh . . . . . . . . . . . . . . . . . . . . . . .

4.6

35

Distribution of final internal energy for the transition: H2 (0,5) + H2 (0,8)
→ H2 (0,9) + H2 (0,4) (From Strongly Agree Class). The vertical axis
is the probability and the horizontal axis is the final internal energy
in Eh . Figures a1 to a3 are for impact parameter range 0.0 - 1.0, 1.0
- 2.0, and 2.0 - 3.0 a0 at translational energy 0.102 Eh respectively.
Similarly figures b1 to b3 are for translational energy 0.202 Eh , figures
c1 to c3 are for translational energy 0.302 Eh , and figures d1 to d3 are
for translational energy 0.402 Eh . . . . . . . . . . . . . . . . . . . . .

viii

36

4.7

Distribution of final internal energy for the transition: H2 (0,2) + H2 (0,9)
→ H2 (0,6) + H2 (0,5) (From Somewhat Agree Class). The vertical axis
is the probability and the horizontal axis is the final internal energy
in Eh . Figures a1 and a2 are for impact parameter range 0.0 - 1.0 and
1.0 - 2.0 a0 at translational energy 0.100 Eh respectively. Figures b1
to b3 are for impact parameter range 0.0 - 1.0, 1.0 - 2.0, and 2.0 - 3.0
a0 at translational energy 0.200 Eh respectively. Similarly figures c1
to c3 are for translational energy 0.300 Eh and figures d1 to d3 are for
translational energy 0.396 Eh . . . . . . . . . . . . . . . . . . . . . . .

4.8

37

Distribution of final internal energy for the transition: H2 (0,5) + H2 (0,6)
→ H2 (0,9) + H2 (0,2) (From Somewhat Agree Class). The vertical axis
is the probability and the horizontal axis is the final internal energy
in Eh . Figures a1 to a3 are for impact parameter range 0.0 - 1.0, 1.0
- 2.0, and 2.0 - 3.0 a0 at translational energy 0.105 Eh respectively.
Similarly figures b1 to b3 are for translational energy 0.206 Eh , figures
c1 to c3 are for translational energy 0.305 Eh , and figures d1 to d3 are
for translational energy 0.401 Eh . . . . . . . . . . . . . . . . . . . . .

4.9

38

Distribution of final internal energy for the transition: H2 (0,2) + H2 (1,5)
→ H2 (0,8) + H2 (0,3) (From Somewhat Agree Class). The vertical axis
is the probability and the horizontal axis is the final internal energy in
Eh . Figures a1 to a3 are for impact parameter range 0.0 - 1.0, 1.0 - 2.0,
and 2.0 - 3.0 a0 at translational energy 0.200 Eh respectively. Similarly
figures b1 to b3 are for translational energy 0.300 Eh and figures c1 to
c3 are for translational energy 0.400 Eh . . . . . . . . . . . . . . . . . .

ix

39

4.10 Distribution of final internal energy for the transition: H2 (0,3) + H2 (0,8)
→ H2 (1,5) + H2 (0,2) (From Somewhat Agree Class). The vertical axis
is the probability and the horizontal axis is the final internal energy
in Eh . Figures a1 and a2 are for impact parameter range 0.0 - 1.0 and
1.0 - 2.0 a0 at translational energy 0.206 Eh respectively. Similarly figures b1 and b2 are for translational energy 0.306 Eh . Figures c1 to c2
are for impact parameter range 0.0 - 1.0, 1.0 - 2.0, and 2.0 - 3.0 a0 at
translational energy 0.407 Eh respectively. . . . . . . . . . . . . . . .

40

4.11 Distribution of final internal energy for the transition: H2 (0,1) + H2 (0,4)
→ H2 (1, 3) + H2 (1, 2) (From Somewhat Agree Class). The vertical
axis is the probability and the horizontal axis is the final internal energy in Eh . Figures a1 and a2 are for impact parameter range 0.0 - 1.0
and 1.0 - 2.0 a0 at translational energy 0.200 Eh respectively. Similarly
figures b1 and b2 are for translational energy 0.300 Eh and figures c1
and c2 are for translational energy 0.398 Eh . . . . . . . . . . . . . . .

41

4.12 Distribution of final internal energy for the transition: H2 (1,2) + H2 (1,3)
→ H2 (0,4) + H2 (0,1) (From Somewhat Agree Class). The vertical axis
is the probability and the horizontal axis is the final internal energy in
Eh . Figures a1 and a2 are for impact parameter range 0.0 - 1.0 and 1.0
- 2.0 a0 at translational energy 0.162 Eh respectively. Similarly figures
b1 and b2 are for translational energy 0.264 Eh and figures c1 and c2 are
for translational energy 0.362 Eh . . . . . . . . . . . . . . . . . . . . .

x

42

4.13 Distribution of final internal energy for the transition: H2 (0,1) + H2 (0,2)
→ H2 (0,7) + H2 (0,4) (From Disagree Class). The vertical axis is the
probability and the horizontal axis is the final internal energy in Eh .
Figures a1 and a2 are for impact parameter range 0.0 - 1.0 and 1.0 - 2.0
a0 at translational energy 0.100 Eh respectively. Figures b1 to b3 are for
impact parameter range 0.0 - 1.0, 1.0 - 2.0, and 2.0 - 3.0 a0 at translational energy 0.200 Eh respectively. Similarly figures c1 to c3 are for
translational energy 0.300 Eh and figures d1 to d3 are for translational
energy 0.402 Eh . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

43

4.14 Distribution of final internal energy for the transition: H2 (0,4) + H2 (0,7)
→ H2 (0,2) + H2 (0,1) (From Disagree Class). The vertical axis is the
probability and the horizontal axis is the final internal energy in Eh .
Figures a1 and a2 are for impact parameter range 0.0 - 1.0 and 1.0 - 2.0
a0 at translational energy 0.082 Eh respectively. Figures b1 to b3 are for
impact parameter range 0.0 - 1.0, 1.0 - 2.0, and 2.0 - 3.0 a0 at translational energy 0.182 Eh respectively. Similarly figures c1 to c3 are for
translational energy 0.282 Eh and figures d1 to d3 are for translational
energy 0.384 Eh . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

44

4.15 Distribution of final internal energy for the transition: H2 (0,3) + H2 (1,5)
→ H2 (0,7) + H2 (0,9) (From Disagree Class). The vertical axis is the
probability and the horizontal axis is the final internal energy in Eh .
Figures a1 and a2 are for impact parameter range 0.0 - 1.0 and 1.0 - 2.0
a0 at translational energy 0.099 Eh respectively. Figures b1 to b3 are for
impact parameter range 0.0 - 1.0, 1.0 - 2.0, and 2.0 - 3.0 a0 at translational energy 0.200 Eh respectively. Similarly figures c1 to c3 are for
translational energy 0.300 Eh and figures d1 to d3 are for translational
energy 0.399 Eh . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

xi

45

4.16 Distribution of final internal energy for the transition: H2 (0,7) + H2 (0,9)
→ H2 (0,3) + H2 (1,5) (From Disagree Class). The vertical axis is the
probability and the horizontal axis is the final internal energy in Eh .
Figures a1 to a3 are for impact parameter range 0.0 - 1.0, 1.0 - 2.0,
and 2.0 - 3.0 a0 at translational energy 0.192 Eh respectively. Similarly
figures b1 to b3 are for translational energy 0.292 Eh and figures c1 to
c3 are for translational energy 0.391 Eh . . . . . . . . . . . . . . . . . .

46

4.17 Distribution of final internal energy for the transition: H2 (0,3) + H2 (0,7)
→ H2 (0,3) + H2 (1,3) (From Disagree Class). The vertical axis is the
probability and the horizontal axis is the final internal energy in Eh .
Figures a1 to a3 are for impact parameter range 0.0 - 1.0, 1.0 - 2.0,
and 2.0 - 3.0 a0 at translational energy 0.200 Eh respectively. Similarly
figures b1 to b3 are for translational energy 0.300 Eh and figures c1 to
c3 are for translational energy 0.402 Eh . . . . . . . . . . . . . . . . . .

47

4.18 Distribution of final internal energy for the transition: H2 (0,3) + H2 (1,3)
→ H2 (0,3) + H2 (0,7) (From Disagree Class). The vertical axis is the
probability and the horizontal axis is the final internal energy in Eh .
Figure a1 is for impact parameter range 0.0 - 1.0 and 1.0 - 2.0 a0 at
translational energy 0.093 Eh respectively. Figures b1 to b3 are for impact parameter range 0.0 - 1.0, 1.0 - 2.0, and 2.0 - 3.0 a0 at translational
energy 0.192 Eh respectively. Similarly figures c1 to c3 are for translational energy 0.292 Eh and figures d1 to d3 are for translational energy
0.395 Eh . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

xii

48

Acknowledgements
First of all, I would like to give thanks to my supervisor Dr. Margot E. Mandy,
who has been more than just a supervisor for me. Her support, encouragement,
patience, guidance in both my academic and practical life was what made this thesis
possible. I have learned an enormous amount from working with her that I hope to
follow in my future. I would like to thank my supervisory committee members Dr.
Peter Jackson and Dr. Erik Jensen for their guidance and encouragement.
I would like to thank to the University of Northern British Columbia for providing support throughout my study. I would like to thank the Chemistry Department,
Jahangirnagor University, Dhaka, Bangladesh and Bangladesh Agricultural University, Mymensingh, Bangladesh for encouragement and support.
Thanks to all my friends and Bangladeshi community in Prince George. All of
you made living here amazing.
I am deeply grateful to my mother, Rabeya Khatun, who has always been a
source of strength and motivation in my life and specially in this work. My father,
Abdul Helim, who always inspired and encouraged me throughout my study. And my
brother, Md. Kamal Hossain, who made my life here smooth by giving me valuable
and important suggestions.
Finally, I would like to thank my lovely husband, Md. Ferdous Alam. It was
difficult for me to complete this work without all his support. And my daughter,
Fiona, who was my constant companion and source of my main strength throughout
this study.

xiii

Chapter 1
Introduction
1.1

The H2 + H2 system represents a useful prototype to understand interaction of two molecules
Collisions between two hydrogen molecules – H2 being lightest molecule in na-

ture – are among the simplest elementary chemical processes. Quantum effects are
potentially important in H2 + H2 system since the vibrational – rotational energy levels within the molecule are more widely spaced than in the case of heavier molecules
(Figure 1.1). By using this simple system of four nuclei and four electrons (the
H2 + H2 system) a high level of precision can be attained in the theoretical work of
computational chemistry which has direct applications in astrophysics. In spacecraft
modelling 1 and in combustion, 2,3 the H2 + H2 system has also been of great interest.
The rotational and vibrational transitions in H2 induced by collisions are of practical
importance in models of astrophysical environments.

1.2

H2 is an important species in interstellar medium
In the universe, hydrogen is the most abundant element. In the interstellar

medium, hydrogen is present in the dominant chemical species, mostly as H and H2 .
Shock fronts and photodissociation regions in giant molecular clouds have stimulated
considerable interest in the role of molecular energy transfer of H2 . 4,5 If collisional
cross sections and rate coefficients of H2 + H2 system can be calculated accurately
then these calculations can be used to build a model that can explain the cooling
mechanism in the interstellar medium where star formation occurs.

1

A number of collisional outcomes result from H2 + H2 collisions (Table 1.1) and
these outcomes can be distinguished in classical calculations and permit the detailed
exploration of dynamic effects. 6 The H2 + H2 system allows for theoretical treatments
that are very transparent and can be compared directly with other theoretical and
experimental results.

Figure 1.1: (v, j) levels of the ground state hydrogen molecule, H2 (1 Σ+
g ). The dashed
line is the minimum energy required for dissociation. Energy levels above this line
are quasibound. The dotted line is the internal energy boundary for the states that
are studied in this research work.

2

Table 1.1: Collisional outcomes of H2 + H2
Nonreactivea :
0
0
000
000
0000 0000
H2 (vab , jab ) + H2 (vcd
, jcd
) → H2 (vab
, jab
) + H2 (vcd
, jcd )
Reactive:
0
0
000
000
0000 0000
H2 (vab , jab ) + H2 (vcd
, jcd
) → H2 (vac
, jac
) + H2 (vbd
, jbd )
0000 0000
000
000
→ H2 (vad , jad ) + H2 (vbc , jbc )
Nonreactive Dissociative:
0
0
0000 0000
H2 (vab , jab ) + H2 (vcd
, jcd
) → H + H + H2 (vcd
, jcd )
000
000
→ H2 (vab , jab ) + H + H
Reactive Dissociative:
000
000
0
0
)+H+H
, jac
) → H2 (vac
, jcd
H2 (vab , jab ) + H2 (vcd
000
000
)+H+H
→ H2 (vad , jad
000 000
→ H2 (vbc , jbc ) + H + H
000 000
→ H2 (vbd
, jbd ) + H + H
Doubly Dissociative:
0
0
) → H+H+H+H
, jcd
H2 (vab , jab ) + H2 (vcd

a

1.3

a, b, c, and d identifies the four atoms

Computational methods in the study of molecular collisions
There are three classes of methods in the computational study of molecular

collisions. These are quantum mechanical, semiclassical, and quasiclassical trajectory
methods. All the methods require a potential energy surface (PES) for computational
study in molecule-molecule interaction.

1.3.1

Quantum mechanical

Quantum mechanical scattering calculations on four-atom reactions have undergone significant progress from the computational perspective. All degrees of freedom
are treated quantum mechanically. Like atoms are indistinguishable. Tremendous
challenge has remained for the full-dimensional quantum mechanical treatment of

3

H2 + H2 systems due to the need to treat all bonds as reactive ones. 7–13 The size
of quantum calculations rises with the number of channels required for convergence
which in turn rises with increases of total energy. Therefore, quantum calculations
can be possible in principle, but not practically, at higher energies for the H2 + H2
system and its isotopic analogs.

1.3.2

Semiclassical

Semiclassical methods describe some degrees of freedom of a system quantum
mechanically whereas the other degrees of freedom are treated classically. 14,15 These
methods are particularly useful for the description of the large-angle elastic scattering of chemical reactive systems. Examples of semiclassical methods include WKB
approximation and wave packet methods. The theory of elastic scattering of atoms
has been one of the most successful applications of semiclassical techniques. Semiclassical methods are expensive for the calculation of large number of state-to-state
cross sections in H2 + H2 system. 16

1.3.3

Quasiclassical trajectory (QCT)

Using classical mechanics, the quasiclassical trajectory method describes the
collisions of atoms and molecules. All degrees of freedom are treated classically. Like
atoms are distinguishable. This method does not increase in difficulty with increasing
energy. At the present this is the most practical method of obtaining large number
of state-to-state cross sections and rate coefficients for a system like H2 + H2 .

1.3.4

Advantages and disadvantages of the quasiclassical trajectory method

The quasiclassical trajectory method uses the classical equations of motion to
describe the collision between a target molecule and a collisional partner. Since the
trajectory is classical, there is a need to assign the trajectory result to a final quantum
4

state. Results from the QCT method can be rotationally and vibrationally hot 17 if the
trajectory results are assigned to final rotational and vibrational states by using the
bin histogram method. The cross sections of endoergic (upward) transitions calculated
by using QCT method can be shifted to higher rotational or vibrational level compared
to quantum and experimental results.
The QCT method limits the validity of rate coefficients at temperature between
600 and 10000 K. The reason for the lower limit is because QCT does not include
tunneling. There is also a consideration of zero point energy leak in the QCT method
when the potential energy barrier to exchange is traversed classically with less than
the zero point energy that would be required to do this quantum mechanically. The
H4 surface extrapolates to the H3 surface when one H atom is removed to infinite
distance. There is a conical intersection between the ground state and an electronically
excited state of H3 at energy of 100 mEh (63 kcal/mol). 18 Above 10000 K, there is a
possibility of electronic excitation and then the system is not described adequately by
the potential energy surface for the ground electronic state. 19 For H4 PES Boothroyd
et al. 20 identified the conical intersection with the first electronically excited state but
did not include it in their potential.

1.4

Motivation for this study
H2 molecules are the dominant molecular species in interstellar medium. Quan-

tum calculations are expensive for the determination of large numbers of state-to-state
cross sections and rate coefficients. The QCT method is comparatively inexpensive.
Rate coefficients must obey microscopic reversibility. The rate coefficients calculated by using QCT obeyed microscopic reversibility in H + H2 system except at
low temperatures. 19 The H2 + H2 system is a system of light atoms similar to H + H2
including one more atom and more degrees of freedom with a different potential.
Whether detailed balance is obeyed for calculation of state-to-state cross sections and
rate coefficients for H2 + H2 is important to assess the feasibility of using QCT to
5

calculate state-to-state cross sections.

1.5

Research objectives
The goal of this study is to do a survey of detailed balance of state-to-state

cross sections for the transitions resulting from the H2 + H2 system by using the
QCT method. The specific objectives of this study are:
1. Use the QCT method to determine state-to-state cross sections for the collisions
in H2 + H2 system for those states of H2 with internal energy below 37 mEh (1 eV).
2. Survey how well those cross sections follow detailed balance.
3. Assess the feasibility of the QCT method for the calculation of state-to-state rate
coefficients in H2 + H2 system.

1.6

Thesis outline
The objectives of this research are addressed in the following chapters of this

thesis. Chapter 2 is a survey of H2 + H2 literature including a survey of experiments
on the H2 + H2 system and of H4 potential energy surfaces and calculations thereon.
Chapter 3 describes the QCT method including issues for the binning method.
There is also a discussion about numerical integration and calculation of cross section
and rate coefficients in Chapter 3. Microscopic reversibility and application of microscopic reversibility in the calculation of cross sections and rate coefficients are also
discussed in this chapter. The assessment of detailed balance with statistical analysis
are included in Chapter 3.
Chapter 4 reports on the examination of the microscopic reversibility calculations for H2 + H2 system and results with the survey of internal energy distribution.
The results of internal energy distributions and their significance of differences are
also included in Chapter 4. Chapter 5 concludes with the major findings of this study
and provides suggestions for the future research directions.
6

Chapter 2
Survey of H2 + H2 literature
Collisions that occur between two H2 molecules in interstellar molecular clouds
are of great relevance in astrophysics. In the interstellar medium, H2 may act as a
coolant in giant molecular clouds of low densities where star formation occurs. Within
these clouds, strong shock waves can cause rotational and vibrational excitation of
the H2 molecules which can lead to dissociation of H2 into free H atoms 21 or to the
emission of photons. McCaffery and Marsh used a computational model of energy flow
to study translation-to-internal energy conversion in a gas ensemble consisting of H2
in a bath of H atoms. 22 Their study suggested that the translation-to-internal energy
conversion within H2 could enhance the cooling mechanism in interstellar clouds. The
presence of H atoms has the potential to change the kinetics of the system because H
atoms are more efficient colliders with respect to energy transfer and dissociation than
are H2 molecules. Dissociative cooling arises from the collisions leading to dissociation
of the molecules. When a molecule dissociates, some of the kinetic energy and internal
energy are used in dissociation and this energy can never again be available because
of the low probability of recombination. The cloud cools after losing some energy in
dissociation.
The cooling processes can also involve the conversion of kinetic energy to radiant
energy as photons which can escape from the system. This typically occurs through
collisional excitation, followed by radiative decay through quadrupole emission. Cooling is necessary for star formation. Information about these phenomena is necessary
to understand the cooling mechanisms in the interstellar medium.

7

2.1

H2 + H2 experimental
In 1975, the vibrational relaxation rate of ortho- and para-H2 between 40 and

500 K was measured by Audibert et al . 23 Nearly 12 years after this experiment, the
rotational relaxation in vibrationally excited H2 and D2 in collisions with H2 , D2
and He at 300 K was studied by Meier et al. 24 in 1986. From their studies it was
found that, in general, the pure rotational energy transfer rate coefficients are larger
than the corresponding vibrational energy transition rate coefficients. For the v = 2
vibrational level, the direct measurement of energy transfer of vibration-vibration and
vibration-translation rate coefficients near 300 K were reported by Kreutz et al. 25 in
1988.
Farrow and Chandler 26 presented experimental results in the same year for the
rovibrational energy transfer rate coefficients between H2 molecules at 295 K colliding
in their ground and first excited vibrational levels. Theoretical and experimental rate
coefficients for rotational excitations in p-H2 – p-H2 collisions between 2 and 110 K
were reported by Maté et al. 27 in 2005. Quantum mechanical calculations based on the
coupled channel method were used for their theoretical study. From their experimental
studies, it was found that the pure rotational energy transfer rate coefficients are larger
than the corresponding vibrational energy transfer rate coefficients.
Experimental rate coefficients for vibrational-vibrational transitions for v=0-5
at 300 K were reported by Ahn et al. 28 in 2007. They studied both non-resonant
and resonant VV process. For non-resonant VV process they studied H2 (v=1) +
H2 (v=1) → H2 (v=2) + H2 (v=0) and for resonant VV process, H2 (v=1) + H2 (v=0)
→ H2 (v=0) + H2 (v=1). For both processes, the measured rate coefficients were found
to be comparable to previous experimental results of Kreutz et al. 25 and Farrow and
Chandler. 26
Some significant differences exist between theory and experiment for both the
non-resonant and resonant vibration-vibration transitions discussed above. The semi-

8

classical results for the resonant vibration-vibration process of Ahn et al. 28 are found
to be about a factor of 30 smaller than the experimental value. The non-resonant
vibration-vibration process (1,1 → 2,0) is more efficient than the resonant process predicted by semiclassical calculation while experiment showed the latter rate coefficient
is 2.5 times larger. 28
In 2009, Kelley 29 studied vibration-vibration and vibration-translation energy
exchange in H2 + H2 collisions. He reported results of VV and VT using a semiclassical approach in which relative translational motion is treated classically and both
rotation and vibration were treated quantum mechanically. The dependence of the
H4 interaction potential on the H-H separations was modeled based on a He-H2 potential and transition probabilities were calculated by using first order perturbation
theory. For the resonant vibration-vibration process described above, this approach
produced improved agreement with experiments but it is less satisfactory for the other
non-resonant vibration-vibration processes for the diatom-diatom system.

2.2

H4 potential Energy Surface and Calculations
A potential energy surface is generally developed within the adiabatic and Born-

Oppenheimer approximation to model chemical reactions and interactions in simple
chemical and physical systems. In the test case of molecule-molecule interactions, the
potential energy surface of H4 is of great importance for quantum chemistry since it is
needed for quantum calculations. Being the simplest test case, the H4 energy surface is
necessary for the basic understanding of intermolecular energy transfer and chemical
reactions between molecules. It is also important for the comparison of theoretical
predictions against experimental results.
For a detailed study of the calculation of recombination rate coefficients for
hydrogen containing compounds in the gas phase, Schwenke 2 constructed a new potential hypersurface for H4 in 1988. This new potential was a faithful representation
for ab initio electronic structure calculations at that time. They used high accuracy
9

level ab initio energies at 92 geometries with C2v symmetry. This potential is unchanged under the exchange of H atoms and reproduces the H3 potential when one
H atom is removed. They used resonance complex theory and an energy transfer
mechanism over the temperature range 100-5000 K to estimate the rate coefficients
of three-body recombination by using quasiclassical trajectory calculations. Over a
broad temperature range, the predicted rate coefficients were approximately a factor
of two smaller than experimental results.
In 1991 Boothroyd et al. 30 reported an extensive study of the H4 PES presenting ab initio energies at 6101 geometries. They used the multiple reference configuration interaction (MRD-CI) program with a large basis set. In 1992, Wind and
Roeggen calculated the potential energy surface of H4 within the rigid rotor approximation by using an extended geminal model. 31 The total electronic energies for 16
different internuclear distances from 3 to 12 au, and 7 relative orientations for each
distance were calculated by using a [8s, 4p, 2d] subcontracted Gaussian type basis
set and the numerical model EXRHF3. The semiclassical calculations of Billing and
co-workers 32 16 33 34 were the most widely cited results of that time. They used an
angle-averaged potential based on the functional form of the potential energy surface
suggested by Schwenke. 2
Using an extension of Schwenke’s 2 approach to fitting the H4 potential this
was carried out both by Keogh 35 and by Aguado 36 et al.. The initial 6101 points
of Boothroyd et al. 30 were fitted to a six-dimensional many-body expansion form by
Aguado, Suarez, and Paniagua (ASP) 36 in 1994. They performed a four-body generalization of the global fitting procedure proposed by Aguado et al. 37 for triatomic
systems. With respect to the permutations of the H atoms this global H4 potential energy surface is totally symmetric and satisfied the strategies that are used in quantum
scattering calculations.
A slightly-improved version of the PES of Keogh 35 using an extension of
Schwenke’s approach to fitting H4 was developed by Martin et al. 38 and used by

10

Mandy et al. 6 for quasiclassical trajectory calculations of H2 + H2 collisions.
Another potential energy surface has been calculated by Diep and Johnson for
the rigid monomer model of (H2 )2 39 in 2000. They have employed the complete basis
set limit using coupled-cluster theory with single, double and perturbational triple excitations [CCSD(T)]. A four-term spherical harmonics expansion model was selected
to fit the surface. This potential energy surface generated the quadrupole moment to
within 0.58% and the experimental Van der Waals well depth to within 1%. From
the fitted potential energy surface the second virial coefficient has been computed.
In the temperature range of 100-500 K, the semiclassical treatment of quantum mechanical effect on the second virial coefficient was utilized. By combining Feynman’s
path integral formulation and Monte Carlo integration, they have constructed a new
technique for computing the quantum second virial coefficient and compared it with
published experimental measurements. Diep and Johnson calculated the second virial
coefficient from this PES and it showed good agreement with the experimental results
from 15-500 K.
In 2002, a chemically accurate potential energy surface was calculated by
Boothroyd et al. 20 which demonstrated significant inaccuracies in other analytic H4
surfaces available at that time. They had computed 42079 new ab initio H4 energies
and the original 6101 ab initio energies to improve the coverage of the H4 conformation space. A large portion of conical intersection of the ground state with the first
excited state and both the ground state energy and first few electronically excited
state energies were calculated. The 6101 ab initio H4 energies from Boothroyd et al. 30
provided information to fit an H4 PES that could approach chemical accuracy.
Boothroyd, Martin, Keogh, and Peterson (BMKP2) 20 reported a more elaborate fit of the 48180 ab initio points which is of chemical accuracy over a large
configuration space. The PES of Boothroyd et al. 30 provided for the first time full
coverage of the entire six-dimensional conformation space of H4 that can be reached
by a pair of H2 molecules colliding with sufficient energy to dissociate one or both

11

of the molecules. They demonstrated that the PES of Schwenke 2 had errors of more
than an order of magnitude in the repulsive wall of H4 interaction region.
A six-dimensional potential energy surface for the (H2 )2 dimer was constructed
by Hinde 40 based on coupled cluster electronic structure calculations in 2008. Large
atom-centered Gaussian basis sets and a small set of midbond functions were employed at the center of mass of the dimer. The bound and quasibound states of the
dimer are described by this surface. A close-coupled approach was used to compute
the energies of these bound and quasibound dimer states. They have compared the
computed energies for infrared and Raman transitions involving rovibrational levels
(v, j)= (0,0), (0,2), (1,0) and (1,2) with experimentally measured transition energies.
They have used four of the experimentally measured dimer transition energies to make
two empirical adjustments to the ab initio potential energy surface.
Among these potential energy surfaces, only those of Aguado et al., Boothroyd
et al., and Hinde 40 are global potential surfaces. Aguado et al. calculated a global
potential energy surface of 6101 ab initio energies with root mean square error (rms)
value of 2 mEh . An accurate global potential energy surface was also calculated by
Boothroyd et al. with an rms error of 1.43 mEh . Based on coupled-cluster electronic
structure calculations a six-dimensional potential energy surface for the H4 was constructed by Robert J. Hinde with the objective of achieving spectroscopic accuracy.
For quantum calculations, the ASP PES is more efficient because integration
over the angular coordinates can be done analytically. In QCT calculations, derivatives of the potential are needed to obtain the forces and to do numerical integration
of the equations of motion. Analytical derivatives are not available for the Hinde
or ASP PES. For the BMKP2 PES analytical derivatives are available. Calculating
a trajectory that involves numerical evaluations of the derivatives requires about 10
times as much computer time as does a trajectory where analytical derivatives are
available.
In this work, transitions of the states below 37 mEh (1 eV) energy of two hydro-

12

gen molecules are studied by using BMKP2 PES (Table 2.1). The state-to-state cross
sections for transitions among the vibrational and rotational (v, j) states of (0,0),
(0,1), (0,2), (0,3), (0,4), (0,5), (0,6), (0,7), (0, 8), (0, 9), (1,0), (1,1), (1,2), (1,3),
(1,4), (1,5) of each H2 molecule as the result of collisions in H2 + H2 system were
determined by using quasiclassical trajectory method.
Table 2.1: States below 37 mEh (1eV)
State

Para

Ortho

v
0
0
0
0
0
1
1
1
0
0
0
0
0
1
1
1

j
0
2
4
6
8
0
2
4
1
3
5
7
9
1
3
5

Eint (mEh )
9.892770
11.50692
15.21663
20.89275
28.34991
28.85506
30.38861
33.91218
10.43247
13.10638
17.81936
24.41199
32.67760
29.36784
31.90798
36.38355

13

Eint (eV)
0.267929187
0.311645695
0.412116981
0.565845312
0.767810062
0.781491207
0.823024961
0.918454793
0.282546005
0.354964485
0.482607774
0.661158062
0.885018332
0.79537891
0.864174376
0.985387811

Chapter 3
Methods of calculation
3.1

An overview of the quasiclassical trajectory
method

3.1.1

Quasiclassical trajectory method

Using classical mechanics, it is possible to simulate collisions of chemical species
if the interaction energies are described by a well-behaved potential energy surface.
The atoms are moving under the influence of an interaction potential and are assumed
to be point masses moving according to Newton’s equations of motion. The first
derivatives of the potential are used to determine the forces involved. By integrating
the equations of motion, a trajectory of a collision can be determined.
By using the BMKP2 20 potential energy surface, fully six-dimensional quasiclassical trajectories have been run. For each hydrogen molecule, the initial rotational
and vibrational energy were assigned precisely to the corresponding quantum states
while the initial translational energy was fixed at the desired value. For the collisions
of H2 + H2 , the maximum value of total energy is 480 mEh (13 eV). Above this energy
the electronic ground state potential energy surface is no longer an adequate description of the system. From this total energy, the translational energy is decreased in
intervals of 8 mEh (0.220 eV) until 30 mEh (0.866 eV) above that required for dissociation of both H2 molecules is reached. Translational energies were then decreased
by 2 mEh (0.0433 eV) from that total energy down to ensure adequate characterization of the threshold for dissociation and energy transfer. Translational energy was
decreased until all trajectories were observed to be elastic. For the impact parameter,
stratified sampling is used.
14

3.1.2

Numerical integration

After initialization, the coordinate system has been transformed from that used
for selection of initial conditions to one more appropriate for integration. The forces
are determined by the derivatives of the interaction potential and the motion of the
atoms was treated as completely classical. Using a variable step Runge-Kutta method
the equations of motion are integrated. Throughout the integration, conservation of
energy and of momentum is monitored.

3.1.2.1

Runge-Kutta Method

The Runge-Kutta technique of any order is one of the general methods for
constructing stable and efficient reaction integrators. For solving ordinary differential
equations (ODEs), Runge-Kutta (RK) methods are among of the most efficient classes
of methods. Some of these Runge-Kutta methods have excellent stability properties
and they are used widely to solve stiff ODEs. 41

3.1.2.2

Predictor-Corrector Method

In numerical analysis, a predictor-corrector method is an algorithm that proceeds in two steps. First, the prediction step calculates a rough approximation of the
desired quantity. Second, the corrector step refines the initial approximation using
another means. Predictor-corrector (PC) methods have been used as one of the major classes of methods on parallel computers for solving non-stiff ODEs. A desirable
feature of a multistep method is that the local truncation error can be determined
and a correction term can be included, which improves the accuracy of the answer at
each step.
Predictor-corrector methods need a well-behaved function for integration. The
potential for the H2 + H2 system was shown not to be sufficiently well-behaved for
integration by this method. 35 Therefore, the DVERK program is used in this present
research.
15

3.1.2.3

Discrete Variable Explicit Runge-Kutta (DVERK) Program

The DVERK Program is a double precision subroutine for solving systems
of first-order ordinary differential equations with initial conditions based on RungeKutta formulas of order 5 and 6. This program attempts to keep the global error
within a tolerance specified by the user. The magnitude of the relative error depends
on three factors. These are: the kind of error control that is used, the function being
integrated, and the range of integration i.e. step size.
There are also various options in this program including different kinds of error
control, restrictions on step sizes, and interrupts that permit the examination of the
state of the calculation as automatic modifications are possible during intermediate
stages. The equations of motion are integrated using a variable step Runge-Kutta
integrator. In this program, any options can be either initiated or altered from a previous selection prior to re-entry to the subroutine. It does not need user intervention
to successfully complete a trajectory.

3.1.3

Analysis of the quasiclassical trajectory results

The behaviour of a quasiclassical trajectory is completely described by classical
mechanics. Therefore some conventions need to be used in assigning the trajectory
result to the final quantum states. The atoms are paired for analysis of final momenta
and energy on the basis of the shortest interatomic distance. Angular momentum is
determined for each pair and then within the pair, motion along the line of centers is
considered and examined to determine whether the pair is bound. If the pair is bound,
then quasiclassical quantum numbers are calculated. This produces continuously
valued (v 00 , j 00 ) quantum numbers, which are then assigned to the discrete (v 000 , j 000 )
and (v 0000 , j 0000 ) states by using a binning method (Figure 3.1). In this study, four
dimensional bins are used, based on the final vibrational and rotational quantum
numbers of both H2 molecules.

16

3.1.3.1

Binning Method

The cross sections can be calculated by using the bin histogram method. 42 In
this method, the continuously-valued (v 00 , j 00 ) state values are assigned to a product
(v 000 , j 000 ) quantum state in accordance with the following criteria:
(v 000 − 1/2) ≤ v 00 < (v 000 + 1/2)
(j 000 − 1) ≤ j 00 < (j 000 + 1) (when j 000 6= 0)
0 < j 00 < 1 (when j 000 = 0)
where j 000 is odd for ortho and j 000 is even for para.
The v 000 value was taken as the nearest integer to v 00 .
In non-reactive collisions, ortho-para interconversion is forbidden and when the
initial values of j are even, then only the even j 000 values are permitted. Similarly,
when the initial j values are odd, then only the odd j 000 values are permitted in nonreactive collisions. Thus from a para-para initial state, non-reactive trajectories lead
only to para-para final states. Similarly non-reactive trajectories from an ortho-ortho
initial state lead only ortho-ortho products and from an ortho-para state to ortho-para
products.
The exchange reaction produces both para and ortho products, so the results
of exchange trajectories must be assigned appropriately to the para and ortho states.
The classical calculation does not include spin but the 3:1 ratio for nuclear spin is
considered in the calculation of cross sections. The reactive trajectories are binned
1
is used to calculate
four times. When the final state is ortho-ortho, a weighting of 16
9
the cross sections. When the final state is para-para, a weighting of 16
is used. When
3
the final state is ortho-para, a weighting of 16
is used to calculate the cross sections.
3
When the final state is para-ortho, a weighting of 16
is used.

17

Figure 3.1: Schematic of binning method. The axes delineate (v 0 , j 0 ) space. The ×
indicates the trajectory result (v 00 , j 00 ). The dotted lines represent the boundaries of
the bins common to all methods.

18

3.1.3.2

Calculation of Cross Sections

When a batch of trajectories has been converged with respect to impact parameter the integral cross sections, σ, can be determined by using following formula: 19
v


uP
u
1
1
2 Pi (1−Pi )
+ N2 + N3
u i Ai
X
Ni
u
i
i


σ ± dσ =
Pi A i ± u
t
3
2 2

(1 + Ni )(1 + Ni )

i

By summing over the annuli of area Ai , the cross sections and their errors are obtained.
The number of trajectories run in annulus i is Ni . Here ni represents the number of
trajectory events that lead to the outcome of interest. The probability is Pi = Nnii .
For non-negligible probabilities, Pi >> N1i and the formula for the error reduces to
that of a binomial distribution.

3.1.3.3

Calculation of Rate Coefficients

For a gas at equilibrium, the distribution of thermal energy is described by
the Maxwell-Boltzmann distribution. To determine rate coefficients, the MaxwellBoltzmann distribution is convolved with the excitation function, σ(E). The thermal
rate coefficient, γ(T ), is: 19
8kT
γ(T ) =
πµ


1/2 Z +∞
E◦

E
−E
E
σ(E)exp
d
kT
kT
kT


 



The value of γ(T ) is sensitive to E◦ which must be determined. For endoergic processes E◦ is the energy difference between the final and initial states and, for
exoergic processes, E◦ is zero. At lower temperatures, the most significant contributions to γ(T ) come from the cross sections near threshold. Cross sections need to be
determined at finely-spaced translational energies near threshold.

19

3.2

An overview of microscopic reversibility

3.2.1

Microscopic reversibility

Microscopic reversibility or detailed balance relates the cross section for a transition from state 1 to state 2 to that for the transition from state 2 to state 1. The
equation is:
σ(1 → 2; E1,rel )g1 E1,rel = σ(2 → 1; E2,rel )g2 E2,rel
where E2,rel = E1,rel − ∆E, ∆E is E2,int − E1,int , that is the difference of internal
energies in the states, g is the statistical weight for each state, and E1,rel and E2,rel
represent translational energies for respective states. Therefore, for each transition
connecting states 1 and 2, two independent estimates of the cross section can be
obtained.
The cross sections which are obtained from trajectories originating in state 1
are defined as direct cross sections and those obtained from the trajectories originating in state 2 are defined as indirect cross sections. In QCT calculations, microscopic
reversibility is not expected to be obeyed rigorously because of the way of the quantization is imposed. Direct trajectories start precisely at the initial (v, j) and (v 0 , j 0 )
bin and they end anywhere of the final (v 000 , j 000 ) and (v 0000 , j 0000 ) bin. Indirect trajectories are the reverse and may be thought of as starting precisely at the (v 000 , j 000 ) and
(v 0000 , j 0000 ) bin and ending somewhere in (v, j) and (v 0 , j 0 ) bin.

3.2.2

Microscopic reversibility and rate coefficients

Microscopic reversibility can be applied with a thermal distribution of translational energy. The rate coefficients obey:
∆E
γ(1 → 2)g1 = γ(2 → 1)g2 exp −
KT


20



It is then possible to obtain two estimates of rate coefficients for a particular transition.
The first is directly from state 1 to state 2. The second one is to apply microscopic
reversibility or detailed balance to the rate coefficient from state 2 to state 1 to obtain
a rate coefficient from state 1 to state 2. These two estimates may be combined to give
improved rate coefficients that obey detailed balance. Microscopic reversibility must
be obeyed if equilibrium is to be attained. If it is not obeyed, the rate coefficients will
not be able to be used in a meaningful way.

3.2.3

Microscopic reversibility in H + H2 system

It has been shown that for H + H2 system, the rate coefficients may differ when
microscopic reversibility is applied. 17 At high temperature, microscopic reversibility
was obeyed well for rate coefficients but at low temperature it showed discrepancies.
Hence, after applying microscopic reversibility the rate coefficients determined are not
precisely equal to those generated directly from cross sections with a piecewise linear
excitation function convolved with the Maxwell-Boltzmann distribution. At high
temperature, the forward and reverse rate coefficients may be combined to produce
rate coefficients that rigorously obey microscopic reversibility. This fact suggests
that the application of microscopic reversibility can be used to seek to improve the
evaluation of cross sections and rate coefficients in the H2 + H2 system.
The quasiclassical cross sections for a transition in the upward or endoergic
direction are susceptible to a systematic error in cross sections near threshold leading
to an over-estimation of their value. This is due to trajectory results being assigned to
a particular (v, j) box when the center of the box is not energetically accessible (See
figure 3.1). This is known as the barely-in-the-box phenomenon. This leads ”hotness”
in calculating cross sections by using QCT.

21

3.2.4

Microscopic reversibility in H2 + H2 system
The H2 + H2 system is more complicated system than the H + H2 system.

For H2 + H2 system, four quantum numbers are required to describe each of the
initial and final states. This means that a four-dimensional bin is used instead of a
two-dimensional bin to assign the trajectory outcome to a state. In this system, the
transition of one molecule may be downward in energy and other one is upward in
energy. Therefore, the applicability of microscopic reversibility may be a challenge
for this system with more degrees of freedom. The rate coefficients of a system must
obey microscopic reversibility. The applicability of detailed balance to cross sections
for the H2 + H2 system is assessed in this study.

3.3

Enumeration of transitions considered
There were 16 states considered for each hydrogen molecule below 37 mEh (1

eV) energy (Table 2.1). Among these 16 states, there are 136 possible combinations
of quantum numbers for two hydrogen molecules. There are 36 para-para, 36 orthoortho, and 64 ortho-para or para-ortho combinations for two hydrogen molecules to
be considered.
For inelastic non-reactive collisions, transitions are restricted to para-para, orthoortho, and ortho-para or para-ortho states. There are total of 666 [S=36(36+1)/2=666]
transitions for the 36 para-para states. There are total of 666 [S=36(36+1)/2=666]
transitions for the 36 ortho-ortho states. There are total of 2080 [S=64(64+1)/2=2080]
transitions for the 64 ortho-para or para-ortho states. In one direction, the total transitions are 3276 (3412-136=3276) with exclusion of 136 elastic transitions. For both
directions (endoergic and exoergic), with the exclusion of elastic transitions, a total
of 6552 (3276×2=6552) transitions are considered.

22

3.4

Assessment of detailed balance

3.4.1

Statistical analysis
Since a total of 6552 transitions is considered, there is a need for a robust

method to classify whether the cross sections for a transition are following detailed
balance. The translational energy grids of endoergic and exoergic cross sections do
not match each other exactly. In order for detailed balance to be applicable at the
cross section level, the total energy must match.

Etot = E1,int + E1,rel = E2,int + E2,rel

E2,rel = E1,rel + E1,int − E2,int
The translational energy grid E1,rel of one set of cross sections is taken as the
reference set and the other set of cross sections is interpolated to obtain cross sections
and errors at the corresponding E2,rel .

3.4.2

Interpolation of cross sections and errors

The following mathematical calculations were used to evaluate the indirect cross
section at the desired energy to assess the applicability of detailed balance to the cross
sections. The cross section, σ, is assumed to be related to translational energy, E, by
the following piecewise linear relationship where E1 < E2,rel < E2 :
σ= slope×E + intercept


Slope =

σ2 −σ1
E2 −E1



Intercept = σ1 − E1 × Slope
= σ 1 − E1 ×


=


=



σ2 −σ1
E2 −E1

E2 × σ1 −E1 × σ1
E2 −E1

E2 σ1 −E1 σ2
E2 −E1





−



E1 × σ2 −E1 × σ1
E2 −E1





23

Therefore, σ ∗ will be:
σ∗ =



σ2 −σ1
E2 −E1



× E∗ +



E2 σ1 −E1 σ2
E2 −E1



=F (σ1 , σ2 )
Therefore to get dσ ∗ , σ ∗ as a function of σ1 and σ2 is taken into account.
F (σ1 , σ2 ) = σ ∗


∂F
∂σ1



=




E2 −E ∗
E2 −E1

∂F
∂σ2



=



=



−E ∗
E2 −E1

E ∗ −E1
E2 −E1



+

E2
E2 −E1







=



E∗
E2 −E1



−



E1
E2 −E1





∗ 2

Therefore (dσ ) =



E2 −E ∗
E2 −E1

2



2

(dσ1 ) +

E ∗ −E1
E2 −E1

2

(dσ2 )2

After obtaining cross sections and errors for corresponding energies, the Z-test
is used. The Z-test statistic is:

Z=

σ − σ∗

!

q

(dσ)2 + (dσ ∗ )2

If σ and σ∗ are two independent measurements of cross sections for the same transition
at the same total energy and they agree then 68% of the time Z will be between +1
and -1 and 95% of the time between +2 and -2. For large batches of trajectories
Ni is large and the binomial distribution approaches the normal distribution when
Pi >> N1i (Section 3.1.3.2). A typical data set consists of between 180 and 300
points. Results for each data set are classified as follows:
-Direct and indirect cross sections for a transition that agree with detailed
balance within one standard deviation 68% of the time and within the
two standard deviations 95% of the time are classified as Strongly Agree (STA).
-Direct and indirect cross sections for a transition that agree with detailed
balance within one standard deviation 68% of the time or within the two
24

standard deviations 95% of the time are classified as Somewhat Agree (SWA).
-Direct and indirect cross sections for a transition that neither agree with
detailed balance within one standard deviation 68% of the time nor within
the two standard deviations 95% of the time are classified as Disagree (DIA).

3.4.3

Histograms of internal energy distribution
Three classes of agreement with detailed balance were categorized from the

assessment of agreement of direct and indirect cross sections. There may be a possibility for similar internal energy distributions within each class. If there are any
distinct energy distributions among the three classes of agreement then that could
explain why some of the cross sections followed detailed balance and others did not.
To investigate the distributions of internal energy among these three classes of agreement with detailed balance, the distribution of final internal energy is examined for
selected transitions (Figure 3.2).
In this test, there are 20 bins of internal energy that used to generate histograms
of energy distributions. The total energy range used is 66.63 mEh to cover the internal
energy span of the four-dimensional bin associated with the four quantum numbers
of two hydrogen molecules. Four translational energies (approximately 100, 200, 300,
and 400 mEh ) were examined.

25

Figure 3.2: Selected cases for three classes of agreement of direct and indirect cross
sections. The vertical axis is the cross section in a0 2 and the horizontal axis is the
translational energy in Eh . The square indicates direct cross section and the circle
is for indirect cross section. First three cases from Strongly Agree Class for the
transitions H2 (0,1) + H2 (1,4) → H2 (1,1) + H2 (0,6), H2 (0,3) + H2 (0,6) → H2 (0,7)
+ H2 (0,2), and H2 (0,4) + H2 (0,9) → H2 (0,8) + H2 (0,5) (Top row a1 − a3 ). Second
three from Somewhat Agree Class for the transitions H2 (0,2) + H2 (0,9) → H2 (0,6) +
H2 (0,5), H2 (0,2) + H2 (1,5) → H2 (0,8) + H2 (0,3), and H2 (0,4) + H2 (0,1) → H2 (1,2)
+ H2 (1,3) (Middle row b1 − b3 ). Last three from Disagree Class for the transitions
H2 (0,1) + H2 (0,2) → H2 (0,7) + H2 (0,4), H2 (0,3) + H2 (1,5) → H2 (0,7) + H2 (0,9), and
H2 (0,3) + H2 (0,7) → H2 (0,3) + H2 (1,3) (Bottom row c1 − c3 ).
26

3.4.4

χ2 –test for the histograms of energy distribution

To investigate the significance of differences in distributions of internal energy
distributions among the three classes of agreement of direct and indirect cross sections,
the χ2 -test is used. The a priori distribution is taken as the Gaussian distribution
centered on exact quantum energy of the final state. χ2 values were found from the
differences of observed and expected values of the histogram bins. Selected cases from
the three classes of agreement of direct and indirect cross sections were examined to
assess the significance of differences in energy distributions.

27

Chapter 4
Insights from examination of microscopic reversibility calculation
4.1

Comparison of direct and indirect cross sections by using Z-test

4.1.1

Z-test results

A total of 6552 transitions were examined of which 472 (242 in the endoergic
and 230 in the exoergic direction) were classified as in the STA Class, 590 (301 in
the endoergic and 289 in the exoergic direction) as in the SWA Class, and 5490
(2733 in the endoergic and 2757 in the exoergic direction) as in the DIA Class (Table
4.1). The number of points of the translational energy grids for the transition in one
direction are not equal to the number of points of the translational energy grid in the
other direction for the transition. The justification for this energy grid lies in how
the rate coefficients are calculated (Section 3.2.1). As a result of this algorithm the
translational energy grids do not match precisely for endoergic and exoergic direction.
This means that a transition in one direction may be classified differently than the
same transition in the other direction.
Table 4.1: Detailed balance results for transitions of direct and indirect cross sections
with
respect to Z-test
```
```
``` Exoergic strongly agree somewhat agree disagree total
```
Endoergic
``
`
strongly agree
170
50
10
230
somewhat agree
61
166
62
289
disagree
11
85
2661
2757
total
242
301
2733
6552
28

4.2

Survey of histograms of energy distribution results

4.2.1

Analysis of internal energy distribution

From each of the three classes of agreement three cases (for a total of nine cases)
were selected based on the similar relative errors of the direct and indirect cross
sections for the examination of histograms of internal energy distribution (Section
3.4.3). The final internal energy distribution for the transition H2 (0,1) + H2 (1,4) →
H2 (1,1) + H2 (0,6) was examined from the STA Class (Figure 4.1) and was found to
be normally distributed according to the χ2 -test. Additional cases for the transitions
H2 (0,6) + H2 (1,1) → H2 (1,4) + H2 (0,1) (Figure 4.2), H2 (0,3) + H2 (0,6) → H2 (0,7)
+ H2 (0,2) (Figure 4.3), H2 (0,2) + H2 (0,7) → H2 (0,6) + H2 (0,3) (Figure 4.4), H2 (0,4)
+ H2 (0,9) → H2 (0,8) + H2 (0,5) (Figure 4.5), and H2 (0,5) + H2 (0,8) → H2 (0,9) +
H2 (0,4) (Figure 4.6) from the STA Class were examined and were also found to be
normally distributed.
Similarly, the final internal energy distribution for the transition H2 (0,2) +
H2 (0,9) → H2 (0,6) + H2 (0,5) (Figure 4.7) was examined from SWA Class. The
histograms for the internal energy distribution were found to be normally distributed
for this case. The internal energy distribution of the additional cases for the transitions
H2 (0,5) + H2 (0,6) → H2 (0,9) + H2 (0,2) (Figure 4.8), H2 (0,2) + H2 (1,5) → H2 (0,8) +
H2 (0,3) (Figure 4.9), H2 (0,3) + H2 (0,8) → H2 (1,5) + H2 (0,2) (Figure 4.10), H2 (0,4) +
H2 (0,1) → H2 (1,2) + H2 (1,3) (Figure 4.11), H2 (1,2) + H2 (1,3) → H2 (0,4) + H2 (0,1)
(Figure 4.12) from the SWA class were examined. Significant deviations from the
normal distribution were found at 200 mEh for the transition H2 (0,2) + H2 (1,5) →
H2 (0,8) + H2 (0,3) (Figure 4.9). Examination of figure 4.9 (a1 ) showed that at lower
energy, the observed values for number of events was lower than was expected from
the normal distribution. At higher energy, the number of events was greater than

29

was expected from the normal distribution. There were no statistically significant
differences from the normal distribution among the other cases in internal energy
distribution in the SWA class.
From the DIA Class, the final internal energy distribution for the transition
H2 (0,1) + H2 (0,2) → H2 (0,7) + H2 (0,4) (Figure 4.13) was examined. Significant
differences in internal energy distributions at 100, 200, 300, and 402 mEh were found
for this case. Additional cases for the analysis of internal energy distribution in DIA
Class for the transitions H2 (0,4) + H2 (0,7) → H2 (0,2) + H2 (0,1) (Figure 4.14), H2 (0,3)
+ H2 (1,5) → H2 (0,7) + H2 (0,9) (Figure 4.15), H2 (0,7) + H2 (0,9) → H2 (0,3) + H2 (1,5)
(Figure 4.16), H2 (0,3) + H2 (0,7) → H2 (0,3) + H2 (1,3) (Figure 4.17), and H2 (0,3) +
H2 (1,3) → H2 (0,3) + H2 (0,7) (Figure 4.18) were examined. There were no statistically
significant differences from normal distribution among these cases.
The results for the histograms of internal energy distributions were found to
be similar for all three classes of agreement with detailed balance. There were no
statistically significant differences in internal energy distributions among the three
classes of agreement except for the following cases:
H2 (0,2) + H2 (1,5) → H2 (0,8) + H2 (0,3) (From SWA Class, Figure 4.9)
H2 (0,1) + H2 (0,2) → H2 (0,7) + H2 (0,4) (From DIA Class, Figure 4.13)

30

Figure 4.1: Distribution of final internal energy for the transition: H2 (0,1) + H2 (1,4)
→ H2 (1,1) + H2 (0,6) (From Strongly Agree Class). The vertical axis is the probability
and the horizontal axis is the final internal energy in Eh . Figures a1 and a2 are for
impact parameter range 0.0 - 1.0 and 1.0 - 2.0 a0 at translational energy 0.200 Eh
respectively. Similarly figures b1 and b2 are for translational energy 0.300 Eh and
figures c1 and c2 are for translational energy 0.396 Eh .

31

Figure 4.2: Distribution of final internal energy for the transition: H2 (0,6) + H2 (1,1)
→ H2 (1,4) + H2 (0,1) (From Strongly Agree Class). The vertical axis is the probability
and the horizontal axis is the final internal energy in Eh . Figures a1 and a2 are for
impact parameter range 0.0 - 1.0 and 1.0 - 2.0 a0 at translational energy 0.194 Eh
respectively. Similarly figures b1 and b2 are for translational energy 0.294 Eh and
figures c1 and c2 are for translational energy 0.390 Eh .

32

Figure 4.3: Distribution of final internal energy for the transition: H2 (0,3) + H2 (0,6)
→ H2 (0,7) + H2 (0,2) (From Strongly Agree Class). The vertical axis is the probability
and the horizontal axis is the final internal energy in Eh . Figures a1 and a2 are for
impact parameter range 0.0 - 1.0 and 1.0 - 2.0 a0 at translational energy 0.100 Eh
respectively. Figures b1 to b3 are for impact parameter range 0.0 - 1.0, 1.0 - 2.0, and
2.0 - 3.0 a0 at translational energy 0.200 Eh respectively. Similarly figures c1 to c3
are for translational energy 0.300 Eh and figures d1 to d3 are for translational energy
0.398 Eh .

33

Figure 4.4: Distribution of final internal energy for the transition: H2 (0,2) + H2 (0,7)
→ H2 (0,6) + H2 (0,3) (From Strongly Agree Class). The vertical axis is the probability
and the horizontal axis is the final internal energy in Eh . Figures a1 and a2 are for
impact parameter range 0.0 - 1.0 and 1.0 - 2.0 a0 at translational energy 0.098 Eh
respectively. Figures b1 to b3 are for impact parameter range 0.0 - 1.0, 1.0 - 2.0, and
2.0 - 3.0 a0 at translational energy 0.198 Eh respectively. Similarly figures c1 to c3
are for translational energy 0.298 Eh and figures d1 to d3 are for translational energy
0.396 Eh .

34

Figure 4.5: Distribution of final internal energy for the transition: H2 (0,4) + H2 (0,9)
→ H2 (0,8) + H2 (0,5) (From Strongly Agree Class). The vertical axis is the probability
and the horizontal axis is the final internal energy in Eh . Figures a1 and a2 are for
impact parameter range 0.0 - 1.0 and 1.0 - 2.0 a0 at translational energy 0.100 Eh
respectively. Figures b1 to b3 are for impact parameter range 0.0 - 1.0, 1.0 - 2.0, and
2.0 - 3.0 a0 at translational energy 0.200 Eh respectively. Similarly figures c1 to c3
are for translational energy 0.300 Eh and figures d1 to d3 are for translational energy
0.400 Eh .

35

Figure 4.6: Distribution of final internal energy for the transition: H2 (0,5) + H2 (0,8)
→ H2 (0,9) + H2 (0,4) (From Strongly Agree Class). The vertical axis is the probability
and the horizontal axis is the final internal energy in Eh . Figures a1 to a3 are for impact
parameter range 0.0 - 1.0, 1.0 - 2.0, and 2.0 - 3.0 a0 at translational energy 0.102 Eh
respectively. Similarly figures b1 to b3 are for translational energy 0.202 Eh , figures c1
to c3 are for translational energy 0.302 Eh , and figures d1 to d3 are for translational
energy 0.402 Eh .

36

Figure 4.7: Distribution of final internal energy for the transition: H2 (0,2) + H2 (0,9)
→ H2 (0,6) + H2 (0,5) (From Somewhat Agree Class). The vertical axis is the probability and the horizontal axis is the final internal energy in Eh . Figures a1 and a2 are
for impact parameter range 0.0 - 1.0 and 1.0 - 2.0 a0 at translational energy 0.100 Eh
respectively. Figures b1 to b3 are for impact parameter range 0.0 - 1.0, 1.0 - 2.0, and
2.0 - 3.0 a0 at translational energy 0.200 Eh respectively. Similarly figures c1 to c3
are for translational energy 0.300 Eh and figures d1 to d3 are for translational energy
0.396 Eh .

37

Figure 4.8: Distribution of final internal energy for the transition: H2 (0,5) + H2 (0,6)
→ H2 (0,9) + H2 (0,2) (From Somewhat Agree Class). The vertical axis is the probability and the horizontal axis is the final internal energy in Eh . Figures a1 to a3 are
for impact parameter range 0.0 - 1.0, 1.0 - 2.0, and 2.0 - 3.0 a0 at translational energy
0.105 Eh respectively. Similarly figures b1 to b3 are for translational energy 0.206 Eh ,
figures c1 to c3 are for translational energy 0.305 Eh , and figures d1 to d3 are for
translational energy 0.401 Eh .

38

Figure 4.9: Distribution of final internal energy for the transition: H2 (0,2) + H2 (1,5)
→ H2 (0,8) + H2 (0,3) (From Somewhat Agree Class). The vertical axis is the probability and the horizontal axis is the final internal energy in Eh . Figures a1 to a3 are
for impact parameter range 0.0 - 1.0, 1.0 - 2.0, and 2.0 - 3.0 a0 at translational energy
0.200 Eh respectively. Similarly figures b1 to b3 are for translational energy 0.300 Eh
and figures c1 to c3 are for translational energy 0.400 Eh .

39

Figure 4.10: Distribution of final internal energy for the transition: H2 (0,3) + H2 (0,8)
→ H2 (1,5) + H2 (0,2) (From Somewhat Agree Class). The vertical axis is the probability and the horizontal axis is the final internal energy in Eh . Figures a1 and a2 are
for impact parameter range 0.0 - 1.0 and 1.0 - 2.0 a0 at translational energy 0.206 Eh
respectively. Similarly figures b1 and b2 are for translational energy 0.306 Eh . Figures c1 to c2 are for impact parameter range 0.0 - 1.0, 1.0 - 2.0, and 2.0 - 3.0 a0 at
translational energy 0.407 Eh respectively.

40

Figure 4.11: Distribution of final internal energy for the transition: H2 (0,1) + H2 (0,4)
→ H2 (1, 3) + H2 (1, 2) (From Somewhat Agree Class). The vertical axis is the
probability and the horizontal axis is the final internal energy in Eh . Figures a1 and
a2 are for impact parameter range 0.0 - 1.0 and 1.0 - 2.0 a0 at translational energy
0.200 Eh respectively. Similarly figures b1 and b2 are for translational energy 0.300 Eh
and figures c1 and c2 are for translational energy 0.398 Eh .

41

Figure 4.12: Distribution of final internal energy for the transition: H2 (1,2) + H2 (1,3)
→ H2 (0,4) + H2 (0,1) (From Somewhat Agree Class). The vertical axis is the probability and the horizontal axis is the final internal energy in Eh . Figures a1 and a2
are for impact parameter range 0.0 - 1.0 and 1.0 - 2.0 a0 at translational energy 0.162
Eh respectively. Similarly figures b1 and b2 are for translational energy 0.264 Eh and
figures c1 and c2 are for translational energy 0.362 Eh .

42

Figure 4.13: Distribution of final internal energy for the transition: H2 (0,1) + H2 (0,2)
→ H2 (0,7) + H2 (0,4) (From Disagree Class). The vertical axis is the probability
and the horizontal axis is the final internal energy in Eh . Figures a1 and a2 are for
impact parameter range 0.0 - 1.0 and 1.0 - 2.0 a0 at translational energy 0.100 Eh
respectively. Figures b1 to b3 are for impact parameter range 0.0 - 1.0, 1.0 - 2.0, and
2.0 - 3.0 a0 at translational energy 0.200 Eh respectively. Similarly figures c1 to c3
are for translational energy 0.300 Eh and figures d1 to d3 are for translational energy
0.402 Eh .

43

Figure 4.14: Distribution of final internal energy for the transition: H2 (0,4) + H2 (0,7)
→ H2 (0,2) + H2 (0,1) (From Disagree Class). The vertical axis is the probability
and the horizontal axis is the final internal energy in Eh . Figures a1 and a2 are for
impact parameter range 0.0 - 1.0 and 1.0 - 2.0 a0 at translational energy 0.082 Eh
respectively. Figures b1 to b3 are for impact parameter range 0.0 - 1.0, 1.0 - 2.0, and
2.0 - 3.0 a0 at translational energy 0.182 Eh respectively. Similarly figures c1 to c3
are for translational energy 0.282 Eh and figures d1 to d3 are for translational energy
0.384 Eh .

44

Figure 4.15: Distribution of final internal energy for the transition: H2 (0,3) + H2 (1,5)
→ H2 (0,7) + H2 (0,9) (From Disagree Class). The vertical axis is the probability
and the horizontal axis is the final internal energy in Eh . Figures a1 and a2 are for
impact parameter range 0.0 - 1.0 and 1.0 - 2.0 a0 at translational energy 0.099 Eh
respectively. Figures b1 to b3 are for impact parameter range 0.0 - 1.0, 1.0 - 2.0, and
2.0 - 3.0 a0 at translational energy 0.200 Eh respectively. Similarly figures c1 to c3
are for translational energy 0.300 Eh and figures d1 to d3 are for translational energy
0.399 Eh .

45

Figure 4.16: Distribution of final internal energy for the transition: H2 (0,7) + H2 (0,9)
→ H2 (0,3) + H2 (1,5) (From Disagree Class). The vertical axis is the probability and
the horizontal axis is the final internal energy in Eh . Figures a1 to a3 are for impact
parameter range 0.0 - 1.0, 1.0 - 2.0, and 2.0 - 3.0 a0 at translational energy 0.192 Eh
respectively. Similarly figures b1 to b3 are for translational energy 0.292 Eh and figures
c1 to c3 are for translational energy 0.391 Eh .

46

Figure 4.17: Distribution of final internal energy for the transition: H2 (0,3) + H2 (0,7)
→ H2 (0,3) + H2 (1,3) (From Disagree Class). The vertical axis is the probability and
the horizontal axis is the final internal energy in Eh . Figures a1 to a3 are for impact
parameter range 0.0 - 1.0, 1.0 - 2.0, and 2.0 - 3.0 a0 at translational energy 0.200 Eh
respectively. Similarly figures b1 to b3 are for translational energy 0.300 Eh and figures
c1 to c3 are for translational energy 0.402 Eh .

47

Figure 4.18: Distribution of final internal energy for the transition: H2 (0,3) + H2 (1,3)
→ H2 (0,3) + H2 (0,7) (From Disagree Class). The vertical axis is the probability and
the horizontal axis is the final internal energy in Eh . Figure a1 is for impact parameter
range 0.0 - 1.0 and 1.0 - 2.0 a0 at translational energy 0.093 Eh respectively. Figures b1
to b3 are for impact parameter range 0.0 - 1.0, 1.0 - 2.0, and 2.0 - 3.0 a0 at translational
energy 0.192 Eh respectively. Similarly figures c1 to c3 are for translational energy
0.292 Eh and figures d1 to d3 are for translational energy 0.395 Eh .

48

4.2.2

Application of χ2 –test to histograms of energy distribution

From the χ2 –test, it is found that there were no significant differences in energy
distributions for almost all the cases. The cross sections for the transition of H2 (0,6)
+ H2 (1,1) → H2 (1,4) + H2 (0,1) showed significant differences at 95% confidence level
(Table 4.2).
The cross sections for the transition of H2 (0,2) + H2 (1,5) → H2 (0,8) + H2 (0,3)
exhibited significant differences at 99.95% and 95% confidence levels (Table 4.3). The
cross sections for the transition of H2 (0,1) + H2 (0,2) → H2 (0,7) + H2 (0,4) showed
significant differences at 99.95% confidence level (Table 4.4).

49

Table 4.2: Results for the cases of the Strongly Agree Class with respect to χ2 –test
Energy,
Cross section
Number of
a
b
Initial state Final state
mEh
with error, a0 2
Trajectoriesc χ2 valued
0114
1106
100
(0 ± 2.24) × 10−4
0
−3
200
(5.520 ± 0.481) × 10
148
6.50
−3
300
(9.794 ± 1.50) × 10
208
5.84
−3
396
(7.629 ± 1.76) × 10
237
4.71
1106
0114
94
(0 ± 2.24) × 10−4
0
−3
194
(3.343 ± 0.384) × 10
102
40.77
−3
294
(7.391 ± 1.13) × 10
160
31.97
390
(5.834 ± 0.164) × 10−3
177
32.31
−1
0306
0702
100
(1.061 ± 0.007) × 10
195
16.20
−1
200
(1.354 ± 0.002) × 10
1872
10.81
−1
300
(1.091 ± 0.006) × 10
208
8.50
398
(1.061 ± 0.008) × 10−1
70
10.00
0702
0306
98
(1.263 ± 0.008) × 10−1
207
18.81
−1
198
(1.656 ± 0.001) × 10
11203
29.10
298
(1.416 ± 0.007) × 10−1
224
6.65
−1
396
(1.281 ± 0.009) × 10
93
3.93
−2
0409
0805
100
(9.649 ± 0.072) × 10
153
26.40
−1
200
(1.286 ± 0.002) × 10
1923
5.32
300
(1.267 ± 0.007) × 10−1
198
5.54
−1
400
(1.023 ± 0.008) × 10
81
7.63
−2
0805
0409
102
(8.347 ± 0.068) × 10
128
15.01
202
(1.091 ± 0.003) × 10−1
754
13.23
−1
302
(1.083 ± 0.006) × 10
157
10.14
−1
402
(1.027 ± 0.008) × 10
69
12.33

Initial states for two hydrogen molecules vjv 0 j 0
Final states for two hydrogen molecules v 000 j 000 v 0000 j 0000
c
Number of Trajectories in impact parameter range 0 - 1 a0 that contribute to the cross section for
a transition
d
Italic type indicates differences are statistically significant at 95% confidence level

a

b

50

Table 4.3: Results for the cases of the Somewhat Agree Class with respect to χ2 –test
Energy,
Cross section
Number
a
b
Initial state Final state
mEh
with error,a0 2
of Trajectoriesc χ2 valuede
0209
0605
100
(1.182 ± 0.007) × 10−1
206
21.55
−1
200
(1.659 ± 0.003) × 10
2461
8.24
−1
300
(1.547 ± 0.006) × 10
2819
7.27
−1
396
(1.341 ± 0.009) × 10
2915
7.34
0605
0209
105
(8.347 ± 0.699) × 10−2
132
8.45
−1
206
(1.126 ± 0.003) × 10
750
11.76
−1
305
(1.005 ± 0.008) × 10
839
12.35
401
(9.177 ± 0.809) × 10−2
903
19.06
−3
0215
0803
100
(1.795 ± 1.22) × 10
2
4.82
−2
200
(2.088 ± 0.009) × 10
466
81.61
−2
300
(2.701 ± 0.339) × 10
535
37.80
400
(2.580 ± 0.421) × 10−2
571
23.66
−4
0803
0215
107
(2.243 ± 3.17) × 10
1
−2
206
(1.052 ± 0.110) × 10
102
5.37
306
(1.272 ± 0.207) × 10−2
144
8.78
−2
407
(1.144 ± 0.295) × 10
167
9.44
−4
0104
1312
100
(0 ± 2.24) × 10
0
−3
200
(1.165 ± 0.112) × 10
149
13.90
300
(8.168 ± 1.97) × 10−3
182
22.57
−3
398
(8.527 ± 2.54) × 10
202
25.72
−4
1312
0104
64
(0 ± 2.24) × 10
0
162
(1.366 ± 0.273) × 10−3
45
24.63
264
(4.241 ± 1.02) × 10−3
70
34.94
−3
362
(7.629 ± 2.58) × 10
83
35.95

Initial states for two hydrogen molecules vjv 0 j 0
Final states for two hydrogen molecules v 000 j 000 v 0000 j 0000
c
Number of Trajectories in impact parameter range 0 - 1 a0 that contribute to the cross section for
a transition
d
Italic type indicates differences are statistically significant at 95% confidence level
e
Bold type indicates differences are statistically significant at 99.95% confidence level

a

b

51

Table 4.4: Results for cases of the Disagree Class with respect to χ2 –test
Energy,
Cross section
Number
a
b
Initial state Final state
mEh
with error,a0 2
of Trajectoriesc
0102
0704
100
(7.584 ± 0.618) × 10−2
147
−1
200
(1.937 ± 0.001) × 10
12726
−1
300
(1.899 ± 0.008) × 10
13089
−1
402
(1.638 ± 0.102) × 10
13236
0704
0102
82
(1.054 ± 0.236) × 10−2
20
−2
182
(2.625 ± 0.120) × 10
433
−2
282
(2.261 ± 0.311) × 10
486
384
(1.638 ± 0.341) × 10−2
509
−2
0315
0709
99
(1.346 ± 0.249) × 10
33
−2
200
(6.137 ± 0.182) × 10
1089
−2
300
(6.518 ± 0.517) × 10
1235
399
(5.767 ± 0.626) × 10−2
1308
−3
0709
0315
91
(1.795 ± 1.22) × 10
2
−2
192
(1.640 ± 0.093) × 10
321
292
(2.356 ± 0.317) × 10−2
379
−2
391
(2.513 ± 0.444) × 10
404
−4
0307
0313
100
(8.975 ± 5.01) × 10
4
−2
200
(2.737 ± 0.116) × 10
606
300
(3.298 ± 0.360) × 10−2
703
−2
402
(3.635 ± 0.515) × 10
741
−3
0313
0307
93
(3.814 ± 1.78) × 10
12
192
(7.324 ± 0.198) × 10−2
1429
292
(1.116 ± 0.007) × 10−1
1627
−1
395
(1.124 ± 0.009) × 10
1722

χ2 valuede
95.11
166.93
166.15
165.65
8.40
10.47
15.63
19.57
2.48
21.69
17.28
21.64
9.33
11.85
12.08
13.51
5.18
26.65
29.02
34.42
6.01
10.64
12.14
27.77

Initial states for two hydrogen molecules vjv 0 j 0
Final states for two hydrogen molecules v 000 j 000 v 0000 j 0000
c
Number of Trajectories in impact parameter range 0 - 1 a0 that contribute to the cross section for
a transition
d
Italic type indicates differences are statistically significant at 95% confidence level
e
Bold type indicates differences are statistically significant at 99.95% confidence level

a

b

52

Chapter 5
Conclusions and Future Directions
5.1

Microscopic reversibility results in H2 + H2
system
The results of the examination of detailed balance for the collisions between two

hydrogen molecules below 37 mEh (1 eV) energy were categorized into three classes.
A total of 472 transitions of direct and indirect cross sections were found to follow
detailed balance very strongly, 590 to somewhat follow detailed balance, and 5490 not
to follow detailed balance. The internal energy distributions among the three classes
of agreement show similar normal distributions of energy except in two cases where
the differences from the normal distributions were found to be statistically significant
but were not correlated with the classification.
The agreement of detailed balance of state-to-state cross sections in H2 + H2
system appeared to be random and was not found to be predictable. In QCT calculations, microscopic reversibility is not expected to be obeyed rigorously because of the
way of the quantization is imposed. Direct trajectories start precisely at the initial
(v, j) states and they end anywhere of the final (v 0 , j 0 ) bin. Indirect trajectories are the
reverse and may be thought of as starting somewhere in the (v 0 , j 0 ) bin and precisely
ending at (v, j). The objective of this work was to survey detailed balance for the
state-to-state cross sections in H2 + H2 system by using QCT method. The survey
of detailed balance was hoped to provide an important foundation for the calculation
of cross sections and rate coefficients in H2 + H2 system. By using these calculations
a model system of H2 + H2 can be built which can further used to understand the
cooling mechanism in interstellar medium where H2 molecules are dominant chemical
species. Appraisal of the results of this work demonstrated that the QCT method
53

was not to be a feasible method for the calculation of the state-to-state cross sections
in the H2 + H2 system. Since this work began there has been further development in
quantum methods.
In 2014, the state-to-state rate coefficients for selected transitions involving
the (v=0, j≤ 10), (v=1, j≤ 8), and (v=2, j≤ 6) states were calculated by Bohr et
al. 7 quantum mechanically. They have used the coupled-states approximation and
set the maximum energy at 23 mEh (5000 cm−1 ) above the internal energy for the
initial states of two H2 molecules. They have reduced the number of coupled channels
required by using the coupled-states approximation for the initial states of two H2
molecules. They have compared their method with the close-coupling formulation for
54 p-H2 - p-H2 transitions and have found good agreement.
To build a model system of H2 + H2 in the astrophysical environment, the stateto-state rate coefficients of H2 + H2 among all (v, j) states need to be calculated. A
complete model would need to consider over 60,000 combinations of initial states
and over 35 billion possible transitions. Total energy in excess of 170 mEh would be
needed. Quantum mechanical methods could be used to achieve this if computational
power increases and improved techniques are developed.

5.2

Future directions
Additional cases with varying relative errors of direct and indirect cross sections

can be investigated to assess the significance of differences in energy distributions.
Marsh and McCaffery proposed a statistical mechanical model for the collisions in
diatom-diatom systems based on the assumption of angular momentum conservation
in rotational energy transfer. 43 They applied it to the N2 + H2 system and have
found good agreement with experiments. They have examined the final rotational
distributions associated with vibrational energy transfer for several systems. 44 Angular momentum constraints on the formation of products were studied by McCaffery et
al. 45 They have suggested strategies for the examination of exit routes from the tran54

sition state subject to angular momentum constraints. This examination suggested
that the reduction of angular momentum constraints in the exit channels from the
transition state may act as a form of catalysis. In the collisions of different diatomic
molecules, there is an asymmetry in the rotational energy levels that could promote
relaxation over excitation. The quadratic dependence of the rotational energy on
quantum number could be the source of this effect. They examined competitive
partitioning of rotational and vibrational energy in equilibration of a gas ensemble
consisting of diatomic molecules in a bath gas of different diatomic molecules. 46 They
have found that the rotational distributions of two different types of molecules were
distinct.
Angular momentum conservation effects on inelastic transitions in H2 + H2
system were examined by dos Santos et al. 47 quantum mechanically. They have
found that both energy and angular momentum performed fundamental roles in the
rovibrational transitions. For inelastic processes, the transitions that conserved the
internal rotational angular momentum and involved small changes of internal energy
of the two molecules were found to be favored.
This suggests that the angular momentum distributions for the three classes of
agreement with detailed balance could be explored further. This could be investigated
by rerunning selected batches of trajectories and retaining additional information
about angular momentum. This is beyond the scope of the current study.

5.3

Concluding statement
Detailed balance is not obeyed well for cross sections determined by QCT for

the H2 + H2 system. Therefore, it can be concluded that the QCT is not a robust
method for the evaluation of energy transfer cross sections and rate coefficients in
H2 + H2 system.

55

Bibliography
[1] C. Park, J. Thermophys. Heat Transfer, 1993, 7, 385.
[2] D. W. Schwenke, J. Chem. Phys., 1988, 89, 2076.
[3] D. W. Schwenke, J. Chem. Phys., 1990, 92, 7267.
[4] A. Sternberg and A. Dalgarno, Astrophys. J., 1989, 338, 197.
[5] P. G. Martin, D. H. Schwarz and M. E. Mandy, Astrophys. J., 1996, 461, 265.
[6] M. E. Mandy, P. G. Martin and W. J. Keogh, J. Chem. Phys., 1998, 108, 492–
497.
[7] A. Bohr, S. Paolini, R. C. Forrey, N. Balakrishnan and P. Stancil, J. Chem.
Phys., 2014, 140, 064308.
[8] S. F. dos Santos, N. Balakrishnan, R. C. Forrey and P. Stancil, J. Chem. Phys.,
2013, 138, 104302.
[9] S. K. Pogrebnya and D. C. Clary, Chem. Phys. Lett., 2002, 363, 523–528.
[10] R. A. Sultanov and D. Guster, Chem. Phys., 2006, 326, 641.
[11] M. Hankel, S. C. Smith, R. J. Allan, S. K. Gray and G. G. Balint-Kurti, J. Chem.
Phys., 2006, 125, 164303.
[12] F. Otto, F. Gatti and H. D. Meyer, J. Chem. Phys., 2008, 128, 064305.
[13] S. Y. Lin and H. Guo, J. Chem. Phys., 2002, 117, 5183.
[14] A. Panda, F. Otto, F. Gatti and H. Meyer, J. Chem. Phys., 2007, 127, 114310.
[15] R. Santiago, O. Aĺvarez-Bajo, J. M. Arias, J. Goḿez-Camacho and R. Lemus,
Mole. Phys., 2012, 110, 2003.
[16] M. Cacciatore and G. D. Billing, J. Chem. Phys., 1992, 96, 217.
56

[17] M. E. Mandy, P. G. Martin and W. J. Keogh, J. Chem. Phys., 1994, 100, 2671.
[18] A. J. C. Varandas, F. B. Brown, C. A. Mead, D. G. Truhlar and N. C. Blais, J.
Chem. Phys., 1987, 86, 6258.
[19] M. E. Mandy and P. G. Martin, J. Chem. Phys., 1999, 110, 7811.
[20] A. I. Boothroyd, P. G. Martin, W. J. Keogh and M. J. Peterson, J. Chem. Phys.,
2002, 116, 666.
[21] J. E. Dove, A. C. M. Rusk, P. H. Cribb and P. G. Martin, Astrophys. J., 1987,
318, 379.
[22] A. J. McCaffery and R. J. Marsh, J. Chem. Phys., 2013, 139, 234310.
[23] M. -M.Audibert, R.Vilaseca, J. Lukasik and J. Ducuing, Chem. Phys. Lett., 1975,
31, 232.
[24] W. Meier, G. Ahlers and H. Zacharias, J. Chem. Phys., 1986, 85, 2599.
[25] T. G. Kreutz, J. Gelfand, R. B. Miles and H. Rabitz, Chem. Phys., 1988, 124,
359.
[26] R. L. Farrow and D. W. Chandler, J. Chem. Phys., 1988, 89, 1994.
[27] B. Maté, F. Thibault, G. Tejeda, J. M. Fernández and S. Montero, J. Chem.
Phys., 2005, 122, 064313.
[28] T. Ahn, I. Adamovich and W. R. Lempert, Chem. Phys., 2007, 335, 55.
[29] J. D. Kelley, J. Phys. Chem. A, 2009, 113, 1995.
[30] A. I. Boothroyd, W. J. Keogh, P. G. Martin and M. J. Peterson, J. Chem. Phys.,
1991, 95, 4331.
[31] P. Wind and I. Roeggen, J. Chem. Phys., 1992, 167, 263–275.
[32] G. D. Billing and E. R. Fisher, Chem. Phys., 1976, 18, 225.
57

[33] R. E. Kolesnick and G. D. Billing, Chem. Phys., 1993, 170, 201.
[34] G. D. Billing and R. E. Kolesnick, Chem. Phys. Lett., 1993, 215, 571.
[35] W. J. Keogh, 1992.
[36] A. Aguado, C. Suarez and M. Paniagua, J. Chem. Phys., 1994, 101, 4004.
[37] A. Aguado and M. Paniagua, J. Chem. Phys., 1992, 96, 1265.
[38] P. G. Martin, W. J. Keogh and M. E. Mandy, Astrophys. J., 1998, 499, 793.
[39] P. Diep and J. K. Johnson, J. Chem. Phys., 2000, 112, 4465–4473.
[40] R. J. Hinde, J. Chem. Phys., 2008, 128, 154308.
[41] W. H. Press, S. A. Teukolsky, W. T. Vetterling and B. P. Flannery, Numerical
Recipes in C, Cambridge University Press, 1992, p. 710.
[42] J. E. Dove, S. Raynor and H. Teitelbaum, Chem. Phys., 1980, 50, 175.
[43] R. J. Marsh and A. J. McCaffery, Chem. Phys. Lett., 2001, 335, 134.
[44] A. J. McCaffery and R. J. Marsh, J. Phys. Chem. A, 2000, 104, 10442.
[45] A. J. McCaffery, M. A. Osborne and R. J. Marsh, J. Phys. Chem. A, 2005, 109,
5005.
[46] A. J. McCaffery and R. J. Marsh, J. Chem. Phys., 2012, 136, 024307.
[47] S. F. dos Santos, N. Balakrishnan, S. Lepp, G. Quéméner, R. C. Forrey, R. J.
Hinde and P. C. Stancil, J. Chem. Phys., 2011, 134, 214303.

58

Appendix
TERMINOLOGY AND ABBREVIATIONS USED IN THESIS
QCT-Quasiclassical Trajectory
PES-Potential Energy Surface
ASP-Aguado, Suarez, and Paniagua (J. Chem. Phys. 101, 4004, 1994)
BMKP2-Boothroyd, Martin, Keogh, and Peterson Potential Energy Surface in 2002
(J. Chem. Phys. 116, 666, 2002)
VV- Vibration-vibration
VT- Vibration-transition
DVERK- Discrete Variable Explicit Runge-Kutta
Bin Histogram Method- The continuously-valued quantum numbers is assigned to
a product state by using following criteria(v 000 − 1/2) ≤ v 00 < (v 000 + 1/2)
(j 000 − 1) ≤ j 00 < (j 000 + 1) (When j 000 6= 0)
0 < j 00 < 1 (When j 000 = 0)
j 000 is odd for ortho state and j 000 is even for para state
Endoergic Transition- A Transition for which the initial state is of lower energy
than the final state
Exoergic Transition- A Transition for which the initial state is of higher energy
than the final state
Direct Cross Sections- Cross sections for a transition from an initial state to final
state generated from trajectories originating in the initial state
Indirect Cross Sections- Cross sections for a transition from an initial state to final
state calculated from trajectories originating in the final state
DB-Detailed Balance
First Impact Parameter Stratum- Impact parameter range is 0 - 1a0 .
STA Class- Strongly Agree with Detailed Balance Class (Direct and indirect cross

59

sections for a transitions that agree within the standard deviation 68% of the time
AND within the standard deviation 95% of the time are classified as strongly agree.)
SWA Class- Somewhat Agree with Detailed Balance Class (Direct and indirect cross
sections for the transitions that agree within the standard deviation 68% of the time
OR within the standard deviation 95% of the time are classified as somewhat agree.)
DIA Class- Disagree with Detailed Balance Class (Direct and indirect cross sections
for the transitions that NEITHER agree within the standard deviation 68% of the
time NOR within the standard deviation 95% of the time are classified as disagree.)

60