Open Collections

UBC Theses and Dissertations

UBC Theses Logo

UBC Theses and Dissertations

A study of ergodicity and chaos in the Bianchi cosmologies Bhushan, Vikas 1995

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

Item Metadata

Download

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

Full Text

A STUDY OF ERGODICITY AND CHAOS IN T H E BIANCHI COSMOLOGIES By Vikas Bhushan Sc. (Mathematics and Physics) University of Toronto, 1992 A T H E S I S S U B M I T T E D I N P A R T I A L F U L F I L L M E N T ' O F T H E R E Q U I R E M E N T S F O R T H E D E G R E E O F M A S T E R O F S C I E N C E •• • in T H E F A C U L T Y O F G R A D U A T E S T U D I E S - D E P A R T M E N T O F P H Y S I C S . We accept this thesis as conforming to the required standard T H E U N I V E R S I T Y O F B R I T I S H C O L U M B I A April 1995 © Vikas Bhushan, 1995 In presenting this thesis in partial fulfilment of the requirements for an advanced degree at the University of British Columbia, I agree that the Library shall make it freely available for reference and study. I further agree that permission for extensive copying of this thesis for scholarly purposes may be granted by the head of my department or by his or her representatives. It is understood that copying or publication of this thesis for financial gain shall not be allowed without my written permission. Department The University of British Columbia Vancouver, Canada Date Aprt] I^ST DE-6 (2/88) Abstract The flow on the phase space of gravitational degrees of freedom for the spatially homoge-neous Bianchi I cosmology is studied in detail. Problems associated with analyzing the dynamics numerically are discussed and a rigorous analysis of the phase space dynamics for Bianchi I is presented. The global topology of the phase space is shown to be non-trivial and to have an important effect on the dynamics of the system. Although this results in complicated dynamics on the phase space, which is contrary to what has pre-viously been thought, the system is nevertheless shown to be non-ergodic. Implications for the study of the dynamics of other Bianchi models is discussed. 11 Table of Contents Abstract ii Table of Contents iii List of Tables v List of Figures vii Acknowledgements xvi 1 Introduction 1 2 General Relativity as a Hamiltonian System 7 2.1 3-fl Spacetime Decomposition 8 2.2 The Initial Value Formulation and A D M Formalism 11 -2.3 Hamiltonian Formulation of General Relativity . 15 2.4 The Phase Space for General Relativity 19 3 Homogeneous Cosmologies 23 3.1 Background . . . 23 3.2 The Bianchi Classification . 29 3.3 Phase Space for the Homogeneous Cosmologies 32 4 Curvature and Einstein Tensors for the Bianchi Models 34 4.1 Riemann Tensor Components 36 in 4.2 Einstein Tensor Components 45 5 Dynamical Systems 54 5.1 Concepts and Theorems 54 5.2 An Illustrative Example 68 5.3 Ergodicity and Chaos 78 6 Bianchi I : A Case Study 99 6.1 The Torus Example 100 6:2 Canonical Variables 102 6.3 Phase Space 105 6.4 Dynamics . . . 112 7 Bianchi II : A Brief Look 128 8 Discussion and Concluding Remarks 130 Bibliography 134 A Listing of Computational Program P O I N R E C . C 137 B Listing of Computational Program S E Q U E N C E . C 151 i v List of Tables 4.1 The Structure Constants Cj f c and the corresponding values of dl-k for each of the Bianchi Classes 47 4.2 The Spatial Parts of the Goo component of the Einstein tensor for each of the Bianchi Classes, in terms of dl-k and (fi, /3+, /?_). The full expression is Goo = 3 ( f i 2 - fa - - \ 3 R . . . . . 49 4.3 The G 0 3 components of the Einstein Tensor and the corresponding Mo-mentum Constraints for each of the Bianchi Classes 50 4.4 The- Spatial Parts ^Gij of the Gij components of the Einstein Tensor for each of the Bianchi Classes. 53 5.1 Average return time for 10 molecules over a varying number of runs with a 25x50 box and maximum speed for each molecule of 100 gridpoints per timestep. Different initial data was used for each run, and if the return time for any given run exceeded 100000 iterations, then that run was discarded • '. 62 5.2 The average return time rav over 1000 runs. If the maximum return time of one million iterations was exceeded for a run, or if a run had a return time of one iteration, then that run was-discarded. A value of X was. calculated according to the relation rav = XN 64 v 5.3 The minimum estimated return times Tm{n based on the data in table 5.2, calculated by assuming that all of the discarded runs have return times of one million iterations. X was calculated using r m j n according to the relation r^in = X N . ' . ' 64 5.4 The number and frequency of occurrences of each of the digits in the first 200 elements of the sequence S N , calculated using the computer program listed in the appendix running on a SUN SparclO computer . 70 5.5 The frequency of occurrences of each of the digits in the sequence S N , calculated using equation ( 5.40 ) . . . 76 6.1 Bounces for trajectories approaching the point P from region I. 90 is the bounce angle at point P. A l l angles are measured from the normal to the wall. , 118 6.2 Bounces for trajectories approaching the point P from region II. 80 is the bounce angle at point P. A l l angles are measured from the normal to the wal l . ' 119 6.3 Bounces for trajectories approaching the point P from region III. 9Q is the bounce angle at point P. A l l angles'are measured from the normal to the wall 119 6.4 The different values of the slopes of the line segments making up the trajectories approaching P at an angle 6Q 119 v i List of Figures 2.1 A pictorial description of the 3+1 representation of a generic spacetime (M.,gab). na is a unit normal vector to the Cauchy surface Et. The vector t a , which generates the flow of "time", is decomposed into the lapse N and the shift Sa. 13 2.2 The spacetime interval between two points with coordinates (x\ t) and (V + dx% t + dt) can be expressed in terms of the lapse N and shift Sl by ds2 = -N2dt2 + hi^S'dt-r dx^Sidt + dxj). In the diagram, rj is a unit normal vector to the surface of constant t 14 5.1 The molecules contained in the box are initially confined to side A by a . partition. In my simplified model, the molecules may move in any one of eight directions as indicated by the arrows 59 5.2 The frequency distribution of return times over 1000 runs of the simulation for a box with 5 molecules . . 65 5.3 The frequency distribution of return times over 1000 runs of the simulation for a box with 10 molecules. Data is not shown for 7 of these runs for which the maximum return time of 1000000 iterations was exceeded. 65 5.4 The frequency distribution of return times over 1000 runs of the simulation for a box with 15 molecules. Data is not shown for 43 of these runs for which the maximum return time of 1000000 iterations was exceeded. . . . 66 vn 5.5 The frequency distribution of return times over 1000 runs of the simulation for a box with 20 molecules. Data is not shown for 424 of these runs for which the maximum return time of 1000000 iterations was exceeded. . . . 66 5.6 The frequency distribution of return times over 1000 runs of the simulation for a box with 25 molecules. Data is not shown for 944 of these runs for which the maximum return time of 1000000 iterations was exceeded. . . . 67 5.7 The frequency distribution of return times over 1000 runs of the simulation for a box with 30 molecules. Data is not shown for 998 of these runs for which the maximum return time of 1000000 iterations was exceeded. . . . 67 5.8 Circle Sl represented as the unit interval with endpoints identified . . . . 72 5.9 The difference between the frequency of occurrence of the digit 1 in the sequence Sn calculated using the program S E Q U E N C E . C and the actual probability of occurrence given by equation ( 5.40 ). This was calculated for the first n elements of Sn, with n ranging from 0 to one billion. The diamonds represent the data on the actual scale, while the other symbols represent the same data on smaller scales (e.g. data shown as stars, on the 102 scale, should be read as the axis value times 10 - 2 ) 79 5.10 The difference between the frequency of occurrence of the digit 2 in the sequence Sn calculated using the program S E Q U E N C E . C and the actual probability of occurrence given by equation ( 5.40 ). This was calculated for the first n elements of Sn, with n ranging from 0 to one billion. The diamonds represent the data on the actual scale, while the other symbols represent the same data on smaller scales (e.g. data shown as stars, on the 102 scale, should be read as the axis value times 10 - 2 ) . 80 Vlll 5.11 The difference between the frequency of occurrence of the digit 3 in the sequence Sn calculated using the program S E Q U E N C E . C and the actual probability of occurrence given by equation ( 5.40 ). This was calculated for the first n elements of Sn, with n ranging from 0 to one billion. The diamonds represent the data on the actual scale, while the other symbols represent the same data on smaller scales (e.g. data shown as stars, on ' the 102 scale, should be read as the axis value times 10 - 2 ) 5.12 The difference between the frequency of occurrence of the digit 4 in the sequence Sn calculated using the program S E Q U E N C E . C and the actual probability of occurrence given by equation ( 5.40 ). This was calculated for the first n elements of Sn, with n ranging from 0 to one billion. The diamonds represent the data on the actual scale, while the other symbols represent the same data on smaller scales (e.g.' data shown as stars, on the 102 scale, should be read as the axis value times 10 - 2 ) 5.13 The difference between the frequency of occurrence of-the digit 5 in the sequence Sn calculated using the program S E Q U E N C E . C and the actual probability of occurrence given by equation ( 5.40 ). This was calculated for the first n elements of Sn/ with.n ranging from 0 to one billion. The diamonds represent the data on the actual scale, while the other symbols represent the same data on smaller scales (e.g. data shown as stars, on the 102 scale, should be read as the axis value times 10~2) ix 5.14 The difference between the frequency of occurrence of the digit 6 in the sequence Sn calculated using the program S E Q U E N C E . C and the actual probability of occurrence given by equation ( 5.40 ). This was calculated for the first n elements of Sn, with n ranging from 0 to one billion. The diamonds represent the data on the actual scale, while the other symbols represent the same data on smaller scales (e.g. data shown as stars, on the 102 scale, should be read as the axis value times 10 - 2 ) 5.15 The difference between the frequency of occurrence of the digit 7 in the sequence Sn calculated using the program S E Q U E N C E . C and the actual probability of occurrence given by equation ( 5.40 ). This was calculated for the first n elements of Sn, with n ranging from 0 to one billion. The diamonds represent the data on the actual scale, while the other symbols represent the same data on smaller scales (e.g. data shown as stars, on the 102 scale, should be read as the axis value times 10~2) 5.16 The difference between the frequency of occurrence of the digit 8 in the sequence Sn calculated using the program S E Q U E N C E . C and the actual probability of occurrence given by equation ( 5.40 ). This was calculated for the first n elements of Sn, with n ranging from 0 to one billion. The diamonds represent the data on the actual scale, while the other symbols represent the same data on smaller scales (e.g. data shown as stars, on the 102 scale, should be read as the axis value times 10~2). . . . . . . . . 5.1.7 The difference between the frequency of occurrence of the digit 9 in the sequence Sn calculated using the program S E Q U E N C E . C and the actual probability of occurrence given by equation ( 5.40 ). This was calculated for the first n elements of Sn, with n ranging from 0 to one billion. The diamonds represent the data on the actual scale, while the other symbols represent the same data on smaller scales (e.g. data shown as stars, on the 102 scale, should be read as the axis value times 10~2) 5.18 The difference between the frequency of occurrence of the digit 1 in the sequence Sn calculated using the program S E Q U E N C E . C and the actual probability of occurrence given by equation ( 5.40 ). This was calculated for one million elements of Sn starting with element 1000. The diamonds represent the data oh the actual scale, while the other symbols represent the same data on smaller scales (e.g. data shown as stars, on the 102 scale, should be read as the axis value times 10 - 2 ) . 5.19 The difference between the frequency of occurrence of the digit 2 in the sequence Sn calculated using the program S E Q U E N C E . C and the actual probability of occurrence given by equation ( 5.40 ). This was calculated for one million elements of Sn starting with element 1000. The diamonds represent the data on the actual scale, while the other symbols represent the same data on smaller scales (e.g. data shown as diamonds, on the 102 scale, should be read as the axis value times 10 - 2 ) . x i 5.20 The difference between the frequency of occurrence of the digit 3 in the sequence Sn calculated using the program S E Q U E N C E . C and the actual probability of occurrence given by equation ( 5.40 ). This was calculated for one million elements of Sn starting with element 1000. The diamonds represent the data on the actual scale, while the other symbols represent the same data on smaller scales (e.g. data shown as diamonds, on the 102 scale, shouldbe read as the axis value times 10 - 2 ) 5.21 The difference between the frequency of occurrence of the digit 4 in the sequence Sn calculated using the program S E Q U E N C E . C and the actual probability of occurrence given by equation ( 5.40 ). This was calculated for one million elements of Sn starting with element 1000. The diamonds represent the data on the actual scale, while the other symbols represent the same data on smaller scales (e.g. data shown as stars, on the 102 scale, should be read as the axis value times 10~2). 5.22 The difference between the frequency of occurrence of the digit 5 in the sequence Sn calculated using the program S E Q U E N C E . C and the actual probability of occurrence given by equation ( 5.40 ). This was calculated for one million elements of Sn starting with element 1000. The diamonds represent the data on the actual scale, while the other symbols represent the same data on smaller scales (e.g. data shown as stars, on the 102 scale, should be read as the axis value times 10 - 2 ) xu 5.23 The difference between the frequency of occurrence of the digit 6 in the sequence Sn calculated using the program S E Q U E N C E . C and the actual probability of occurrence given by equation ( 5.40 ). This was calculated for one million elements of Sn starting with element 1000. The diamonds represent the data on the actual scale, while the other symbols represent the same data on smaller scales (e.g. data shown as stars, on the 102 scale, should be read as the axis value times 10~2). .• 5.24 The difference between the frequency of occurrence of the digit 7 in the sequence Sn calculated using the program S E Q U E N C E . C and the actual probability of occurrence given by equation ( 5.40 ). This was calculated for one million elements of Sn starting with element 1000. The diamonds represent the data on the actual scale, while.the other symbols represent the same data on smaller scales (e.g. data shown as stars, on the 102 scale, should be read as the axis value times 10 - 2 ) . ; . . 5.25 The difference between the frequency of occurrence of the digit 8 in the sequence Sn calculated using the program S E Q U E N C E . C and the actual probability of occurrence given by equation ( 5.40 ). This was calculated for one million elements of Sn starting with element 1000. The diamonds represent the data on the actual scale, while the other symbols represent the same data on smaller scales (e.g. data shown as stars, on the 102 scale, should be read as the axis value times 10 - 2 ) xm 5.26 The difference between the frequency of occurrence of the digit 9 in the sequence Sn calculated using the program S E Q U E N C E . C and the actual probability of occurrence given by equation ( 5.40 ). This was calculated for one million elements of Sn starting with element 1000. The diamonds represent the data on the actual scale, while the other symbols represent the same data on smaller scales (e.g. data shown as stars, on the 102 scale, should be read as the axis value times 10 - 2 ) 96 6.1 A 2-torus is constructed from R 2 by identifying the two sides labeled a and also the two sides labeled b. The coordinate q1 is made L a-periodic and q2 is made Z/{,-periodic (La and Lb are the lengths of sides a and b, respectively) by defining an equaivalence relation on R 2 101 6.2 Identification of points in the /?_) plane as a result of the group action of 7 T 3 . The points marked by the same symbol are identified with each other. 109 6.3 The bouncing of a trajectory off one of the walls amounts to a reflection off that wall. It is clear from this figure that a,- = « / I l l 6.4 The two dimensional momentum space W consists of two cones with tips touching at the origin 113 6.5 A radial trajectory 0o — constant in configuration space is reflected at the origin to another radial line given by | — 90 116 6.6 A sample of trajectories in configuration space which bounce off the point P at various angles " 117 6.7 The shaded region represents the set of points Sc in configuration space covered by the evolution of all the trajectories which approach the point P from an angle in the range 70° < 60 < 80° 121 xiv 6.8 The shaded region represents the full set of points Sm in momentum space projected onto the plane corresponding to the slopes of all the line segments comprising the trajectories which approach the point P from an angle in the range 70° < 80 < 80°. The ra,-. are the slopes of the indicated . lines. . . .. . . . . . . . . . . . . ; 122 6.9 If we construct an invariant set .starting with: the region Sc x Sm, then since it contains a non-radial trajectory, such as the one represented by the solid lines 1-4, it will contain the whole configuration space 124 6.10 If we construct an invariant set starting with a product set such as the region Sc x Sm, then if the set contains any radial trajectory other than 6 = | (line 4), such as the one represented by lines 1 and 2, it will contain the whole configuration space. . . . . . 125 6.11 The evolution of the f3+(3_ wedge, represented by the shaded region, for one particular momentum value yields four strips. This is shown for the maximum and minimum values of the ratio p^/p+. These place.bounds on the slopes of the strips corresponding to all the other values of the ratio • p~/p+. in the invariant set Y. 127 xv Acknowledgements I would like to gratefully acknowledge the tremendous help and guidance provided to me by my supervisor, Dr. Kristin Schleich, during the course of this research. I would also like to give special thanks to Dr. Don Witt for teaching and guiding me through much of the background material and for his many useful comments. It has truly been a pleasure to work with both of you. You have helped me to realize the value of being rigourous in one's research and of questioning assumptions and results that others may take for granted. For this and for helping to make my experience at U B C a pleasant and rewarding one, I am in your debt. I must offer my most heartfelt thanks and everlasting appreciation to my parents for always and unfailingly supporting me through my ups and downs. Thank you also to my brother Vipul for all your advice and encouragement. I love you all very dearly. I appreciate the helpful comments made Dr. Frank Curzon, and I thank him for taking the time to carefully read this thesis. I also thank Dr. B i l l Unruh for reading this thesis and for all his comments. I give my compliments to the physics department at U B C . It is really a warm, open, and friendly place to learn, study, and teach. This is thanks in large part to the faculty, and staff of the department. I especially thank Dr. Malcolm McMillan for helping me to get off ou the right foot in my graduate work. My experience here would not have been as enjoyable and stimulating as it has been were it not for my fellow graduate students in the department - we came together in physics; we stay together in friendship. Particular thanks go to Lori and Tiago for the interesting and always helpful discussions we had. Finally, I would like to thanks Calhoun's at 3035 West Broadway for staying open all night and providing a place with an excellent atmosphere and great mochas - a place where I did a substantial portion of this work! xvi Chapter 1 Introduction Gravity remains the most mysterious of the known forces of nature. Newton's law of gravitation does a very good job of describing how bodies interact under the influence of gravity, but it does not tell us anything about what causes the interaction. It is an empirical law. In 1915, Albert Einstein introduced his theory of gravitation, also known as the theory of general relativity. It succeeded in resolving certain discrepancies which had arisen between Newton's theory and observations, and firmly established itself as the "correct" theory of gravitation. Perhaps more importantly though, it provides us with an explanation of what gravity is - a physical manifestation of the curvature of spacetime. This was a great achievement, but it does not answer all our questions by any stretch of the imagination. Researchers have, for nearly three quarters of a century now, been exploring the im-plications and applications of general relativity. Nevertheless, there remain significant difficulties which have yet to be overcome and which limit our ability to usefully exploit Einstein's work. If we wish to use general relativity to construct a model of the universe, we must be able to solve Einstein's equation. Herein lies the source of the problem. Gen-eral relativity is a nonlinear theory containing an infinite number of degrees of freedom. This makes looking for a completely general solution to Einstein's equation a formidable if not impossible task. We are thus forced to seek solutions for special cases in which the degrees of freedom have been reduced to a finite number by imposing some symmetry. Many such solutions have been compiled, for example in [16] and [18]. 1 Chapter 1. Introduction 2 O n c e w e h a v e s u c c e e d e d i n s o l v i n g E i n s t e i n ' s e q u a t i o n f o r a g i v e n s p e c i a l c a s e , w e h a v e a c o s m o l o g i c a l m o d e l w h i c h d e s c r i b e s t h e e v o l u t i o n o f s o m e u n i v e r s e - n o t n e c e s s a r i l y o u r s . N o w , b y w h i c h c r i t e r i a m a y w e d i s c a r d a p a r t i c u l a r s o l u t i o n ? I n o t h e r w o r d s , h o w m a y w e d e t e r m i n e w h e t h e r o r n o t i t is p h y s i c a l l y r e a l i s t i c b a s e d o n o b s e r v a t i o n s o f o u r U n i v e r s e ? T h e p r o b l e m is t h a t e v e n i f t h i s m o d e l u n i v e r s e d o e s n o t i n i t i a l l y l o o k l i k e o u r s , i t m a y e v o l v e t o e v e n t u a l l y l o o k l i k e o u r s . A f t e r a l l , n o o n e k n o w s f o r s u r e w h a t o u r U n i v e r s e l o o k e d l i k e i n i t s i n f a n c y . A s w e l l , w e a r e a b l e t o o b s e r v e o n l y a v e r y m i n i s c u l e p o r t i o n o f a l l t h e s p a c e t i m e p o i n t s o f o u r U n i v e r s e . . M a n y m o d e l s c a n b e r u l e d o u t b e c a u s e w e c a n p r e c i s e l y t r a c e t h e i r e v o l u t i o n a n d s h o w t h a t t h e y w i l l n e v e r e v o l v e t o l o o k a n y t h i n g l i k e o u r U n i v e r s e , f o r a n y c h o i c e o f i n i t i a l c o n d i t i o n s . B u t t h i s is n o t . a l w a y s t h e c a s e . I n o r d e r t o e v o l v e a p a r t i c u l a r m o d e l , a set o f i n i t i a l c o n d i t i o n s m u s t b e s p e c i f i e d . I n m a n y case s , i t m a y b e t h a t t h e e v o l u t i o n o f t h e s y s t e m d e p e n d s s e n s i t i v e l y u p o n t h e i n i t i a l c o n d i t i o n s t h a t a r e c h o s e n . T h u s , i t c a n b e c o m e d i f f i c u l t t o d e t e r m i n e i f s u c h a s y s t e m c a n e v o l v e t o l o o k l i k e o u r U n i v e r s e i f t h e i n i t i a l c o n d i t i o n s a r e c h o s e n c a r e f u l l y e n o u g h . W h e t h e r o r n o t c o s m o l o g i c a l m o d e l s d i s p l a y i n g s u c h s e n s i t i v i t y t o i n i t i a l c o n d i t i o n s o c c u r g e n e r i c a l l y i n t h e s p a c e o f a l l s o l u t i o n s o f E i n s t e i n ' s e q u a t i o n is a n o p e n q u e s t i o n . I n o r d e r t o i m p r e s s u p o n y o u t h e i m p o r t a n c e o f u n d e r s t a n d i n g h o w a n d t o w h a t e x t e n t a c h o i c e o f i n i t i a l c o n d i t i o n s w i l l af fect t h e e v o l u t i o n o f a p h y s i c a l s y s t e m , l e t m e r e c o u n t t o y o u t h e s t o r y o f t h e s o - c a l l e d " b u t t e r f l y e f fect" , a s t o r y w h i c h h a s n o w f a l l e n i n t o t h e r e a l m o f f o l k l o r e f o r m a n y w h o s t u d y c h a o s a n d d y n a m i c a l s y s t e m s 1 . It t e l l s o f h o w a s l i g h t m o d i f i c a t i o n i n a s y s t e m ' s i n i t i a l c o n d i t i o n s c a u s e d d r a s t i c c h a n g e s i n t h e s y s t e m ' s f u t u r e e v o l u t i o n , a n d a l s o i n o u r v i e w o f d e t e r m i n i s m i n p h y s i c s . T h e r e w a s a l o n g h e l d b e l i e f a m o n g s t p r a c t i t i o n e r s o f p h y s i c s t h a t i t is p o s s i b l e t o d e t e r m i n e t h e b e h a v i o u r o f a p h y s i c a l s y s t e m t o a v e r y g o o d a p p r o x i m a t i o n w i t h o u t 1 James Gleick tells a m u c h more detailed version of this story in his book [12]. Chapter 1. Introduction 3 knowing the initial conditions exactly. In other words, a very small deviation in the initial state of the system will produce a correspondingly very small change in its evolution. I can remember a physics lab instructor some years ago telling us that we do not have to account for the effect of Pluto's gravitational interaction with the Earth in our experimental results. Its influence is so negligibly small that it may be ignored. This kind of talk was common, a product of the era of determinism - the great Newtonian paradigm that nature operates according to a fixed set of rules. If there is any phenomena occurring in nature which we cannot predict or do not understand, then it is simply because we have not yet unraveled the underlying fundamental laws governing it or because we lack the computational power to solve the equations. Understanding nature was like solving a very complicated puzzle; once we know how to put all the pieces together, we have the ability to predict how the world around us will evolve. As Richard Feynman once put it, "Physicists like to think that all you have to do is say, these are the conditions, now what happens next?". And why not? This recipe has worked very well for us for the last few centuries. Sure we can't measure the moon's orbit around the Earth or the course that Halley's Comet is following exactly. But we can measure these things well enough to know with certainty when high tide will occur tomorrow anywhere on the planet and that Halley's Comet will return every 76 years. Meteorologist Edward Lorenz must have been working under these same preconcep-tions when in 1960 he set up his simplified model atmosphere running on his new and, by today's standards, primitive Royal McBee computer at MIT. The weather in his little digital universe followed a precise set of mathematical rules - a set of three coupled differ-ential equations which are known today as the Lorenz system. These equations contain a number of parameters which Lorenz could fix as he wished, and then by specifying some initial values he set his universe in motion. One day in 1961, he decided that he wanted to examine a particular run in greater detail. But instead of starting the entire Chapter 1. Introduction 4 run over again from the beginning, he took a shortcut. He restarted the previous run part way through by typing in the initial conditions straight from the earlier printout. Then he walked away to get a cup of coffee and to let his computer crunch away. When he returned an hour later, he saw something very unexpected. There was no doubt in his mind that this run should have duplicated the previous one exactly. After all, the numerical rules governing the evolution were the same, and so were the initial conditions that he typed in. And yet what he saw on the printout was that the weather pattern in this latest run rapidly diverged from the previous run's pattern and soon the two patterns bore no resemblance at all to each other. His first thought was that his computer may have malfunctioned, but he soon realized the cause of the discrepancy. His computer's memory stored six decimal places, but only three appeared on his printout. So Lorenz had entered the initial conditions rounded off to three decimal places and believed that the difference, being less than one part in a thousand, would have no consequence on , the computer's output. Any one of his colleagues at the time would have likely agreed that this was a perfectly reasonable assumption. Such a small deviation in the initial conditions could be compared to the change caused by a butterfly flapping its wings. Surely this should have no effect on the large-scale features of next month's weather! And yet in Lorenz's model - a system much less complex than the real atmosphere - it had a very dramatic effect. The lesson to be learned from this story is that knowing what the phase space portrait of a system looks like is very important when determining the evolution of the system. If the system exhibits seemingly random behaviour as did Lorenz's system (behaviour which we characterize as chaotic, a term which I will discuss later), then the predictive power of that model is severely limited. It becomes impossible to accurately predict the state of the system at some arbitrary time in the future (or past) even if the initial conditions and evolution equations are known. Uncertainty in the initial conditions and Chapter 1. Introduction 5 round-off errors in the computations will eventually render meaningless the results of numerical simulations. Thus it is important to know whether the trajectories in phase space are generically chaotic or not; otherwise there is in general no way of knowing whether or not we can trust the results of a numerical simulation. If the final state of the physical system is not known beforehand (as is case with weather forecasting, cosmology, and most other things in life), there is nothing against which to check the results of a numerical simulation. In this study, I will examine the issues associated with determining whether or not the evolution of a.given cosmological model is ergodic (a necessary ingredient of chaos). This includes defining and constructing an appropriate phase space on which to study the dynamics of the model. Methods from dynamical systems theory are then applied, taking into account the potentially nontrivial topology of the phase space. I will focus on the homogeneous Bianchi cosmologies, and in particular on the Bianchi type I model as a concrete example. The homogeneous cosmologies are a natural class of models to choose as the basis for this study, for a number of reasons. The most obvious reason is the one that was suggested earlier, namely that it is a class of solutions that we can study analytically; it has a finite number of degrees of freedom. There are, however, other considerations which provide a strong motivation for studying this special class of solutions. The distribution of galaxies on very large scales appears to be homogeneous. Also, a homogeneous Universe is consistent with observations of the cosmic microwave background radiation [32]. While this does not conclusively establish that we live in a homogeneous Universe, it strongly suggests that homogeneous cosmologies may at least be useful models which are good approximations. A case in point is the homogeneous and isotropic Friedman-Robertson-Walker model, which has been the basis of most of the studies carried out by physical cosmologists. Finally, a homogeneous Universe is in accord with the Copernican notion Chapter 1. Introduction 6 that we do not occupy a privileged place in our Universe. In other words, if our location in the. Universe were changed, our surroundings would still look essentially the same. While this is not in itself a valid argument in favour of homogeneous cosmologies, it may sit well with many people's prejudices. Much work, both numerical and analytical, has already been done on various aspects of the homogeneous solutions of the Einstein equatons. Many people have carried out numerical studies of the evolution of certain homogeneous models (for example, [4], [5], [6], [13], [15], [23], [25], [27], and others). My goal is to demonstrate that such numer-ical work alone is not sufficient to" provide a complete understanding of the dynamical behaviour of these models. I will show how an analysis using analytic tools from dynam-ical systems theory can serve to complement numerical studies by providing additional information about the evolution of these systems. Such information can be essential in order to correctly perform and interpret the results of numerical studies. Chapter 2 General Relativity as a Hamiltonian System The vacuum Einstein field equation Gab = 0 (2.1) expresses the complete dynamics of general relativity in an empty spacetime, and in many instances it is convenient to work directly with this equation. However, it is often more useful to work with a Lagrangian or a Hamiltonian formulation. The initial value formulation of general relativity, which I will discuss shortly, is best motivated by the latter. Of more immediate interest to us though, a Hamiltonian formulation of general relativity casts the theory into a form which is most suitable for applying dynamical systems analysis. . . In this chapter, I will develop the Hamiltonian formulation for the dynamics of general relativity using what is often called the A D M 1 formalism. Before proceeding to this task, I will need to present a considerable amount of terminology and background material. Specifically, I will discuss the causal structure of spacetime and define what a globally hyperbolic spacetime is. This material will be necessary in developing the initial value formulation of general relativity. I will then present a derivation of the action for general relativity in canonical form in the A D M formalism. This will lead into a discussion of the gauge invariance of general relativity under coordinate transformations and the con-struction of a phase space of gravitational degrees of freedom defined by the canonically 1Arnowitt, Deser, and Misner 7 Chapter 2. General Relativity as a Hamiltonian System 8 conjugate variables in the Hamiltonian formulation. 2.1 3+1 Spacetime Decomposi t ion In the Hamiltonian formulation, it is necessary to decompose the spacetime into space and time (a "3+1 decomposition"). Therefore, it is useful to have conditions which will tell us when such a spacetime decomposition is possible. To identify what these conditions are, we must examine the causal structure of spacetime. In this section I will present some terminology which will enable us to define the "domain of dependence" of a hypersurface embedded in a spacetime. This will lead us to a condition which guarantees the existence of a 3+1 spacetime decomposition. A spacetime ( M , ^ ) can possess a well defined causal structure only if it is possible to define, on it the notions of future and past events. In Minkowski space, this is done by means of a lightcone. Since the tangent space Tp .of M at any point p is isomorphic to Minkowski space, we may define the lightcone of M at p to be the lightcone passing through the origin of Tp. We may thus make a pointwise definition of past and future by calling one half of the lightcone "past" and the other half "future" for each p £ M . If such an identification can be made in a continuous way for all p G M (this is not always possible if M is not simply connected),' then M is called time orientable. Since it is essential to have notions of future and. past on M , only time orientable spacetimes will be considered in the following discussion. A future directed timelike or null vector at p G M is one which lies in the future half of the lightcone at p. One which lies in the past half of the lightcone is past directed. There exist smooth non-vanishing timelike vector fields on ( M , gab) if and only if M is time orientable. If it is, then a differentiable curve \{t) is called a future directed timelike Chapter 2. General Relativity as a Hamiltonian System 9 (resp. causal) curve if at each point p £ A(i), the tangent vector ta is a future directed timelike (resp. future directed timelike or null) vector: Past directed timelike and causal curves are defined analogously. The chronological future of p £ M , denoted I+ (p)r, is the set of events that can be reached by a future directed timelike curve starting from p. For any subset S C M , we can define the chronological future of 5" as ' I+(S) = UpeSi+(P) • j (2.2)' If q £ 7+(p), then q there is a future directed timelike curve connecting p and q. Now it is always possible to slightly deform the endpoint q of this curve while preserving its timelike nature. Thus, there exists an open ball B around q such that B C I+(p)-Therefore,'7+(p) is an open set in M for all p and consequently so is I+(S). Analogous statements can be made for the chronological past of a point p, I~(p). The causal future of p £ M , denoted J + ( p ) , i s defined in the same way as I+(p) except that each q £ I+(p) can be connected to p by a future directed causal, rather than timelike, curve. Note that in general, for any p, p £ J+(p) but p ^I+(p). Also, J+(p) may or may not be closed in a general spacetime, but it must be closed in any globally hyperbolic2 spacetime. As we did for 7+(p), we define J+(S) = UpeSJ+(p) , (2.3) and analogous statements can be made about the causal past of any point p, J~{p). A-set S C M is achronal if there are.no two points p and q in S such that q lies in the " chronological future of p. The boundary of the chronological future of S C M , dI+(S), is (if it is non-empty) a 3-dimensional achronal C°-submanifold embedded in M , as long as M is a time-orient able spacetime. 2 This will be defined shortly. Chapter 2. General Relativity as a Hamiltonian System 10 We next define the notion of extendibility of a continuous curve. Let X(t) be a future directed causal curve, p G M is a future endpoint of A if for every neighbourhood B of p, there exists a t0 such that X(t) £ B for all t > t0. The curve A is said to be future inextendible if it has no future endpoint. Past inextendibility is defined similarly. Let S be closed and achronal. We define the edge of S as the set of points p £ S such that every open neighbourhood B of p contains a point q £ I+{p), a point r £ I~(p)-, and a timelike curve A from r to q which not intersect S. If S has no edge, it is called a slice. The future domain of dependence of S, denoted D+(S), is the collection of all points p £ M such that every past inextendible causal curve through p intersects S. The past domain of dependence D~(S) is defined analogously, by interchanging "future" and "past". The full domain of dependence of S is D(S) = D+ (S)+D~(S) . (2.4) D(S) represents the set of all events for which conditions are completely determined by a knowledge of conditions on S. A closed achronal set E for which -D(S) = M is called a Cauchy surface. It follows from the definitions that a Cauchy surface is a slice. A globally hyperbolic spacetime (M,gab) is one which possesses a Cauchy surface S. Since E is achronal, we can think of it as representing an instant in "time" in M. In a glob-ally hyperbolic spacetime, the entire future and past histories can be determined from conditions at any given instant of time represented by the Cauchy surface E t . On the other hand, if the spacetime is not globally hyperbolic then knowledge of conditions on a single Cauchy surface is not enough to determine the full history (past and future) of the spacetime. While it has been argued that all physically realistic spacetimes must be globally hyperbolic, there is nothing in the theory of general relativity that excludes ones Chapter 2. General Relativity as a Hamiltonian System 11 that are not. However, the following theorem 3 asserts that globally hyperbolic space-times permit a 3+1 spacetime decomposition, a necessary condition for the Hamiltonian formulation. Theorem 2.1 Let (M.,gab) be a globally hyperbolic spacetime. Then, a global time func-tion t can be chosen such that each surface S t of constant t is a Cauchy surface. Thus, M can be foliated by Cauchy surfaces and the topology o / M is I x E , where £ denotes any Cauchy surface. . . . 2.2 The In i t i a l Value Formula t ion and A D M Fo rma l i sm In classical physical theories, we can generally specify a set of initial conditions, possibly subject to some set of constraints, and then allow the system to evolve according to a set of evolution equations. Its behaviour thereafter is completely specified by the choice of initial conditions. Such is the case for Newtonian mechanics and electromagnetism. In the case of general relativity we are solving for the spacetime itself rather than for the evolution of some quantity in a given spacetime. Thus, in order to construct an initial value formulation, we must first determine which quantities need to be prescribed initially in order for the spacetime structure to be completely determined. We then view general relativity as describing the time evolution of these quantities. We will consider only globally hyperbolic spacetimes, so that we may foliate (M, gab) by Cauchy surfaces S< which are parametrized by a global time function t. Otherwise, it is not possible to construct an initial value formulation. Let na be the unit normal vector field to the hypersurfaces E 4 , as shown in figure 2.1. The spacetime metric gab induces a 3-dimensional Riemannian metric hab on each Cauchy surface according to Kb = gab + nanb . (2.5) 3 Theorem 8.3.14 (page 208) of [3i]/ Chapter 2. General Relativity as a Hamiltonian System 12 We define a vector field ta on M satisfying t a V a i = 1. Then we can decompose ta into two parts, normal and tangential to £ 4 as follows: N = -tana = (n'Vaty1 • (2.6) . Sa = habtb . (2.7) These components are called the lapse and the shift, respectively. The vector field ta can be viewed as generating the flow of "time" throughout M. To move forward in parameter time t, starting at the surface E 0 and ending up at the surface £ t , we map E 0 into £ 4 by a diffeomorphism. This diffeomorphism results from moving along the integral curves of ta. We may view this mapping as effectively transforming the spatial metric on a 3-dimensional manifold E from hab(0) to hab(t)- Viewed this way, the spacetime represents the development of a Riemannian metric on a some fixed 3-dimensional manifold. This suggests that the spatial metric on a Cauchy surface is a dynamical variable in general relativity. The lapse and the shift are not dynamical variables. They only prescribe how to move forward in parameter time t. Once they have been specified on E 0 , they become thereafter free variables. By analogy with classical mechanics then, we would expect an initial data set to consist of a Riemannian metric hab and its "time" derivative on a Cauchy surface E. The "time" derivative of on £ is a well-defined quantity, called the extrinsic curvature Kab- It is given by Kab = KhdhVcnd (2.8) = ^Lfihab , (2.9) where na is a unit timelike vector field normal to E. The following theorem4 asserts that (E, hab, Kab) forms a complete initial data set for general relativity. 4This theorem follows from the discussion in section 10.2 of [31], starting on page 256. Chapter 2. General Relativity as a Hamiltonian System 13 M + At 1 N n a / t a / X t —=3*- / Figure 2.1: A pictorial description of the 3+1 representation of a generic spacetime (M,</a;,). n a is a unit normal vector to the Cauchy surface S*. The vector r a , which generates the flow of "time", is decomposed into the lapse N and the shift Sa. Theorem 2.2 Let E be a 3-dimensional manifold, hab a Riemannian metric on E , and Kab a symmetric tensor field on E , all subject to certain initial value constraints. Then, there exists a globally hyperbolic spacetime (M,^) satisfying Einstein's equation which possesses a Cauchy surface diffeomorphic to E on which the induced metric is and the induced extrinsic curvature is Kab-Given the above discussion of the initial value formulation, it is very straightforward to summarize the A D M formalism. The spacetime line element may be expressed in terms of an arbitrary set of coordinates {x^} as ds2 = g^dx^dx1. (2.10) Chapter 2. General Relativity as a Hamiltonian System 14 x1 x1 + di Figure 2.2: The .spacetime interval between two points with coordinates (xl, t) and (V + dx1', t + dt) can be expressed in terms of the lapse JV and shift Sl by ds2 — —N2dt2 + hij(Sldt + dxl)(S^dt + dx3). In the diagram, n1 is a unit normal vector to the surface of constant t. A foliation of the spacetime M by spacelike hypersurfaces is parametrized by a function t, and each hypersurface is parametrized by a set of cooordinates x%. The spacetime metric may then be decomposed in terms of the lapse N(t,xk), the shift'S l(t,x k), and the 3-metric hij(t,xk) on the hypersurfaces, as is evident by examining figure 2.2, in the following way: ds2 = -N2dt2 + hi, (STdt + dx*) (S3dt + dxj) . (2.11) By comparing the previous two expressions, we are led to the following expressions relating the lapse and shift to the spacetime metric: <7oo = hij&Si-N2 (2.12) Chapter 2. General Relativity as a Hamiltonian System 15 </u, = h,rV . (2.13) 2.3 H a m i l t o n i a n Formula t ion of Genera l R e l a t i v i t y The action in the Lagrangian formulation can be written as J drR, _ (2.14) I where dr is the covariant volume on the spacetime M. The A D M formalism enables us to transform the Lagrangian formulation into a Hamiltonian formulation, in which the action is written in a canonical form. In this section I will present a derivation of this Hamiltonian formulation. If M is globally hyperbolic, the action can be split into spatial and time parts, as I discussed in the last section: ' l = L d t L d V N R ' (2J5> where dV is the covariant 3-volume on E ( . The "time" vector ta can be given in terms of the lapse and shift by ta = -Nna - Sa . (2.16) Then, the expression for the extrinsic curvature becomes: = -7^{Kb + D(aSb)) , (2.17) where. hab is the Lie derivative of hab with respect to the the time ta and D is the covariant derivative with respect tp hab. '.• The Gauss-Codazzi equation of differential geometry relates the Riemann curvature Rabcd on the full spacetime to the curvature on any embedded hypersurface. In this case Chapter 2. General Relativity as a Hamiltonian System 16 the embedded hypersurface ig a Cauchy surface, and the equation reads as follows: nabcd = Kh(h3X Refgh + KadKbc ~ KacKbd , (2.18) where lZabcd is the Riemann curvature for the induced spatial metric on E t given by ( 2.5 Contracting the Gauss-Codazzi equation gives: bd TZ = (Kh{hlhhd ReJgh) (gac + nanc) (gbd + nbnd) + KadKbchachbd - KacKbdhachb = R + Rbdnbnd + Racnanc + nancnbnd (heah(h9chhd Refgh) + KcbKbc - KccKd = R + 2Rabnanb + KabKab - K2 . (2.19) The 2nd term on the right hand side of ( 2.19 ) can be written as: 2Rabnanb = 2{-K a c K c a + K2 + V a (ncVcna - naVcnc)\ . (2.20) So, R = 11 + KabKab - K2 - 2 V a ( n c V c n a - n a V c n c ) . (2.21) Thus, the Lagrangian density may be written as follows: C = NR = N{n+ KabKab - K2) - 2NVa (ncVcna - naVcnc) . (2.22) Therefore, the action is: 1=1 dt f dV N(K + KabKab - K2) - 2 f dt I dV NVa (ncX7cna - naVcnc) . • (2.23) On the full spacetime M , the last integral in ( 2.23 ) is the integral of a divergence, and so we can write it as a surface integral. This means that we may drop it when M has no boundary as in the case of cosmology, and write the Lagrangian density as C = N (n + KabKab -K2) . (2.24) Chapter 2. General Relativity as a Hamiltonian System 17 This can now be,written in terms of hab and hab: C = N •—hab + D(aSb)^ 2N Kb + D(aSb))— (h \,ab 1 . \ \ 2 Kb .+,D(aSb) 2N (2.25) We will use this to construct the Hamiltonian density. The first step is to find the conjugate momentum, defined by _ a 6 7T = 8C Using the chain rule, Now, • Shab 6C SC 6Kab Shab SKab Shab N \2Kab -2Kha and Thus, SKab = 1 Shab 2N irab = -Kab + Khab . ab ,a£> (2.26) (2.27) (2.28), (2.29) (2.30) Then £ can be rewritten as: C = N (il - KabKah + K2 + 2Kab-Kab - 2K2) N ~ n + 2 \ 2 N (hab + ) 1 < a b ~ 2 \ ~ ^ (kab + h a b R 2N v-»" ' " V>~»>JJ " ~ v 2N -NH - (habKab + D(aSb)Kab - habhabK - habD{aSb)K) = -NH - (hab(Kab - habK) + Da_ [Sb (Kab - habK) • -SbDa (Kab - hahK) f SaDb [Kah - habK)) = -Nri + hab7rab + 2SaDbirab = ^NH + Trabhab - sana Db Sa. (it** -hahK) (2.31) Chapter 2. General Relativity as a Hamiltonian System 18 where . H = - n + KabKab- K2 (2.32) and Ha = -2DbTrab . (2.33) In the fourth line of ( 2.31 ), I have dropped the two terms which contribute only bound-ary terms to ( 2.15 ). The dynamics of the system are completely.embodied in the Hamiltonian ' H(hij,^) = j.d3x(NH + SaHa) . (2.34) The configuration space A4 is the space of all C°° Riemannian metrics on a Cauchy surface E and the corresponding phase space is its cotangent bundle T*A4. The canonically conjugate variables on the phase space hij and -n13 satisfy the Poisson bracket relation: {hij{x)iirkl{y)} = ^63{x-y).- (2.35) If XH '• T*Ai —*• T(T*A4) is the Hamiltonian vector field corresponding to ( 2.34 ), then dH{h,Tt) = Sl(XH(h,w),-) • '(2.36) where the two-form Q is the symplectic structure on T*M. The time evolution of any quantity F which depends on the canonically conjugate variables and ir13 is given by its Poisson bracket with the Hamiltonian: ftF = {F,H} (2.37) where H is given by ( 2.34 ). This is equivalent to the usual Hamilton equations on T*M: ; ' SH , ... SH ,' ha = -r-rr and ^ = - — . 2.38 3 Sir* Shi3 v ; Chapter 2. General Relativity as a Hamiltonian System 19 2.4 The Phase Space for General Relativity In this section I will give a general discussion of the phase space of gravitational degrees of freedom for general relativity. This will set the framework for constructing the phase spaces for the Bianchi models. I will proceed to carry out this construction for Bianchi I and discuss some important features of its phase space in chapter 6. I mentioned earlier that the lapse and shift are not dynamical variables. In other words, they do not enter into the equations of motion in any way. It thus follows from ( 2.12 ) and ( 2.13 ) that goo and gpi can also be chosen arbitrarily and that all the non-trivial degrees of freedom are contained in the canonically conjugate variables and TT'-7, which parametrize the full phase space. The physical system,though is confined to the region of phase space delimited by the scalar energy constraint . H{hij, =0 (2.39) and the vector momentum constraint rii{hij, 7TL J) = 0 . (2.40) An important remark to make here is that the Hamiltonian ( 2.34 ) is just a combination of the constraints. This is unlike other standard physical theories where the Hamiltonian is the combination of both a constraint and a dynamical part. The labels "energy" and "momentum" for the constraints are used by analogy with classical mechanics, though they do not truly represent physical energy or momentum. The fact that there exist these four constraints on the phase space suggests that there are really at most four gravitational degrees of freedom per spacetime point, even though there are six canonical variables on T*Ai. For any choice of the lapse and the shift on an initial Cauchy surface, the constraints will be maintained by the evolution equations. Generically, the constraint Chapter 2. General Relativity as a Hamiltonian System 20 set' C will be a submanifold of T*M 5. General relativity is invariant under the infinite dimensional group of diffeomorphisms Z>(M) on the spacetime M . This is equivalent to saying that the theory is independent of the choice of coordinate system. This group of diffeomorphisms will act on the phase space, suggesting that the constraint set C is not truly the space of gravitational degrees of freedom. There is an invariance under diffeomorphisms which will be removed by picking a gauge (i.e. by choosing a coordinate system to work in). Let £ denote the set of all spacetime metrics gab which are solutions of the vacuum Einstein equation: £ = {gab | Gab{gab) = 0} (2.41) and. such that ( M , gab) is a maximal development of the Cauchy data on some Cauchy hypersurface. Then, the space of gravitational degrees of freedom is: G = S/V . (2.42) Fischer and Marsden [10] state that near metrics gab in £ with no isometries, Q is a smooth symplectic manifold. Singularities in Q occur in spacetimes with symmetries, and these are of a "conical" nature. In other words, Q is generically a smooth manifold. Unfortunately, we do not know how to find generic solutions of Einstein's equation. As I pointed out earlier, we must impose some symmetry in order to be able to solve Einstein's equation. It is precisely in such non-generic cases that we will expect to find singularities in Q. Indeed, as .we will see in chapter 6, this is the case for the homogeneous cosmologies. As a result, the gravitational phase spaces for the Bianchi cosmologies will be shown to have non-trivial topologies. Let us now examine how the diffeomorphism group V affects the constraint set C. We begin with a globally hyperbolic spacetime M with spacetime metric g which we write 5Some technical conditions must be satisfied for this to be true always. See [10] for more details. Chapter 2. General Relativity as a Hamiltonian System 21 in terms of- some as yet unspecified set of coordinates {x^}. Let S be a compact oriented three-dimensional manifold and i : S M (2.43) an embedding of S such that the embedded manifold i(S) — £ is spacelike, and the pull-back i*(g). = h is a Riemannian metric on S. Let 0(5*; M , ^ ) be the space of all such spacelike embeddings. Let E 0 = io(S) be some Cauchy slice with initial data (<70)fl"o) which has (M,g) as its maximal development. Then i G C maps (go,TTo) to the ( ( 7 , 7 r ) induced on the hypersurface E = i(S).- The set of all of these ( < 7 , 7 r ) defines an orbit in C. Each of these disjoint orbits develops the same spacetime M and thus they are physically equivalent. The orbits define an equivalence relation ~ in C, and two equivalent orbits, belong to the same equivalence class. Choosing one from amongst all these equivalent orbits corresponds to fixing gauge. Two elements belonging to the same equivalence class describe two distinct choices of initial three-geometries whose maximal developments yield the same spacetime M . In the A D M formalism the spacetime line element is written in the manifestly 3+1 form ( 2.11 ) in an arbitrary basis (t,xl). Fixing gauge corresponds to choosing a lapse and shift. This removes all the gauge freedom associated with making coordinate trans-formations on M and fixes the form of the spacetime metric. Now once the gauge has been fixed, there may yet be some residual gauge freedom remaining which must be removed. A l l the gauge transformations involving the t co-ordinate have been accounted for. So if there is any such freedom remaining it must correspond to isometries of the hypersurface metric hai, since, as we can see from ( 2.11 ), the spacetime metric has been completely specified in all other respects by the choice of a lapse and a shift vector. After this residual gauge freedom has been dealt with, we are left with the "true" space of gravitational degrees of freedom; each point in this space Chapter 2. General Relativity as a Hamiltonian System 22 corresponds to a physically distinct spacetime. Chapter 3 Homogeneous Cosmologies Imagine yourself sitting in the middle of an ocean or large lake, away from any land mass and with no ships or other objects in sight; you would not be able to, differentiate one location on the surface of the water from another. If you were to start swimming in any direction, the view around you would not change. Every point on the water would appear the same as every other point at any given instant of time. The world around you would look homogeneous. Intuitively, we may think of a homogeneous cosmology in a similar way, by imagining ourselves to be floating in interstellar space. If every point around us in space were to look the same, then we would, say that the Universe appears homogeneous. • Let us now make this intuitive idea more precise. Before stating the exact definition, I will need to present some definitions and background material. Following this, I will proceed to describe a classification scheme of all the homogeneous cosmologies, which gives rise to the so-called "Bianchi cosmologies". 3.1 Background I will begin this section with a brief summary of some important concepts that I will make use of. A complete presentation of the ideas presented here may be found in a number of sources; the reader is referred especially to [31], [29], and [21] for further details. The material presented here is drawn primarily from these sources. A vector field specifies a vector at every point p of a manifold M . At each point 23 Chapter 3. Homogeneous Cosmologies 24 p £ M is a tangent space T p M , and a vector field selects one vector from each space T p M . Every curve in M has a tangent vector at each point. Conversely, for an arbitrary C 1 vector field, there exists a curve whose tangent vector at each point is the vector field; such curves are called integral curves of the vector field. The paths of different integral curves in M do not cross, except possibly at a point where the vector field V vanishes. If M is n - dimensional, then there exists an (n-1) - dimensional family of integral curves for each vector field on M that covers M . Such a "manifold-filling" set of curves is called a congruence of the manifold. A foliation F of class Cr (r > 0) and dimension p (p < n) on M is a decomposition of M into disjoint connected subsets {La}a£A (A is some index set) which are called leaves of the foliation: F = {La}aeA , (3.1) such that each point of M has a neighbourhood U and a system of coordinates (x,y) : U —• R p x R n ~ p such that for each leaf La, the components of U fl La are described by the equations: yi = .constant, ... , yn-p = constant . " (3-2) A one parameter group of diffeomorphisms (j>t is a C°° map from R x M —> M such that for any fixed value of t £ R , <f>t : M —• M is a diffeomorphism and for all s, t € R , the group multiplicative operation is defined by <f>t- <t>s = </>t+s • (3.3) The map <f>t=Q is the identity element of the group. For each point p £ M , the curve 4>t(p) • R —> M is called an orbit of 4>t. It passes through the point p at t = 0. Now we can associate with the group </>t a vector field v by defining v \p to be the tangent to the curve (f>t(p) at t = 0. The one-parameter group of transformations (j>t is called the Chapter 3. Homogeneous Cosmologies 25 flow of the vector field v, which can be though of as the infinitesimal generator of these transformations. A map between manifolds induces a map on tensor fields. If <f> : M —> N is a diffeomorphism between manifolds, then it can be used to map any tensor field on M to one on N , or vice-versa. 4> is called a pull-back because it "pulls back" a function / : N -> K on N to the function / • 4> : M —• R on M . 4> also defines a push-forward map fa which takes tangent vectors at the point p G M to tangent vectors at the point fap) G N as follows: . _ (o>)(/) = r(/-o) (3.4) where v G T p , (j)*v £ T^v\, and / : N —> 1R is any smooth function .on N . In coordinate systems {Xv} at p on M and {Y^} at fap) on N , fa can be represented as a matrix in these coordinate bases. This matrix is the Jacobian of the coordinate transformation (j) can similarly be used to "pull back" dual vectors (or 1-forms) on N at fap) to dual vectors on M at p via the map <^>* : T^pj —> T* which is defined as follows 1 : (<f>«u)ava = ua (<j)*v)a (3.6) for all va G Tv. <f>* and c^ , can be extended to mappings of tensors of arbitrary rank. If R is a tensor field on M and <j>t a one-parameter group of diffeomorphisms generated by the vector field va such that fyR = R, then fa is a symmetry transformation for R. A symmetry transformation for the metric gai is called an isometry, and the vector field which generates it is called a Killing vector field. The set of all diffeomorphisms of M is an infinite dimensional group. The multiplicative operation of this group is functional composition. The set of isometries of M with metric gab is a subgroup of this group. 1T* is the cotangent space at p £ M Chapter 3. Homogeneous Cosmologies 26 We can use <f>* to "push forward" any smooth tensor field R by a parameter distance r. We thus define the Lie derivative of R with respect to va at the point p by comparing R with the new tensor field arising from the action oi'<f>-t for small t: { /* r>ai...a„ p«i ...a„ -\ ^ R b ^ t K b ^ j . (3.7) Lv is a linear map between smooth tensor fields of the same rank, and it satisfies the Leibniz rule on outer products of tensors2. Since va is tangent to the integral curves of <f>t, Lv(f) = v{f) f ° r functions / : M —> R . Lie derivatives are used to express the notion of the invariance of a tensor field. R is invariant under va if LVR = 0. Those special vector fields va under which R is invariant will be important for cases in which R represents a physical quantity. An example of particular interest is a Kill ing field, for which Lvg = 0 where g. is the metric tensor. A Lie group G is a smooth manifold which is a group and for which the operations of group multiplication and inversion are smooth-. The elements of an m-dimensional Lie group can be locally characterized by m distinct parameters, with the group operations depending smoothly upon these parameters. We can think of this as a generalization of a one-parameter group of diffeomorphisms. The group of isometries of M is always a Lie group. Leading up to the formal definition of a homogeneous cosmology, we will need to be familiar with some properties of Lie groups. The map Lh '• G —»• G defined by ' Lh(g) = hg , (3.8) called left translation by h, is a diffeomorphism. This follows from the smoothness of the multiplication and inversion operations on G. Let LI denote the pull-back map induced by 4>h- If a vector field va on G satisfies L*hva = va (3.9) 2Lv(Ri ® R2) = LVR\ Cg) R2 -{• Ri ® LVR2 for any two tensor fields Ri and R2 Chapter 3. Homogeneous Cosmologies 27 for all h £ G, then va is called left invariant. The left invariant vector fields form a vector space. A left invariant vector field is determined by its value at T e , the tangent space at the group identity element e, since if va is left invariant, va\h=Ll[va |e] • (3.10) Conversely, if we define a vector field in terms of its value at e by the equation above, we produce a left invariant vector field. Thus, the left invariant vector fields on G form an m-dimensional vector space. If 4> is any diffeomorphism on any manifold and va and wa are any two vector fields, we'have: f([v,w])=lP(v),FW] > (3.11) where [, ] denotes the commutator. So if va and wa are left invariant vector fields on a Lie group, the commutator [u,io}0 is also a left invariant vector field. Since the commutator depends linearly on va and wa, this implies that there exists a left invariant tensor field Cbac such that [v,w]a = Cbacvhwc . (3.12) This tensor field Cbc is called the structure constant tensor of the Lie group. It follows from'the definition that Cl = -C«cb. (3.13) Also, the Jacobi identity for commutators: [u,[v1w}} + [v,[w,u}] + [w,[u,v}} = 0 (3.14) implies that Q « C & ] = 0 . (3.15) A finite dimensional vector space with a tensor satisfying ( 3.13 ) and ( 3.14 ) is called a Lie algebra. The left invariant vector fields of Lie group form a Lie algebra. So, every Chapter 3. Homogeneous Cosmologies 28 Lie group gives rise to a Lie algebra. For example, if G is the Lie group of isometries of a manifold M with metric tensor g, then the Lie algebra of G is that of the Kill ing vector fields of g. Also, every Lie algebra gives rise to a unique simply connected3 Lie group. This fact makes the analysis of Lie groups much simpler because it is a lot easier to deal with a Lie algebra than a Lie group. A left invariant dual vector field aa also satisfies a commutation relation. If va is a left invariant vector field, then aava is a constant. So, 0 = V f c (o-ava) = vaVbaa + aaVbva (3.16) where V(, is any derivative operator on G . From this it follows that 2V[aab] = -acCcab. . (3.17) Right invariant tensor fields are defined analogously to left invariant ones via the right translation map Rh(g)=gh. . (3.18) Right invariant vector fields satisfy the same commutation relations ( 3.11 ) except for a change in sign. The action of a Lie group G on the smooth manifold M is a smooth mapping <j> : G x M' —>• M for which the following conditions hold 4: (i) (f>(e,x) = x V.T £ Ai where e is the identity element of G . (ii) <t> (g, <f>(h, x)) = (j>(gh, x) V<?, / i g G and Va; 6 A l . The action of G on M is transitive on M if for any two points p and q of M , there is some element g of G for which g(p) = q. If for all p and q in M there exists a unique element G such that-g(p) = q, then G said to act simply transitively on M . This implies 3 A topological space U is simply connected if its fundamental group is trivial (i.e. every loop in U is contractible to a point). 4This is actually a left action, strictly speaking (see [1], page 261). Chapter 3. Homogeneous Cosmologies 29 that dim(G) = dim(M). The manifold M is homogeneous if its isometry group G acts transitively on it. This means that the metric properties at each point in the space are the same. If, in addition, there are elements of G which leave some point p of M fixed, then the product of any two such elements also leaves p fixed. Since the identity e G G is one of these elements, they form a subgroup of G called the isotropy group of p. The action of each of the elements in this subgroup corresponds to a rotation about an axis through p. M is said to be isotropic about p if its isotropy group is S0(3). If M is isotropic about every point p, it is said to be isotropic. Defini t ion 3.1 ((Spatial ly) Homogeneous Cosmology) A cosmological model is a homogeneous cosmology if it has a foliation of homogeneous spacelike hypersurfaces, and in addition all fields (scalar, vector, tensor, etc.) which are defined on these hypersurfaces are invariant under the action of the isometry group corresponding to the metric on these hypersurfaces. I should point out that in many references, "spatially" is used in the above definition to differentiate it from a cosmology in which the full 4 - dimensional spacetime is homoge-neous. In the remainder if this document, "homogeneous cosmology" will refer to one which is spatially homogeneous. If the hypersurfaces in the above definition are isotropic and the fields are invariant under the action of the isotropy group, then the model is called a (spatially) isotropic cosmology. 3.2 The Bianch i Classification The first part of this chapter's mission has been accomplished. We now know in precise terms what homogeneous cosmologies are and how to characterize them. The next task is to actually construct these models. To this end, we need to know more about the Chapter 3. Homogeneous Cosmologies 30 group of isometries G . Thus, we shall now examine more closely the nature of the group G . Recall from the definition given above that a homogeneous cosmology (M., gab) has a foliation of spacelike hypersurfaces £ t on which the isonietry group G acts transitively. I will assume that the action of G on each E t is simply transitive. This assumption results in only a small loss of generality. It turns out that the only case in which the action of the isometry group is not simply transitive is the Kantowski-Sachs model [31], which I will not consider. By restricting attention to simply transitive action only, the elements of G and the points of Ej can be put into a one-to-one correspondence by choosing any point p G E t and identifying each g G G with g(p). Thus, the action of each g G G corresponds to a left translation on G by g. The tensor fields on E t which are preserved under the isometries, such as the spatial hypersurface metric hab, correspond exactly to the left invariant tensor fields on G . Therefore, the vector fields on E ( which are preserved under the isometries satisfy the commutation relations ( 3.11 ). Similarly, the invariant dual vector fields satisfy ( 3.17 ). Also, the Killing vector fields of E t correspond to the right invariant vector fields of G , and so they satisfy the commutation relation ( 3.11 ) with a change in sign. Now let's see if we can exploit this information to put the metric of a homogeneous spacetime into a convenient canonical form when G acts simply transitively. We need to construct a coordinate system for the hypersurfaces E t . It would seem that the ideal approach is to choose spatial coordinates which are adapted to the Killing vector fields. The components of the spacetime metric would then be independent of the spatial coordi-nates. Unfortunately, this won't work because the coordinate vector fields must commute; but as I mentioned above, the Kill ing vector fields obey the commutation relation (3.12 ) up to a minus sign and so they do not commute unless Cbc = 0. So the best option is to choose a basis of dual vector fields (crl)a, (&2)a, and (cr3) a which are invariant under the Chapter 3. Homogeneous Cosmologies 31 isometries G . It is always possible to find such a basis on any spatial hypersurface, say Eo, by choosing an arbitrary basis at a point p G E 0 and defining the basis everywhere else on E 0 by left translation. Since the spatial hypersurface metric hab on Eo is left invariant, . . hab = E haf}(aa)a(o^} (3.19) a,P=l where the components of the metric,.hap, are all constant on Eo- Let ta denote the unit normal vector to E 0 at the point p. Then we can define the dual vector fields (cra)a on the whole spacetime in terms of their values on an initial hypersurface E 0 by defining (o-a)a on E t by M « ( 0 = #{* a (0)] o (3.20) where ^denotes the one-parameter group of diffeomorphisms generated by ta. The manifold structure of the spacetime M is R x G as long as G acts simply transitively5. Thus, the full spacetime metric gab takes the form 3 gab = -VatVbt+ haP{aa)a (al3)b , (3.21) a,P=l where (cra)a satisfy the commutation relation ( 3.17 ). We now have a recipe by which we can construct a homogeneous cosmological model with simply transitive action. We choose a three-dimensional Lie group G , a basis of left invariant dual vector fields (o~a)a on G , and the functions hap(t), and define the spacetime metric on M = R x G by ( 3.21 ). Thus, with the exception of the Kantowski-Sachs models for which the the Lie group action is not simply transitive, all the homogeneous cosmologies may be constructed this way by classifying all the three-dimensional Lie algebras. Such a classification was first obtained by Bianchi [7]. He classified all the three-dimensional Lie algebras into nine types. This classification is given in table 4.1 5S.ee [31] section 7.2 for details. Such a decomposition of the spacetime is possible whenever it is globally hyperbolic, as I discussed in chapter 2. Chapter 3. Homogeneous Cosmologies 32 of chapter 4, although the choice of basis vectors for the type VIII model differs from the one Bianchi used. The Bianchi classes are divided into two groups, depending upon whether the trace of the structure constant tensor vanishes. If Cba = 0, then the model is called type A; otherwise, it is called type B . A modified version of Bianchi's classification scheme has also been developed by Ellis and MacCallum [9]. 3.3 Phase Space for the Homogeneous Cosmologies I will now specialize the discussion in section 2.4 about the phase space of gravitational degrees of freedom to spatially homogeneous cosmologies. In this case, the hypersurface metric hab may be written in a manifestly homogeneous form in which its components are completely independent of the spatial coordinates, leaving only a dependence upon t. Further, as I discussed in the previous section, it is possible to construct a global coordinate reference frame on any spatial hypersurface So by picking a local basis of 1-forms at some point p and using left translation to define the basis everywhere else on the hypersurface. This can then be used to construct a globally defined coordinate basis on M by using the, one-parameter group of diffeomorphisms generated by the unit normal vector to E 0 at p to map the coordinate basis on So to the whole spacetime. Thus, the metric can be written in the following form: ds2 = (\2 - SaSa)dt2 + 2Sadaadt'+ hl3(t)da'daj . (3.22) Now by requiring spatial homogeneity, we are restricted to a subset of the set of Riemannian metrics M. On this remaining subset, there still may remain some residual gauge transformations. These transformations must preserve spatial homogeneity; in other words, they must remain within the subset of homogeneous habs. The result gives the actual space of "true" degrees of freedom for the homogeneous cosmologies. This Chapter 3. Homogeneous Cosmologies 33 resulting space is further subdivided by requiring invariance under a particular three-dimensional isometry group. This yields the phase spaces for each of the Bianchi classes. It is these resulting subsets of the full gravitational phase space, from which all the gauge freedom has been removed, which I will examine. Chapter 4 Curvature and Einstein Tensors for the Bianchi Models In the last chapter I described how the spatially homogeneous cosmologies with simply transitive group action have been classified into the nine Bianchi classes. In this chapter I will calculate the components of the Riemann curvature tensor for the Bianchi cos-mologies. I will then use this information to calculate the components of the Riemann curvature tensor for the Bianchi cosmologies. These will be used to determine the Ricci and scalar curvatures and to construct the components of the Einstein tensor. This is largely a lengthy exercise in Cartan calculus. However, as we will see shortly, some com-plications do arise. I will be using a compact notation introduced by Ryan and Shepley [27], and if one carries through the calculations with this notation using the standard method, one obtains incorrect results. Problems arise because certain quantities which appear from the notation to be tensors are in fact not so, and special care must be taken to correctly interpret these quantities. Throughout this document, Roman indices a,b,c,...,i,j,k,... take the values 1,2, or 3, while Greek indices a,/3,7,... take the values 0,1,2, or 3. The Einstein summation convention is also used, unless otherwise noted. In order to carry out a calculation of the curvature tensors, it is necessary to make a definite gauge choice. A convenient choice which I will use for my calculations is where N is the lapse and Sa is the shift which I discussed in the last chapter. This particular choice is called a synchronous gauge. It is defined by the conditions g00 = 1 N = 1 and Sa = 0 , 34 Chapter 4. Curvature and Einstein Tensors for the Bianchi Models 35 and g0a = 0. With this gauge choice, the spacetime line element takes the form ' ds2 ~ dt2 + h^o-'o-' , (4.2) where h{j is a Riemannian 3-metric and {a1} is some as yet unspecified basis of one-forms. Without loss of generality, I can set hij = <-:u (,.•').. . (4.3) where 17 and are arbitrary functions which may in general depend on any of the basis variables. Only nine of these ten functions are independent. In the case of the homogeneous cosmologies, this one-form basis is a, coordinate basis and is constructed as I described in chapter 3. We choose a one-form basis at some point and define the basis everywhere else by left translation. Any basis chosen in this way is an acceptable one. Since {cr4} is a coordinate basis, hij may be chosen to be diagonal. This leaves four arbitrary functions: /?n, / J 2 2 , $33, and O. Only three of these are independent. Thus we may require htj to be traceless. They are functions only of the coordinate t, but not of the spatial coordinates. A particular basis is specified for each of the Bianchi classes by the structure constants, given in table 4.1 Now for convenience in calculations, I will not work with this coordinate basis. Rather, I will use an alternate basis defined as follows: a;0 = dt . (4.4) <J =• <:Q {V') trj : , ' • • (4.5) In this basis, the spacetime line element takes.the form • ds2 =-(u;0)2 + ..(4.6) Chapter 4. Curvature and Einstein Tensors for the Bianchi Models 36 where 6{j is the usual Kronecker delta. The one-forms oj can be expressed in the usual Cartesian basis [27]. They obey the cyclic relation dJ^-C)^ Nak . (4.7) I will adopt the notation /?«=A- • (4-8) Thus, a..' =•. e . - n + ' V (4.9) 4.1 Riemann Tensor Components I will use the method of Cartan 1 to calculate the curvature tensor components. Note that since I am working in the synchronous basis, the spacetime metric in this reference frame is just the flat metric = diag( — l, 1,1,1). Therefore, dg^ = dn^ = 0 (4.10) Thus, it follows from the first Cartan equation dg^ = u v + u>ull (4.H) that the connection two-form u>^ is antisymmetric: ufll/ = (4.12) I now make use of the second Cartan equation dw" = A . w " (4.13) to find these connection two-forms. ^ee for example [21] Chapter 4. Curvature and Einstein Tensors for the Bianchi Models 37 The first step is to calculate the exterior derivatives of the one-forms du° = d(dt) = 0 (4.14) du1 = + ft) t-n+p'u° A cr' + e^'da* = +&)LJ0 ALU'+ e-Q+^C}ka> Aak .= ( - f i +/?•) w° A w»> ^ q f c e n + ^ - ^ - ^ Auk • (4.15) I now adopt two definitions, following the notation used by Ryan and Shepley [27]: Kij; = (-Sl + fySij (4.16) d)k = en+0i-Pi-^qk . . • (4.17) Note that dx-k has the same symmetry properties as C)k\ namely, it obeys the Jacobi relation in three lower indices: 4 A + Vmk + d)md3kl = 0 (4.18) and it is anti-symmetric in the lower two indices: d]k = - 4 j ; - (4.19) With the above definitions, we have dui = KijU0 A LO1 + \d)kLoj A uk . (4.20) It follows from ( 4.12 ) that. cj° = u \ = cu2 = ul = 0 . (4.21) From ( 4.13 ), we find that AJ = -Jv A ui" = -J0 A u° - J, A <J . (4.22) Chapter 4. Curvature and Einstein Tensors for the Bianchi Models 38 Comparison of ( 4.20 ) and ( 4.22 ) yields and -LO\.aJ = ^<ZW A LO (4.23) (4.24) Writing out explicitely the three equations in expression ( 4.24 ), we get -col A to2 — col A LO 1 d^u1 ALO2 + d\3ur ALO3 + d^u2 A LO1. + d\3L02 Aw 3 + d^u3 ALO1 + d\2u3 A LO2 = d^u1 A LO2 + d23u2 ALO3 + d\xL03 A LO1 (4.25) (in the last step, I have made use of ( 4.19 ) and the antisymmetry of the wedge product), and similarly LO2 A LO1 — co2 A LO3 = d\2L0X A LO2 -{7 d\3L02 ALO3 + d231u3 A LO1 (4.26) and - cof A LO1 - u3 A LO2 = [d^to1 ALO2 + d\3L02 A LO3 + d331Lo3 A co1] . (4.27) Now to solve the above system of three equations ( 4.25 ), ( 4.26 ), ( 4.27 ) for the three unknown quantities LO\, LO3, and LO3, let LOo ALO1 + BLO2 + CLO3 . ul = DLO1 + ELO2 + Fco3 cof = GLO1 + Hu2 + Ico3 . Then, inserting these into ( 4.25 ), ( 4.26 ), and ( 4.27 ), we get (4.28) (4.29) (4.30) d\2ul ALO2 + dLu2 A LO3 + dl.LO3 A LO1 - [ALO1 + CLO3) ALO2 + (GLO1 + HLO2) ALO3 = - (BLO2 + CLO3) ALO1- (DLO1 + ELO2) ALO3 =• [d2^1 A LO2 + d223to2 ALO3 + d231u3 A LO1 - [HLO2 + ILO3) A LO1 + (DLO1 + Fco3) Auo2 = [d^LO1 A LO2 + d323L02 ALO3 + d331L03 A LO1] . Chapter 4. Curvature and Einstein Tensors for the Bianchi Models 39 So, identifying terms, it follows that A = ~d\2 (4.31) C + H = d23 (4.32) G = - 4 (4.33) B = -d\2 \ (4.34) C + D = (4.35) E = -d2 u23 (4.36) H + D = d\2 (4.37) I = (4.38) F = ° 2 3 • (4.39) It follows from ( 4.32 ), ( 4.35 ), and ( 4.37 ) that H •= \ (d]3 - 4 + d\2) (4.40) D = | ( - 4 3 + 4 + <4). (4-41) . (• = ~M, -I d^-d\.) . (4.42) Thus, inserting the above values of A through I into ( 4.28 ), ( 4.29 ), and ( 4.30 ), I obtain the following expressions for cuj, and to3: <A = -d\2^ ~ d\2u2 + ^ (4 3 + d231 - dt2) = \(d\k-d\k-d\2)uk (4.43) = \ + <4 + d\2) ux - d223cu2 - d323co3 = \{dlk-d32k-dk23)cok (4.44) Chapter 4. Curvature and Einstein Tensors for the Bianchi Models 40 u3 = - 4 X + l- (d\3 - d231 + d312) u2 - d331u3 = \{d\k - d\k - dk3l) uk . (4.45) The above three expressions may be written in the following compact form:. ^ = \{d)k-d\k-dkJ)uk . (4.46) A l l the components of the connection two-form u^ have now been determined. They are given by expressions ( 4.21 ), ( 4.23 ), and ( 4.46 ). The next step is to calculate the exterior derivative of up*. du°0 = 0 ' ' (4.47) dJu = Ktju° Au3 + K'ijdJ = KijuP A u3 + Ku (K[mu>° A um + l-d!mntom A w") ( using ( 4.20 ) ) = (km + K%lKim)u° Aum+ l-K«dlmnum Aun (4.48) du) = ^ ( 4 - 4 - 4 ) ^ ° A ^ + ^ f e - ^ - 4 ) ^ = \{d)k - 4 - 4 ) , J ) A f \ ~ d\k ~ 4 ) ( W A u' + l-dkmul A c m ) ( using ( 4.20 ) ) " + j ( 4 - d\k - 4 ) Qdfm) J A c m = | [ ( 4 - 4 - 4 ) + ( 4 - 4 - 4 ) * « ] W ° A W ' + \{d)k-d\k-dkJ)dkmul Aum . (4.49) Recall that the metric in my choice of basis is the flat Minkowski metric i]^. Thus, in what follows, I will freely raise and lower indices with no change in the spatial components and an extra factor of -1 in the temporal components. Chapter 4. Curvature and Einstein Tensors for the Bianchi Models 41 The components of the connection two-form can be calculated using the third Cartan equation ' ^ = ^ + U ; A W ? , (4.50) Note that since TV^ is antisymmetric, n°0 = %\ = 1Z22 = ^ 3 = 0 • ( 4 - 5 1 ) It follows from (4.50 ) that 7c° = du° + u°aAu? = dto° + wj A u{ ( since io° = 0 from (4.21) ) = dui + J0/\u{. • . • • (4.52) In the last line, I have raised and lowered some indices and have used.the antisymmetry ( 4.12 ) of io^y. Using ( 4.23 ) and ( 4.46 ), we calculate the wedge product = (K^1) A\(dlk-d)k-dk3)uk = \l<n(d\k-d)k-dk)J Acok . . (4.53) So,, combining ( 4.48 ),( 4.52 ), and (4.53 ) yields 71° = [Kim + K«Klm) u>° A um + \Kjmnum A con + l-K3l (d\k - d)k - dk) J A uk = {Kim + KuKlm) u°Aum+ \^-K«d!mn + l-K3m (dt - d)n - <#)] um A un(4M) (The indices in the last term were relabeled as follows: I —> m, k —» n). Next, we find the VJJ components. It follows from ( 4.50 ) that K) = du)+JaAoo« = dw}+wj AwJ + u4 Acc™. • = duo) + 4 Awj + w i - A w f . (4.55) Chapter 4. Curvature and Einstein Tensors for the Bianchi Models 42 (In the last line, I have again raised and lowered indices and used the antisymmetry ( 4 . 1 2 ) o f u v ) . Using ( 4.23 ) and (.4.46 ), we calculate the following wedge products: = KikKj,ukAul col A tom • fak ~ d?k ~ dl) cufcj A [- {dp, - d?ml - d'mJ) J 7 (<&* - d?k ~ dkm) (dp, - &ml - d'mj) uk A J . (4.56) (4.57) Also, by relabeling the indices in the last term of ( 4.49 ) ( k —»• ra, / —> k, m —»• I), we can rewrite the expression for dco) as H = ^ [ ( 4 - 4 - 4 ) + ( 4 - ^ - 4 ) ^ ° ^ ' '+ \(d)n-dln-d?3)dnkluk Au1 . (4.58) So, combining ( 4.55 ), ( 4.56 ), ( 4,57 )', and ( 4.58 ) yields ^} =' \ [ ( 4 - 4 - 4 ) + {d)k -4k~ 4 ) KH UJ°AU1 + KikKjt + \ (<fmJfc -dfk- dkm) (<% - d?ml - dlmj) + I (d)n - <4 - dl) dnkl u> A to . (4.59) The next step is to use the last of the Cartan equations: = R^a0sua ALO0 ( I a/9 I means a < f3 ) \ r ^ A ^ (4.60) to determine the components of the Riemann curvature tensor R^ap. It follows from ( 4.60 ) that Chapter 4. Curvature and Einstein Tensors for the Bianchi Models 4 3 anc = V A ^ + ^ / A ^ ' - ( 4 , 6 1 ) 1 2' = R)olu° A J + -R)lmJ A um . ( 4 . 6 2 ) Now from ( 4 . 5 4 ) , ( 4 . 5 9 ) , ( 4 . 6 1 ), and ( 4 . 6 2 ), we can pick out the components of R%ap- However, there is a small complication which we must deal with. Note that the second terms in ( 4 . 5 4 ) and ( 4 . 5 9 ) are not antisymmetric in m,n and k, l (respectively), but -R°m n and Rl-ki (resp.) must be antisymmetric in these pairs of indices. Thus, in order to obtain the correct results, we must antisymmetrize these expressions; otherwise, some terms will be missing from the expression. It is important to realize that in doing so, I am not altering the components in question. I am merely ensuring that the notation being used to represent these components is correctly interpreted, and that none of the terms are being inadvertently omitted. This complication does not arise if we write the expressions for the R1- in full detail without any implicit sums and work from there. The moral is that in working with the compact notation of Ryan and Shepley, which I have introduced above, special care must be taken to ensure that the resulting expressions preserve the symmetries of the full expressions, and are not lost as a result of the notation. So, from ( 4 . 5 4 ) and ( 4 . 6 1 ) , we obtain /<;,;,. = Kik-\ Kl,Kjk ( 4 . 6 3 ) and (after antisymmetrizing on the indices k,l) Chapter 4. Curvature and Einstein Tensors for the Bianchi Models 44 = KiAi + ^ >> (4 - 4 - 4 ) - (4 - 4 * - 4 ) j • (4-64) Thus, from (4.63 ) we find: Roo — pa p« ^OaO — -aOiO = —R°;n; -Kit K{j Kji and from ( 4.64 ), we find: •Ok pa pt -"-Oa/c — ^Oik ^iki - K A = K A + h<0k (2d)) K A +-Kjkd)l 4-4)-** ( 4 - 4 - 4 ) using (4.19) ) Now, from ( 4.59 ) and ( 4.62 ), we obtain R)oi - {(4 - 4 - 4 ) t - i ( 4 - 4 . 4 ) K K L and (after antisymmetrizing on the indices l,m) \ {d)k - 4 - 4 ) dL - T fe - 4k - 4 ) <*; (4.65) (4.66) (4.67) + + [ — KimKjl + KnKjm] \ (4m - 4m - C) fe - 4i - 4) + \ fe - 4 - 4 ) fe - 4m - 4 2 (^ 'fc ~~ 4 ~ 4j) 4m — KimKji rf KuKjm + i [- (4m -'4m - C) - 4/ - 4 , ) +.(4/ - 4 - 4 ) (4, - 4m - 4") Chapter 4. Curvature and Einstein Tensors for the Bianchi Models 45 We need to find p pa rtjm — J^jam From ( 4.63 ), we find: = R%m + R]m. (4.69) R]om = Kjm + K3lKlm (4.70) jm and from ( 4.68 ), we find: R)im = \(d)k-dlk-dkJ)dkm-KimK3i + KilK: . - \ (<m ~ din ~ C) [dl - dit - 4 - ) + (d]m - d>nn - d%) (4.71) So, from ( 4.69 ), ( 4.70 ), and ( 4.71 ), we obtain Rjm = kjm + KuK3m + \(d)k-dtk-dk3)dkm ~ \ [dlnm - dln - dZ) {dl - C - dQ + (d]m - dim - dl)) .(4.72) After expanding out and simplifying the last two terms in ( 4.72 ), using Jacobi's relation ( 4.18 ) to rewrite many of the terms, we obtain Rim = K3m + A',,Kjm - (dim + d™) - \d% (<fmk + dkmi) •+ \d\ndZ • .. (4-73) 4.2 E ins te in Tensor Components Now that the components of the Riemann tensor have been determined, I will will use them to determine the components of the Einstein tensor G^ = R^-^g^R (4.74) for the Bianchi cosmologies. I will need the structure constants for the Lie algebras corresponding to each of the Bianchi classes. These, as well as the values of d)k (defined Chapter 4. Curvature and Einstein Tensors for the Bianchi Models 46 in the last section) are given in table 4.1. I mentioned in the previous section that is a diagonal, traceless matrix. Since only two of the parameters /?2, fa are independent, we can rewrite things in terms of the parameters (3+ and /?_, which are defined as follows: P- = 2V3 (fa - AO (4.75) (4.76) or equivalently, fa = (p+ + y/spJ) fa = (p+-Vzp_) fa = -2/3+ . (4.77) (4.78) (4.79) We first work out the scalar curvature. R = R& = Rq + R\ + R\ + R\ = — i?oo + R\\ + R22 + -R33 = —RoO + Rii = \kn + KuKn) + kn + KuKu + \dnnl (4 + 4) - \ d L (4 + 4<) + \d\nd\ = 2kn + A'O Rjt + KuKu + 4 4 - \d\n (4 + 4) + ^ 4 4 = ' -en + i 2 i i 2 + 6 (p\ + /3i) - 4 4 - ^4 (4 + 4) + j 4 4 • We are now ready to begin computing the components of the Einstein tensor. 1 (4.80) 00 (-ki{ - KuK^ + kti + + + 1m% - \d\n (4 + dt) + 1 4 4 1 1 2 v KijKji + KiiKu) 2 L 4 4 + \d\n (dlm + djj) - j 4 4 (4.81) Chapter 4. Curvature and Einstein Tensors for the Bianchi Models 4-7 B ianch i Class Structure Constants Type I 4 = ° A II C23 = — C32 = 1 ^ 3 = - 4 2 = e"e2/3+ 2V3/?- A III Cl3 = ~~ C3 I ' = 1 4 3 = - 4 , = e ne a '+ B IV Cl3 = — = 1 ^ 2 3 = — ' ^ 3 2 = 1 ^ 2 3 ~~ ^ 3 2 — x ^ 3 = _4 2 = ene2/3++2V3f3-4 3 = -42 = e " e ^ B V Cl3 = ~ C 3 1 = 1 ni _ ,02 _ 1 "^23 — ^ 3 2 — 1 C?23 — C?32 = C ^ C 2 ^ " 1 " B VI C\3 = — C g j = 1 C^ = -Ci2 = h(h^ 0,1) d 2 3 = -<i 2 2 = heQe2^ A if h = -1 B : / f c / - l VII ^ 1 3 = — C | l = 1 C23 = — C32 = — 1 C | 3 = - C 3 2 2 = h (h2 < 4) dj3 = ~d231 = e n e 2 / 3 + - 2 V 3 / 3 _ • d\z = -d\2 = -e^e2^+2^-d\3 = -d\2 = heQe2^ A if h = 0 B i / ^ 0 VIII ^ 2 3 = ~~ ^ 3 2 = — 1 C f l = — Cl3 = 1 / ^ 3 _ /^3 _ 1 u 1 2 — — u 2 1 — 1 d\3 = -d\2 = -eneW++2V3P-d231 = -d\3 = e n e 2/3+-2V5/J-d\2 = = e « e - 4 ^ A IX ^ 2 3 ~ — ^ 3 2 = 1 C f l = — Cl3 = 1 ^ 1 2 — ~^21 = 1 d\3 = -d\2 = e " e 2 / ? + + 2 ^ 3 / 3 -d231 = -d\3 = e Q e 2 ^ - 2 ^ -d\2 = -d321 = e"e- 4 ^ A Table 4.1: The Structure Constants C*k and the corresponding values of d%-k for each of the Bianchi Classes. Chapter 4. Curvature and Einstein Tensors for the Bianchi Models 48 The first term in the expression above can be rewritten as follows: l- (••/v0/<;7! + KuKu) = - l - {-tl + A ) {-tl + p3) SijSij + l- {-tl + {-tl + p,) SuSu = -^[dtl2 + p2 + p2 + p23-2tl{p1 + p2-rp3)} + I(-3J7 + ft + p2 + p3f = 3{tt*-p2+ -P2_) . (4.82) Also, we define: 3R = « + \d\n {d\n + d») - id{n4 . (4.83) This is the part of Goo which depends only upon the hypersurface metric 3gi3. Then, our final expression for Goo becomes Goo = 3 {tl2 - P2+ - P2_) - l- 3R . (4.84) The spatial part of this expression, evaluated for each of the Bianchi classes, is given in table 4.2. The results in this table are in agreement with corresponding expressions given in table 11.1 on page 193 of [27], except in the case of Bianchi VIII. Following [27], the expression in the third column of table 4.2 for Bianchi VIII would read: - l-e^ (cosh (4v/3/9_) — l ) - e-2/3+ cosh (2y/3pJ) - (4.85) The reason for this discrepancy is unknown, although having worked through the algebra a number of times barn confident that my expression as given in table 4.2 is correct. Moving on now to the Go; components, Go, = Roi = K3ld!3% + KHd\3 Chapter 4. Curvature and Einstein Tensors for the Bianchi Models 49 Bianch i Class I 0 0 II i d 1 n11 2 U 2 3 U 2 3 _l e 4j3+ e 4v/3/3_ 4 III 2^13^13 IV 2^13^13 + 1^ 23^ 23 + ^d\3d\3 + 2d\3d\\3 -e4 / 3+ (3 + i e 4 ^ - ) V 2d13d13 + 2d\3d\3 + 2d\3d\3 -3e 4 / 3+ V I 2d13d\3 + 2d\3d\3 + 2 ^ 3 ( ^ 2 3 _ e 4/5+ (1 + h + tf) VII 2^13^13 + 2^23^23 + 2d\3d\3 + d31d\2 . -e4^+ (|cosh (4V3 / J - ) -I+-/.2) VIII 1/71 J l , 1 J2 J2 1 J2 11 2 U 2 3 a 2 3 1 2 U 1 3 a 1 3 T a 3 1 a 3 2 - L l / 7 3 J3 , J 2 J 3 1 //I z / 3 ^ 2 12 u12 1 u 1 3 a 1 2 T u 2 3 a 2 1 - |e 4 / 3 + (cosh ( 4 ^ - ) + l ) - e - 2 ^ sinh ( 2 ^ - ) - |e- 8 / J+ IX i/71 Jl _|  1/72 J2 1 J2 J l 2 U 2 3 U 2 3 ^ 2 13 a13 + a 3 1 a 3 2 + 2^12^12 + ^13^12 + ^23^21 - § e 4 / 3 + (cosh (4v /3/3_) - l ) +e~2<3+ cosh (2 v / 3/9-)-- ^ e - 8 ^ Table 4.2: The Spatial Parts of the Goo component of the Einstein tensor for each of the Bianchi Classes, in terms of djk and (ft, /?+, /?_).' The full expression is Goo = 3 ( i i 2 - / 9 2 -P2_)-l3R. Chapter 4. Curvature and Einstein Tensors for the Bianchi Models 50 Bianch i Class G 0 3 Constraint ( G 0 3 = 0) I 0 None II 0 None III (3/3+ + VzpJ) eQe2^ = IV QP+ene2^ p+ = o V 6p+eQe2f3+ p+ = o V I - (-3(1 + h)p+ - V3(l - h)pJ) e n e 2 ^ VII -h (-3/?+ + VSpJ) e P- = Vdp+ifh^o None if h = 0 VIII 0 None i x 0 None Table 4.3: The G03 components of the Einstein Tensor and the corresponding Momentum Constraints for each of the Bianchi Classes. = (-a+A-)^.+ (-ft + /9,)*„4 = (-0 + ^ ) 4 +(-O + A ) 4 = (-& + A)4 • (4-86) Evaluating the above expression, we find that G01 = G 0 2 = 0 (4.87) for all the Bianchi classes. The values of G03 are given in table 4.3, along with the corresponding momentum constraint G 0 3 = 0. Chapter 4. Curvature and Einstein Tensors for the Bianchi Models 51 The next step is to calculate the Gij components. I will derive separate expressions for the diagonal components Gaa and the off-diagonal components Gab(a / b). Please note that in all of the following expressions, there is no implied summation over the index a (although the summation convention continues to apply to all other indices). Also, whenever the indices a and b appear in the same expression, a 7^  b. The diagonal components are given by G a a Baa 2 ^ a a ^ kaa + KuKaa " « " ^4 (4 .+ 4) + J « •1 = 20 - 3ft2 +~Ba- 3flpa - 3 (/32 + 01) " - « ~ \dlan (d'an + dnal) + \dtJl + = 217 - 3 i l 2 + i3a) - mf3a - 3 (/% + (31) + 3 G a a (4.88) where 'Gaa = dldl- U (d}an + •+ J « + l3R (4.89) and 3 R is given by ( 4.83 ). The off-diagonal components are given by Gab = Rab — -8abR -Ui (dfb + dl) -\dlan (d[n + di) + IdfA (4.90) = 3Gab . Table 4.4 lists the spatial parts 3 G a a and 3 G a b of Gij for each of the Bianchi classes. Chapter 4. Curvature and Einstein Tensors for the Bianchi Models 52 B ianch i Class 3r< Q - 2 f t 1 11 ^ n = |e4(/3++v^/3-) 3 Q 3 3 = _le4(!3++V30-) ^ 1 2 = ^ 2 3 = ^ ? 3 1 = 0 III ^ n = 0 ^ = e4^+ ^ 3 3 = 0 . ^ 1 2 = 3G23 = ^ 3 1 = 0 IV 3 G n = e 4 ^ ( l + | e 4 ^ - ) 3 G 2 2 = e 4 ^ ( l - | e 4 ^ - ) 3^ 33 = em ( l - I e ^ - ) 3 G 2 3 = 3 C 3 i = 0 V = e4^+ ^ = e 4 ^ ^33 = e 4 ^ 3Gri2 = 3C 23 — ^Gzx = 0 Chapter 4. • Curvature and Einstein Tensors for the Bianchi Models 53 B ianch i Class v^jj e VI : t7„ = h2e4^ 3G33 = he4^ VII SGu = e 4 ^ ( - | e " 4 ^ - + | e 4 ^ - + h2 - §) ^ 2 2 = e 4 / ? + (_i e4V5/?_ + | g - 4 V 3 / J - _ I ) ^ 3 3 = I e ^ + ( l _ C o S h (4^3/9-)) 3G12 = -he^+e-2^-^J*23 = ^?31 = 0 VIII = £ u = [ e ^ + ( _ I e - ^ - + | e4V5 /J- + I ) _ I e - 8 ^ + + e - 2 / ? + c o s h ( 2 >/3/?-)] 3£ 2 2 = [e4^+ ( _ I e 4 ^ - + | e - 4 v ^ _ + I ) _ 1 e-80+ _ e-20+ c o s h (2^_)] ^ 3 3 = [~\e4^ (cosh (4V3/9-) + l ) + | e - 8 ^ + e " 2 ^ s i n h (2y/38_)] ' ^Cl2 = ^?23 = ^31 = 0 IX V11 = [ e 4 ^ ( - | e - 4 ^ - + | e 4 ^ - - § ) - i e " 8 ^ - e " 2 ^ s inh (2^3/3-)] ^ ? 2 2 = [e^+ ( _ I e 4 v ^ / . - + | e -4^3/3_ _ I ) _ l e - 8 / 3 + + e - 2 / 3 + , s i n h ^ v ^ / ? - ) ] ^ 3 3 = [ - l e 4 ^ (cosh (±y/3B_) - l ) + | e - ^ + - e " 2 ^ cosh (2^3/3-)] ^?12 = ^ 2 3 = ^?31 = 0 Table 4.4: The Spatial Parts ^Gij of the components of the Einstein Tensor for each of the Bianchi Classes. Chapter 5 Dynamical Systems I stated in the introduction that the purpose of this study is to examine the flow on the phase space for the Bianchi models using dynamical systems techniques. In previous chapters, I presented a general discussion of the phase space (of gravitational degrees of freedom) for the homogeneous cosmologies, and I will construct the phase space for the Bianchi I model in the next chapter. Before proceeding to this task though, I feel it would be useful to introduce you to these dynamical systems techniques. Thus my aim in this chapter is to present to you a review of the basic definitions and concepts from dynamical systems theory that I will be using. I will do this more by example than by precise mathematical detail, but all of the detail necessary for the discussion in the next chapter is included and the theorems are stated in their most general forms1. Also, in an attempt to convince you of the power and versatility of dynamical systems techniques, I will apply them to a very simple dynamical system and show how they can lead to unexpected results. 5.1 Concepts and Theorems The setting for a dynamical system is a measure space, which is a special kind of topo-logical space. In order to properly define a measure space, I will first have to present some terminology2. A measure space is equipped with a function called a measure. The 1The theorems presented in this section are as stated in [1], and proofs of these theorems may be found in this reference. 2More details may be found in any standard textbook on real analysis, such as [11]. 54 Chapter 5. Dynamical Systems 55 domain of this function is called a a-algebra, which is defined as follows. Definition 5.1 Let X be a non-empty set. A a-algebra on X is a non-empty collection of sets M C P(X) which is closed under countable unions and complements.3 It is clear from this definition that a a-algebra on X defines a topology on X . Let X be a set with a <r-algebra M. Then, a measure on (X, M) is a function fi : M —> [0, oo] such that the following conditions hold: (1) p(<f>) = 0. 0 0 (2) / / {Ej}^° is a sequence of disjoint sets in M, then ^(U^0 Ej) = ^ ji(Ej). 1 (X, M) is called a measurable space and the sets contained in M are called measurable sets. If we now define a measure a on (X, M), then (X, M, u) is called a measure space. There is one more piece of terminology involving measure spaces which we will need to know. Any mapping / : X —> Y between two sets induces a mapping / _ 1 : P(Y) —> P(X) between their power sets by r1(E) = {x£X\f(x)eE}.. . (5.1) f~x commutes with unions and complements. So, if is a cr-algebra on Y , then {f~l(E) \ E £ A^} is a cr-algebra on X. If ( X , M) and (Y, N) are both measurable spaces, a mapping / : X —> Y. is called a measurable mapping (with respect to M a n d N) if f~l(E) G M for all /•; G .V. The motion of some object in a measure space is described by the flow generated by a vector field. If we imagine our measure space to be a river, then the velocity vectors of the water at each point in the river define a vector field. The integral curves of this vector field, which is a flow on this measure space, describe the motion of the water flowing down the river. 3-P(X) denotes the set of all subsets of X, called the power set of X. Chapter 5. Dynamical Systems 56 Of course, the measure space will usually be more abstract, and the objects in this measure space may not always represent physical objects. But it is often useful never-, theless to think of them as physical objects. Such is the case with the phase space flow of the Bianchi models. The measure, space in this case is a space of functions, and the 'objects' are functions. The flow describes the 'motion' of these functions in the function space, and.the tangent vector fields to this flow describe the 'velocity' with which these functions 'move' in the space. Motion along one of the flow lines describes the evolution of a particular,cosmology.• We have all observed that water flowing down some streams or rivers looks very disordered and chaotic (white water rapids may come to mind) while in others the flow of water appears much more orderly. A similar phenomenon is observed in water coming out of a water faucet. When we turn on the tap, the water usually appears fairly laminar. But as we continue turning the tap to increase the amount of water coming out of the faucet, the flow becomes very disordered. Such flows are often referred to at being chaotic or ergodic, terms which I discuss in greater detail further on. . In order to fre precise about what I mean by the term ergodic, I must first introduce a couple of definitions. A flow on a measure space is called measurable if it is a measurable mapping. Now if every point in some set in a measure space is always carried only to points in this set by the flow mappings, then this set is said to be invariant. Let us go back to the example of water flowing down a river for a moment. If we were to come across a whirlpool in which the water circulates around but never leaves (in other words, water that is in the whirlpool stays in the whirlpool), then this would be an invariant set of the river. The exact definition of an invariant set is as follows. Let S be a measure space and Ft a measurable flow on S. A subset A of S is called invariant if Ft(A) = A V i .€ 1R. .We say that two sets are equal if they are equal almost Chapter 5. Dynamical Systems 57 everywhere4. This now leads to the definition of an ergodic flow on a measure space, which says that there are no measurable proper subsets which are invariant. Def ini t ion 5.2 (Ergodic) A flow Ft on S is called ergodic if the only invariant mea-surable, sets are the empty set <f> and all of S . An equivalent statement of ergodicity can be given in terms of constants of the motion, which are defined as follows. ' ' Def ini t ion '5 .3 A measurable function f: S i—» 1R is called a constant of the motion if and only if f o Ft = f almost everywhere for each t € 1R. As an example, let us think of a rock moving through interstellar space. So long as there are no external forces acting on the rock, it will move in a straight line and its momentum will not change. Thus, the flows on this interstellar measure space for this system are straight lines, and the three components of the momentum are constants of the motion. Now if it so happened that the laws of mechanics were changed so that the total momentum of the rock were not a conserved quantity, then the rock would not necessarily move in a straight line. If fact, we would not know how the rock would move. There would no longer be any constants of the motion. This leads to our second characterization of ergodicity, which is equivalent to the definition above. Theorem' 5^1 A flow Ft on S is ergodic if and only if the only constants of the motion are constant functions almost everywhere. , A symplectic form on a manifold M is a closed5 nondegenerate6 two-form to on M . It defines a symplectic structure on M . A symplectic manifold'(M, CJ) is a manifold M with 4i.e. if they differ by a set of measure zero. 5 A closed forma; satisfies dui = 0. „, 6 A nondegenerate form on a manifold M satisfies: If X is a fixed vector field on M such that u(X, Y) = 0 for any vector field Y in M, then X = 0. Chapter 5. Dynamical Systems 58 a symplectic structure defined by the syfnplectic form to on M . Now if H : M —> IR is a given Cr function on ( M , cu), the vector field XH determined by the condition oo{XH, Y) = <lll • Y (5.2) is called the Hamiltonian vector field with Hamiltonian function H. (M, co, XH) is called a Hamiltonian system. XH is a C r _ 1 vector field, and its existence is guaranteed by the requirement that oo be nondegenerate. The flows of Hamiltonian systems have some very special properties, as we will see in the following theorems and examples. As seen in chapter 2, general relativity can treated as a Hamiltonian system. Thus, all of the results I will present here for Hamiltonian systems may be applied to the dynamics of the Bianchi models. " . A theorem of Liouville states that the flow of a Hamiltonian system will always preserve the volume of the phase space. Here is a precise statement of the theorem. Theorem 5.2 (Liouvi l le ) Let (J\4.,CO,XH) be a Hamiltonian system, and Ft be the flow oj XH • Then for each t, F*to = to; that is, Ft is symplectic on its domain. Thus, Ft also preserves the phase space volume flu- • . • This theorem can be combined with another theorem, due to Poincare, to yield a remark-able result. If the flow of a vector field on a compact space preserves the phase space volume, then any point in the phase space will eventually be carried back to a point arbitrarily close to the original point (although not necessarily back to the point itself). For example, consider an idealized model of a fixed number of ideal gas molecules in a box. AU the molecules are initially confined to one half of the box (side A) by a partition, as illustrated i n figure 5.1. The partition is then removed and the gas molecules are free to move into side B of the box. This is a Hamiltonian system, and Liouville's theorem tells us that the flow of the gas molecules preserves the phase space volume. Thus, if Chapter 5. Dynamical Systems 59 Side B Figure 5.1: The molecules contained in the box are initially confined to side A by a partition. In my simplified model, the molecules may move in any one of eight directions as indicated by the arrows. we wait long enough, then all of the gas molecules will end up back in side A. It may take a very long time for this to occur (possibly longer than the expected lifetime of our Universe), but Poincare's theorem asserts that it will happen. The exact statement of the theorem is as follows. Theorem 5.3 (Poincare Recurrence) Let M be a compact manifold and X a smooth vector field on M with volume-preserving flow Ft. For each open set U in M and T > 0, there is an S > T such that U D Fs(U) ^ (f>. • • Another less general version of this theorem is often more useful in practice 7 . I will state it as a corollary to the theorem stated above. Coro l l a ry 5.3.1 Let g be a volume-preserving continuous one-to-one mapping which maps a bounded region D of Euclidean space onto itself: gD — D. Then, in any neigh-bourhood U of any point of D there is a point x £ U which returns to U. In other words, gnx £ U for some n > 0. 7This is the version given by Arnold in [2] Chapter 5. Dynamical Systems 60 I have written a computer program to simulate the system described above, in an attempt to verify the prediction of this theorem. Since a computer can evolve a simplified model of a system of molecules much faster than the actual system would evolve in nature, I expect to see some positive results before our Universe ceases to exist! In order to reduce complications in the program and to minimize the execution time, I have made a number of simplifying assumptions in my model. As a result, my simplified model does not describe a Hamiltonian system8. However, the flow (being straight lines) is still volume-preserving. Thus, the above theorem still applies and its predictions should hold. This problem belongs to the class of "billiard" problems which consist of a two-dimensional planar domain in which a point particle moves with constant speed along straight line orbits between specular9 bounces. A brief discussion of such systems is given in section 7.5 of [24]; For certain shapes (e.g. a rectangle or a circle) the flow is completely integrable while for other shapes (e.g. a square with an inscribed circle missing or a stadium) it is not. In my simplified model, I assume that the molecules are minimally interacting point particles and that they can pass through each other unperturbed. Two molecules may occupy the same location simultaneously. The molecules may only move in one of eight directions, as indicated in figure 5.1. Collisions of molecules with the box are completely elastic and specular. Including interactions between molecules greatly increases the com-plexity of the dynamics of this system, but this is an unavoidable complication. The dynamics without any such interactions is trivial, as the orbit of each of the molecules forms a closed loop otherwise. I have specified a simple interaction between molecules, which is described in the program listing in the appendix. This is enough to significantly 8The difficulties here stem from the collisions of the molecules with the walls. A satisfactory theory of the motion of a particle in a general potential with singularities (including collisions with walls and with other particles) for Hamiltonian systems is an unsolved problem. Some results related, to this problem are discussed in Section 2.6 of [1] 9angle of incidence equals angle of reflection Chapter 5. Dynamical Systems 61 complicate the dynamics of the system. I set the initial velocities of the N molecules using a pseudo-random number generator, and then let the system evolve in a rectangular box represented by an nx2n grid. The speeds are set such that molecules moving horizontally or vertically move an integral number of grid lengths per timestep and diagonally moving molecules move a multiple of \/2 grid points per timestep. This ensures that there are no round-off errors resulting from the numerical computations.. Now while the above theorem suggests that my program should yield positive results (i.e. all molecules back in side A) for any choice of initial conditions, it may not necessarily arrive at these positive results before I complete this thesis. Thus, it is a matter of trial and error to find initial conditions for which the program will yield such results in a relatively short period of time. I have chosen the box size n to be 25 (i.e. the box is a 25x50 grid) and the maximum speed of the molecules is set to be 100 gridpoints per timestep (or 100\/2 gridpoints per timestep for molecules moving diagonally). I then ran the simulation 1000 times, with different initial conditions each time. The results are shown on the histograms in figures 5.2 - 5.7. A feature common to all of these figures is that there is a big spread in the return times. The maximum time limit of one million iterations was often exceeded, with increasing frequency as the number of molecules was increased. In such cases, there is no way to tell from the simulation alone whether all the molecules will ever end up back oil the same side. But the Poincare recurrence theorem assures us that this will happen. Results for some runs of this simulation averaged over varying numbers of runs for 10 molecules and n=25 are shown in table 5.1. Different initial data was used for each run, and if the return time for any given run exceeded 100000 iterations, then that run was discarded. The number of discarded runs is also given in the table. The average return times in this table are all in the range of 800-900 iterations (except for the first entry). However, these return times would have been greater if it had not been necessary to Chapter 5. Dynamical Systems 62 Number of Runs Average Return T i m e Number of Discarded Runs 10 1612.2 0 100 879.1 0 150 878.3 1 500 816.4 1 1000 862.4 2 5000 848.1 15 10000 863.1 41 50000 872.6 154 100000 870.0 351 Table 5.1: Average return time for 10 molecules over a varying number of runs with a 25x50 box and maximum speed for each molecule of 100 gridpoints per timestep. Different initial data was used for each run, and if the return time for any given run exceeded 100000 iterations, then that run was discarded. discard some of the runs. In fact, minimum estimates for the return times if no runs had been discarded (i.e. assuming the return time for all discarded runs is 100000 iterations) are in the range of 1000-1500 iterations. Table 5.2 shows the average return time for various numbers of molecules averaged over 1000 runs. However, a run was discarded if the return time exceeded one million iterations. Also, return times of 1 were ignored because this indicated that the molecules never left side A of the box and thus "returned" after only 1 interation. Note that while the maximum return time was never exceeded for the simulation with 5 molecules, it was exceeded for most of the runs of the simulation with 25 or 30 molecules. The return time increases rapidly with the number of molecules. This is suggestive of a random walk. For a random walk in one dimension, we would expect a return time of 2^ (of course, this "assumes no walls and no interactions between molecules).1 0 The simulation being discussed is not really a random walk. Molecules bounce off walls, interact with other 1 0See [26] for a discussion of random walks in one dimension. Chapter 5. Dynamical Systems 63 molecules, and are restricted to travel in one of eight discrete directions, even though the box is two dimensional. However, we may still expect the return time T to vary with the number of molecules according to the relation • r = XN ' ' (5.3) where X is some real, constant. The results in table 5.2 appear to be consistent with this hypothesis, with the value of X being around 2. The last three entries in this table have somewhat lower values for X, but a large fraction of the runs were also discarded in these three cases because the maximum return time of one million iterations was exceeded. Thus the average return time shown is lower than it should be in these cases, and so is the value of X. In table 5.3 I have listed the minimum estimated return times rmin based on the data in table 5.2, calculated by assuming that all of the discarded runs have return times of one million iterations. Values of X calculated using these return times r m , n are also shown in this table. They are still consistent with the above hypothesis and a value of X of around 2. The last three entries in table 5.3 are not very reliable because a large fraction of the runs in these three cases were discarded. / I will now present one more theorem which is very useful and important in the world of dynamical systems. It relates spatial averages to time averages, and finds frequent application in the field of statistical mechanics. For example, the classical principle of equipartition of energy may be derived using this theorem. 1 1 . Theorem 5.4 (Mean Ergodic) Let H be a Hilbert space and Ut : H>-> H a strongly continuous one-parameter unitary group (i.e. Ut is unitary for each t, is a flow on H ; and for each x G H, t H-> Utx is continuous). Let the closed subspace Ho be defined by H 0 = {x G H | Utx = xit G R} (5.4) 1 1See [1], section 3.7 for a discussion of this point. Chapter 5. Dynamical Systems 64 Number of Number of runs discarded Number of runs with X molecules N (maximum time exceeded) return time of 1 5 0 111 . 33.2 2.01 10 3 8 883.9 1.97 15 43 0 25241.0 1.97 20 424 • 0 305293.3 1.88 25 944 0 487652.7 1.69 30 998 0 588995.0 1.56 Table 5.2: The average return time rav over 1000 runs. If the maximum return time of one million iterations was exceeded for a run, or if a run had a return time of one iteration, then that run was discarded. A value of X was calculated according to the relation T„V = XN. Number of Train X molecules N 5 33.2 2.01 10 3881.2 2.29 15 67155.6 2.10 20 599848.9 1.94 25 971308.6 1.74 30 999178.0 1.58 Table 5.3: The minimum estimated return times rmin based on the data in table 5.2, calculated by assuming that all of the discarded runs have return times of one million iterations. X was calculated using r T O , n according to the relation r m i n = XN'. Chapter 5. Dynamical Systems 65 Distribution of Return Times for 5 Molecules over 1000 Trials 0 50 100 150 200 250 300 350 400 Return Time Figure 5.2: The frequency distribution of return times over 1000 runs of the simulation for a box with 5 molecules. Distribution of Return Times for 10 Molecules over 1000 Trials 200 • , , , , , , , | , • , | , , , | , • , | • • , 150 : j j I ; -100 fr-10000 12000 Return Time Figure 5.3: The frequency distribution of return times over 1000 runs of the simulation for a box with 10 molecules. Data is not shown for 7 of these runs for which the maximum return time of 1000000 iterations was exceeded. Chapter 5. Dynamical Systems 66 Distribution of Return Times for 15 Molecules over 1000 Trials 150 , i , i i | , , , , | i i , , | , , , , | , , , , | , , , , | , i i i | i i i , j -I • 50000 100000 150000 200000 250000 300000 350000 400000 Return Time Figure 5.4: The frequency distribution of return times over 1000 runs of the simulation for a box with 15 molecules. Data is not shown for 43 of these runs for which the maximum return time of 1000000 iterations was exceeded. Distribution of Return Times for 20 Molecules over 1000 Trials 15 \-0 200000 400000 600000 800000 1000000 Return Time Figure 5.5: The frequency distribution of return times over 1000 runs of the simulation for a box with 20 molecules. Data is not shown for 424 of these runs for which the maximum return time of 1000000 iterations was exceeded. Chapter 5. Dynamical Systems 67 Dis t r ibu t ion of R e t u r n T imes for 25 Molecu les over 1000 T r i a l s 1 1 1 1 1 1 " 1. i i , s i , m iS I, III s 1 n I in I i . I.I I.I •: ii.ii.In. r 0 200000 400000 600000 800000 1000000 R e t u r n T ime Figure 5.6: The frequency distribution of return times over 1000 runs of the simulation for a box with 25 molecules. Data is not shown for 944 of these runs for which the maximum return time of 1000000 iterations was exceeded. e a o O Dis t r ibu t ion of R e t u r n T imes for 30 Molecu les over 1000 T r i a l s 20 I—i—i—i—i—|— i— i—[— 15 10 0 L L -1—i—I—i—i—i—r-_x_ _i_ 200000 300000 400000 500000 600000 700000 800000 900000 1000000 R e t u r n T ime Figure 5.7: The frequency distribution of return times over 1000 runs of the simulation for a box with 30 molecules. Data is not shown for 998 of these runs for which the maximum return time of 1000000 iterations was exceeded. Chapter 5. Dynamical Systems 68 and let-P be the orthogonal projection onto Ho. Then, for any x £ H , • l im - f.Us(x)ds = P(x) . ... (5.5) t - » ± o o t JO • . The limit in this result is called the time average of x, and is denoted x. Thus, x = P(x). This theorem is stated in terms of flow on a Hilbert space. In order to apply it to classical systems, we will need another piece of information to relate the flow Ut on H to the flow Ft on the measure space (M, / / ) . If .Ft is volume-preserving, then it induces a linear one-parameter group of isometries on H = L2(M, u)12 by: Ut(f)=f-F_t. (5.6) From this association follows a useful corollary to the above theorem: Coro l l a ry 5.4.1 Assume that if tn —> t, Ftn(x) —> Ft(x) a.e. for x' £ M and that fi(M) < oo. 1 3 Let Ft be ergodic. Then for f £ L2(M), the limit being in the mean. This result is in a form which can be applied to the dynamics of a classical system. 5.2 A n I l lustrat ive E x a m p l e Now that I have introduced you to many of the basic concepts used in the analysis of dynamical systems, let's explore a concrete, albeit simple example to see how we can actually make use of these results. We will see that even in the case of a simple problem, they "can prove to be invaluable. 12L2(M;';u) denotes the linear, square-integrable.functions on M. 1 3Under these two hypotheses, Ut is a strongly continuous unitary one-parameter group (see [1], section 3.7). Chapter 5. Dynamical Systems 69 The problem which I will examine is the following 1 4. Consider the sequence of num-bers: . Sn = l,2, 4, 8, 1,3, 6 , 1 , 2 , 5 , 1 , 2 , . . . (5.8) .This sequence consists of the first digits of the numbers in the sequence 2 n , n > 0 . (5.9) Now we may ask the following question: What is the relative frequency of occurrence of each number (1 through 9) in this se-quence? One way to approach this problem is to compute numerically a large number of ele-ments of this sequence, and then hope this is a good enough approximation to the whole sequence for the purpose of answering the question posed above. Let's say I decided to follow this approach and I sat down with my calculator and started calculating the elements of this sequence. If I had enough patience to carry out my calculations up to the 200th element, I would end up with the results shown in table 5.4. However, as I will demonstrate, this approach doesn't necessarily work very-well, and the results can be misleading. But first I will show you a very elegant way to solve this problem analytically. The trick is to relate this problem to a physical system which we know is a Hamiltonian dynamical system. We can then make use of the powerful theorems which I described in the last section. The physical system to which I will relate this problem is a single free point particle travelling on a circle S1. The Hamiltonian for this system is ff = £ (5.10) 14 This is a problem posed in [2] Chapter 5. Dynamical Systems 70 Dig i t Occurrences Frequency of Occurrence 1 .60 0.3000000000000000 2 36 0.1800000000000000 3 24 0.1200000000000000 . 4 20 ' 0.1000000000000000 5 16 0.0800000000000000 6 13 0.0650000000000000 7 . 11 0.0550000000000000 8 11 0.0550000000000000 9 9 0.0450000000000000 Table 5.4: The number and frequency of occurrences of each of the digits in the first 200 elements of the sequence Sn, calculated using the computer program listed in the appendix running on a SUN SparclO computer (j) represents the momentum of the particle and m its mass) and the phase space is simply V = S 1 x R . (5.11) The canonical symplectic structure on this space is to = dx A dp (5.12) where x and p are the phase space variables. Since x is a variable on S1, it is periodic. I will choose this period to be 1; thus, any two values of x which differ by an integral amount are equal. The resulting Hamiltonian vector field XH, given by dH = UJ-XH , (5.13) is P_ m (5.14) The flow induced on S1 by this vector field Ft : S1 i-> S1 (5.15) Chapter 5. Dynamical Systems 71 is given by Ft = —t + x. (5.16) m Let us now choose a particular element of this flow by specifying a set of initial conditions. I decide (very judiciously) to choose the following initial condition. - = l o g 1 0 2 . - (5.17) m Notice that now Ft(x) = x + tlogw2 . (5.18) Recall that that the flow is a one-parameter group of diffeomorphisms. I can choose one element of this group by specifying a value for t, as follows. F1{x) = g(x) = x + log102; • (5.19) Upon iterating this mapping n times, we obtain gn(x) = x + ralog102 . (5.20) Now for each n £ Z+, there exists a unique combination of integers jn and kn with 1 < kn ^ 9 such that knWn < 2n < (kn + l)W" .. (5.21) Note that kn will be the first digit of 2 n , which is exactly the [n + l ) ^ n element of the sequence Sn. Keep in mind that the configuration space is the circle S1, which we may represent as the unit line interval with the two endpoints identified with each other, as shown in figure 5.8. Thus two numbers which differ by an integer are equivalent on S1 (i.e. a = b iff a = b mod 1). If we take the logarithm of ( 5.21 ) and drop all additive integers, we obtain log 1 0 K < n log 1 0 2 < l og 1 0 ( f c„ + 1) . (5.22) Chapter 5. Dynamical Systems 72 0 1 Figure 5.8: Circle S1 represented as the unit interval with endpoints identified We may cover the unit interval shown in the figure above with the sets {Ik}9k=1, where Ik = [log10k,\ogw(k + 1)) (5.23) and define an equivalence relation on S1 by: .r ~ k e. Ik • (5.24) We can see from ( 5.22 ) that for some kn, n log 1 0 2 € hn , (5.25) and thus rc log 1 0 2 ~ (5.26) according to the equivalence relation defined above. But Sn = kn as I pointed out earlier, and from ( 5.20 ), . <7"(0) = rcloglo2. (5.27) So, gn(0)e'Ikn^Sn = kn. . (5.28) This last expression provides a concrete link between the Hamiltonian flow of the me-chanical system and the sequence Sn. We are now is a position to apply some of the theorems presented in the previous section to this problem. It is easy to verify that the hypotheses in the statement of the Poincare Recurrence theorem are satisfied: Sl is a compact manifold, X = log 1 0 2 is a Chapter 5. Dynamical Systems 73 smooth vector field on Sl with flow Ft(x) = x + i l o g 1 0 2 , and since the Jacobian (or pull-back) of the flow is unity: f = i , . < * . » ) it is symplectic (or volume-preserving)15. Thus, in any neighbourhood U of any point y £ S1, there exists a point x E U which returns to U. In other words, for any e > 0 there exists some positive integer n such that: | gn(x) - x | < e (5.30) for almost all x 6 S1 1 6 . In particular, for a; = 0, | gn(0) - 0 |< e ' (5.31) for any e > 0 and some integer n. It follows from this that gn(0) will approach any point in S1 arbitrarily closely after some number of iterations. This is clear since gn translates a point on S1 by nlog 2, and this is incommesurable with the circumference of S1 which I have chosen to be unity. So for some n > 0, gn(Q)eh (5.32) for every k G {1,...,9}. Therefore, every number 1 through 9 appears in the sequence As an example, we can ask whether the number 7 appears in the sequence. If x G (log 1 0 7 , log 1 0 8) , (5.33) gn(0) is arbitrarily close to x for some n. Thus, gn(0) G h 1 5 In fact, since Ft is the flow of a Hamiltonian system, Liouville's theorem guarantees us that it is volume-preserving. 16 i.e. for all x £ Sl except possibly for a set of measure zero. Chapter 5. Dynamical Systems 74 <7n(o)~7 ; . Sn = 7 . (5.34) So 7 does indeed appear in the sequence.17 Having shown that each number 1 through 9 appears in the sequence, our next task is to determine with what frequency each of these numbers appears in the sequence. To accomplish this, we will make use of corollary 5.4.1 of the Mean Ergodic Theorem. In order to apply this corollary, we must first verify that the flow is ergodic. It is sufficient to show that the set . 147= {gn(0) | n >• 0} . (5.35) is dense in S1. Then it follows that cf> and S1 are the only invariant sets of the flow, and so by definintion 5.2 it is ergodic. W is dense in S1 if for any two points x and y in W (x < y), there exists a point w £ W such that x < w < y. But from the discussion above, we know that gn(0) will approach any point in S1 arbitrarily closely. In particular, if Bz is an open neighbourhood of the point z £ S1 contained in {p | x < p < y}, then gn(0) £ Bz for some n. So W is dense in S1. Thus it follows from corollary 5.4.1 that: &htf(>iw) = imLf*' (5 36) where f i x ) = Xk{x) = 1 if x £ Ik (5.37) 0 if a; $Ik The left hand side of this expression represents the average "time" <7n(0) spends in the interval Ik, while the right hand side represents the average amount of "space" each 1 7 You may have noticed from ( 5.8 ) that 7 does not appear in the first few elements of the sequence. In fact, 7 does not appear until the 47th element, of Sn. Chapter 5. Dynamical Systems 75 interval Ik occupies in S1. It follows that }™1-n£xk(A0))=P(k) (5.38) i=0 where = (5.39) represents the probability of the the number k appearing in the sequence. Since S1 has measure equal to one, P(k) = u(Ik) = I l °g 1 0 ( fc + l ) - l o g 1 0 f c | = l o g 1 0 ^ ± ^ , fc = l , . . . , 9 . (5.40) Thus, the time the trajectory spends in each interval Ik (or in any set in S1) is propor-tional to the measure of the set. The relative frequency of occurrence of each number in the sequence is given by: P(n) = P(n-l) ) n { , n = 2, 9 . (5.41) l o g i o f e ) . A smaller number will , therefore, appear more often than a larger number: P(n) < P(n - 1) , n = 2, ..., 9 . (5.42) The frequency of occurrence of each digit, calculated using ( 5.40 ), is given in table 5.5. Having obtained exact expressions for the frequencies of occurrence of each digit, we may now return to the question of what kinds of problems we may have encountered had we decided to analyze this problem numerically from the beginning. The first thing to decide is what algorithm to use. I initially tried the most straightforward method: Chapter 5. Dynamical Systems 76 Dig i t Frequency of Occurrence 1 . 0.3010299956639812 2 0.1760912590556812 3 0.1249387366082999 4 0.0969100130080564 5 0.0791812460476248 6 0.0669467896306132 . 7 0.0579919469776867 8 0.0511525224473813 9 0.0457574905606751 Table 5.5: The frequency of occurrences of each of the digits in the sequence Sn, calculated using equation ( 5.40 ) just evaluate 2 n and keep the first digit. Unfortunately, the maximum size of a double precision floating point variable is exceeded on the SUN Sparc 10 workstation that was used for the computations when n — 1024. Thus it is not possible for the computer to compute more than 1024 elements of Sn using this direct approach. Instead, I took ad-vantage of equation ( 5.22 ). The nth element can be found by computing n log 1 0 2 mod 1 and determining in which interval [log 1 0 kn, \og10(kn + 1)) it lies. This algorithm has the advantage of being computationally more efficient than the direct approach, but more importantly it can also evaluate many more elements of the sequence. The C program S E Q U E N C E . C , which implements this algorithm, is listed in the appendix. It was com-piled on a SUN Sparc 10 with gcc, a product of the Open Software Foundation. This program was used to compute the elements of Sn up to the one billionth element. Now in computing the probability of occurrence of each digit, one must decide how many of the sequence elements to use'. One may expect a priori that the more elements one uses, the more accurate will be the result. The number of occurrences of each digit 1 through 9 in the first n elements of Sn (n ranging from one to one billion) were computed using S E Q U E N C E . C , and from this the frequency of occurrence of each digit Chapter 5. Dynamical Systems 77 was calculated. Graphs 5.9 - 5.17 show the difference between the actual probability of occurrence of each digit (calculated from ( 5.40 )') and this frequency of occurrence. Graphs 5.18 - 5.26 display similar quantities, but using elements 1001 through 1001001 of the sequence (i.e. the first one thousand elements of Sn were excluded from the calculations). One feature common to all the graphs is that while the deviations are very large for smaller values of n and get much smaller as more elements of Sn are computed, the deviations continue to oscillate with increasing n. This means that computing a greater number of elements of Sn doesn't necessarily result in a more accurate result for the frequency of occurrence of each digit in Sn. In fact, the result may in some cases be considerably worse when compared to the actual value. Thus, the accuracy of the numerical calculations depends sensitively on the number of sequence elements used in the calculations. So two important questions that one is left with is how to know how many elements of the sequence one needs to calculate in order to get a good approximation, and how good is this approximation? The answer in general is that we simply don't know, unless we have some other information about the system. In other words, we are working in the dark if we rely on numerical computations alone. We need to have some source of qualitative information about the system to shed some light on the results.-Comparing the first set of graphs 5.9 - 5.17 to the second set of graphs 5.18 - 5.26, we see that even if the first 1000 elements of the sequence are left out of the calculations, the main qualitative features of the graphs remain unchanged. In other words, leaving out the first 1000 elements from the calculations does not alter the results. In this fairly simple example of a dynamical system we were able to solve, completely, the problem (ie determine the probability of occurrence of.each digit in Sn) without doing any numerical computations. In general, for more complex systems, we will not be so Chapter 5. Dynamical Systems 78 fortunate. But the qualitative information will at least provide us with a good means to gauge how reliable our numerical calculations are. Numerical results may also depend sensitively upon the exact numerical code used, upon the compiler, and upon the processor. Eor example, the precision of floating point numbers vary between different compilers and different processers. This will result in varying degrees of round-off errors in the calculations. The results one is likely to obtain using a calculator (if one were to calculate the first 200 elements of Sn with a calculator!) are given in table 5.4, and results obtained by using the Mean Ergodic Theorem are given in table 5.5. The deviations between the two sets of results are significant, and we may conclude that the accuracy of the numerical results in this case is not very good. We are able to make this assertion because we already have what we know to be the actual probability of occurrence of each digit for comparison. Were this not the case, we would have no guide to tell us how good the numerical results are. The above discussion demonstrates the advantage of approaching the problem using dynamical systems techniques over numerical approximations. This is a simple example of a question which cannot be answered exactly using numerical techniques, and for which a qualitative analysis using the power of dynamical systems techniques is very useful. 5.3 E rgod ic i ty and Chaos In this section, I will give a formal definition of chaos and discuss its application to general relativity. The definitions which I will present in this section are as stated by Wiggins in [33]. I will begin by defining some relevant terms. Def ini t ion 5.4 A closed invariant set A is said to be topologically transitive if, for any two open sets U,V C A, 3 r £ 1R B ft{U) C\V ^ <f>, where ft is a flow. Chapter 5. Dynamical Systems 79 Digit 1 .o CO _Q O CO o < CD O c CD O O O o c CD 13 O" CD i _ UL 0.5 0.3 0.1 -0.1 -0.3 -0.5 -0.7 -0.9 -1.1 -1.3 •1.5 \1 1*-' Y i \ ' I1 \ I' \ 11 -oeeeosBD—e*-e.oeaaeb — 0 eoeeasB i - ^ / 1 I <t i i i • i \ i \ i -e-ooeeeas- -e-oooemo- -e-o I ,o- -o o-oooo«x> - - o o OOOOOGO o ooooooeo -1 <>- o actual scale B- - - -a 1 o1 scale 102 scale » •* 103 scale 3 104 scale 105 scale « » 108 scale uJ I 10" 10 10 10 10 10 10 10 10 Number of sequence elements counted 10" Figure 5.9: The difference between the frequency of occurrence of the digit 1 in the sequence Sn calculated using the program S E Q U E N C E . C and the actual probability of occurrence given by equation ( 5.40 ). This was calculated for the first n elements of Sn, with n ranging from 0 to one billion. The diamonds represent the data on the actual scale, while the other symbols represent the same data on smaller scales (e.g. data shown as stars, on the 102 scale, should be read as the axis value times 10 - 2 ) . Chapter 5. Dynamical Systems 80 Digit 2 -Q CO JD O * CL "co o < CD O c CD i _ o o o w— o >^  o c CD c r CD 0.5 0.3 0.1 -0.1 -0.3 -0.5 -0.7 -0.9 -1.1 -1.3 -1.5 -1.7 -1.9 -2.1 -2.3 -2.5 TTTTJ-0 <$> — e 9 4 \; >-: . 4 / xd I I _ l i i' i n n i 10 10 10 10" 10' actual scale * . - * 102 scale 0 1> 103 scale < - a 104 scale V V 105 scale A . - A 106 scale Q— — — —Q 107 scale <3 0 108 scale , , . 105 106 107 II1 , llll lit I II. i lU *  1 I-i i -I i 10° 10s Number of sequence elements counted Figure 5.10: The difference between the frequency of occurrence of the digit 2 in the sequence Sn calculated using the program S E Q U E N C E . C and the actual probability of occurrence given by equation ( 5.40 ).' This was calculated for the first n elements of Sn, with n ranging from 0 to one billion. The diamonds represent the data on the actual scale, while the other symbols represent the same data on smaller scales (e.g. data shown as stars, on the 102 scale, should be read as the axis value times 10 - 2 ) . Chapter 5. Dynamical Systems 81 . Q C O _ Q O C O ~5 ••—» o < CD O CD o o o o c. CD 13 cr CD 1.0 0.8 0.6 0.4 0.2 0.0 -0.2 -0.4 -0.6 -0.8 Digit 3 -1.0 "i 1 1— '1 '> TTTTj 1 1 I I I 1 li I hi ' j 1 i» I'M1 1 1 -4 '11' ' 1 $ • !! M | 7 11 1 1 X l l \ 1 1? 1 1 * ,\ 1 1 f\ All; i ' 1 1. V<4 6 V ,° o-o-eecnea>—o-oeeoeeo—o-oeeerao—© ° ° actual scale «----* 102 scale 0 - — > V V & ft 10 scale 104 scale 105 scale 10 scale 107 scale 10 scale 10" 10 10 10 10 10 10 10 10 Number of sequence elements counted 10" Figure 5.11: The difference between the frequency of occurrence of the digit 3 in the sequence Sn calculated using the program S E Q U E N C E . C and the actual probability of occurrence given by equation ( 5.40 ). This was calculated for the first n elements of Sn, with n ranging from 0 to one billion. The diamonds represent the data on the actual scale, while the other symbols represent the same data on smaller scales (e.g. data shown as stars, on the 102 scale, should be read as the axis value times 10~2). Chapter 5. Dynamical Systems 82 Digit 4 C O _ Q O C O ZJ o < CD O c CD ZJ o o o M — o o CD Z5 cr CD 1.2 1.0 0.8 0.6 0.4 0.2 0.0 -0.2 -0.4 -0.6 -0.8 -1.0 -1.2 \ f t i i V V A T A B - — — -a actual scale 10 2 scale 10 3 scale 10 scale 105 scale 10 scale 10 7 scale 10 8 scale I 'I I I I a * | 10" 10 10 10 10 10 10 10 10 Number of sequence elements counted 10" Figure 5.12: The difference between the frequency of occurrence of the digit 4 in the sequence Sn calculated using the program S E Q U E N C E . C and the actual probability of occurrence given by equation ( 5.40 ). This was calculated for the first n elements of Sn, with n ranging from 0 to one billion. The diamonds represent the data on the actual scale, while the other symbols represent the same data on smaller scales (e.g. data shown as stars, on the 102 scale, should be read as the axis value times 10~2). Chapter 5. . Dynamical Systems 83 Digit 5 CO X> O CO 0 < 1 CD O c-CD 13 o o o o >* o c CD ZJ CT CD i _ LL 1.0 0.8 0 :6 0.4 0.2 0.0 -0.2 -0.4 -0.6 -0.8 -1.0 -1.2 -1.4 -1.6 -1.8 df—O-O-0OOG8D -2.0 ™i—7 t M 1 i 1 \ 1 \ i i 1 > i 7 ? »• _o-&oooc«P" -p-e^Geoas^- -o-^eoodsb- W/e-eoda^B— e d-ecra®— i ' A H V I ,1 / I A J ^ - ©- oeecsoes—o-o-eeooBso—o-o-eeooeep O — -- -0 ac tua l s c a l e * — » 102 s c a l e 9- — —> 103 s c a l e — 104 s c a l e V - V 105 s c a l e A —A 10e s c a l e B— — —O 107 s c a l e 0- — —O 108 s c a l e I h in 'ml A 6 9-/' . i 'Li i i i «. i U 'i 10" 10 10 10 10 10 10 10 10 Number of sequence elements counted 10" Figure 5.13: The difference between the frequency of occurrence of the digit 5 in the. sequence Sn calculated using the program S E Q U E N C E . C and the actual probability of occurrence.given by equation ( 5.40 ). This was calculated for the first n elements of Sny with n ranging from 0 to one billion. The diamonds represent the data on the actual scale, while the other symbols represent the same data on smaller scales (e.g. data shown as stars, on the 102 scale, should be read as the axis value times 10~2). Chapter 5. Dynamical Systems 84 Digit 6 CO _Q O CO rs ••—» o < CD O c CD i Z3 O O o o c CD Z5 cr CD 1.0 0.8 0.6 0.4 0.2 0.0 -0.2 -0.4 -0.6 -0.8 •1.0 <p—o-oeeo eeaxeo—o •© 1/ 1/ * © o actual scale * * 102 scale t> > 103 scale < < 104 scale V v 105 scale & & 106 scale D- — a 107 scale O- — o 108 scale I —e e. I i * 7*V 11 II II I  9 II II II Dm / I i |\ A I HM \ f i i i i V di(b .fa i * i' u 11 i i ' 10" 10 10 10 10 10 .10 10 10 Number of sequence elements counted 10" Figure 5.14: The difference between the frequency of occurrence of the digit 6 in the sequence Sn calculated using the program S E Q U E N C E . C and the actual probability of occurrence given by equation ( 5.40 ). This was calculated for the first n elements of Sn, • with n ranging from 0 to one billion. The diamonds represent the data on the actual scale, while the other symbols represent the same data on smaller scales (e.g. data shown as stars, on the 102 scale, should be read as the axis value times l O - 2 ) . Chapter 5. Dynamical Systems 85 Digit 7 _Q CO -Q O CO 13 -*—* O < i CD o c CD i _ ZJ o o o o c CD =5 CT CD LL 1.0 0.8 0.6 0.4 0.2 0.0 -0.2 T 1 I I I I I l | 1 1 I I I I o o actual scale » * 102 scale >••---• 103 scale < < 104 scale v v 10s scale A A 106 scale s Q 107 scale « o 10s scale I I I I III T~ -0-0©60QBffi 0-0-G * f A A A 4 : 0^*0880—o-oeeoosiD—ooeeaeoo—o<3ooecoDD—e oooaB») 10° 101 10 2 103 i o 4 105 10 6 107 108 10 9 Number of sequence elements counted Figure 5.15: The difference between the frequency of occurrence of the digit 7 in the sequence Sn calculated using the program S E Q U E N C E . C and the actual probability of occurrence given by equation ( 5.40 ). This was calculated for the first n elements of Sn, with n ranging from 0 to one billion. The diamonds represent the data on the actual scale, while the other symbols represent the same data on smaller scales (e.g. data shown as stars, on the 102 scale, should be read as the axis value times 10~2). Chapter 5. Dynamical Systems 86 CO .O 0 i _ C L "CO =3 O < 1 CD O c CD i _ 13 O O o o c CD Z3 C r CD 0.2 0.0 -0.2 -0.4 -0.6 -0.8 I1 t 1.4 \ e o actual scale * * 102 scale 0 t* 103 scale V v 104 scale ¥ v 105 scale A A 106 scale B— — H 107 scale o- — o 108 scale Digit 8 \ Y ? n t l y i i 7 I Il I A A A /M II I I I 1 A 10" 10 10' 10 10 10 s 10 10 10a Number of sequence elements counted k i ' t \ i v 10* Figure 5.16: The difference between the frequency of occurrence of the digit 8 in the sequence Sn calculated using the program S E Q U E N C E . C and the actual probability of occurrence given by equation ( 5.40 ). This was calculated for the first n elements of Sn, with n ranging from 0 to one billion. The diamonds represent the data on the actual scale, while the other symbols represent the same data on smaller scales (e.g." data shown as stars, on the 102 scale, should be read as the axis value times 10 - 2 ) . Chapter 5. Dynamical Systems 87 Digit 9 CO .Q 0 L _ C L "co ZJ -t—» O < 1 CD O c CD ZJ o o o o £Z CD =5 c r CD LL 1.3 1.1 0.9 0.7 0.5 0.3 0.1 -0.1 -0.3 -0.5 - - -o actual scale *- •- * 102 scale t>- - *• "IO3 scale - •4- - < 104 scale V - - v 105 scale A— - A "IO6 scale B— - a 107 scale <p—O-O-©60QB8D O-O-'II '•I/ u I 1 - V \ TT y . |7V ^«fe—p\o/eeco88D—o-o-eeeoes)—o'-oeeeaso—© e I i i / \ i KX>eea»9- -e-oooeeeep Vii 10" 10 10 10 10 10 10 10 10 Number of sequence elements counted 10-Figure 5.17: The difference between the frequency of occurrence of the digit 9 in the sequence Sn calculated using the program S E Q U E N C E . C and the actual probability of occurrence given by equation ( 5.40 ). This was calculated for the first n elements of Sn, with n ranging from 0 to one billion. The diamonds represent the data on the actual scale, while the other symbols represent the same data on smaller scales (e.g. data shown as stars, on the 102 scale, should be read as the axis value times 10 - 2 ) . Chapter 5. Dynamical Systems 88 -O co _Q o CO ZS o < CD O CD rs o o o o c CD Z5 C T CD 0.0055 0.0045 0.0035 0.0025 0.0015 0.0005 -0,0005 -0.0015 -0.0025 -0.0035 -0.0045 -0.0055 -0.0065 -0.0075 -0.0085 -0.0095 -0.0105 Digit 1 10' \ / fe ? G •® Ac tua l s c a l e -Q- — -° 101 s c a l e • *- - - * 102 s c a l e A -A 103 s c a l e -0 - o 106 s c a l e _ - o- —e -o—e e o«c|g g— / U /I (I 11 I I I I I I '' ' i •opG| e e-esoe— I. I '• I. I I I I I t 1 • e -e-o-eeeffl i 10 J 10* 10s Number of sequence elements counted 10c Figure 5.18: The difference between the frequency of occurrence of the digit 1 in the sequence Sn calculated using the program S E Q U E N C E . C and the actual probability of occurrence given by equation ( 5.40 ). This was calculated for one million elements of Sn starting with element 1000. The diamonds represent the data on the actual scale, while the other symbols represent the same data on smaller scales (e.g. data shown as stars; on the 102 scale, should be read as the axis value times 10 - 2 ) . Chapter 5. Dynamical Systems 89 Digit 2 C O J O O CL " C O ZS -*—' 0 < 1 C D O c C D i _ ZS O o o >^  o c: C D ZS cr C D 0.0055 0.0045 0.0035 0.0025 0.0015 0.0005 -0.0005 -0.0015 -0.0025 -0.0035 -0.0045 -0.0055 -0.0065 -0.0075 -0.0085 -0.0095 -0.0105 \ \ ' 10 •O-0 \ / t . r rf V o Actual scale -— —a 101 scale (> 102 scale -*- * 103 scale a- A 104 scale --e- - e—o- &-© o o e © G- —e - o - s e eao$ I' \ I" K i \ i \ —I I I I L_ 10 10 10 Number of sequence elements counted 10c Figure 5.19: The difference between the frequency of occurrence of the digit 2 in the sequence Sn calculated using the program S E Q U E N C E . C and the actual probability of occurrence given by equation ( 5.40 ). This was calculated for one million elements of Sn starting with element 1000. The diamonds represent the data on the actual scale, while the other symbols represent the same data on smaller scales (e.g. data shown as diamonds, on the 102 scale, should be read as the axis value times 10~2). Chapter 5. Dynamical Systems 90 Digit 3 -O CO -Q O CO ZJ > O < i CD O d CD i_ =5 o o o o c CD zs D" CD 0.0090 0.0080 0.0070 0.0060 0.0050 0.0040 0.0030 0.0020 0.0010 0.0000 -0.0010 -0.0020 -0.0030 -0.0040 -0.0050 -0.0060 ; \ \ / \ / 'h V [ f \ I • [-/ •/ e - - -<> Actual scale B - - - -a 10 scale o- - - H> 1 o2 scale - - - A 103 scale \ / I I I 1 -o-e-o-oeeo-\ k 1 10' 10' 104 10b Number of sequence elements counted A V * ' r A -Q- -e-o-o-eoeo 10c Figure 5.20: The difference between the frequency of occurrence of the digit 3 in the sequence Sn calculated using the program S E Q U E N C E . C and the actual probability of occurrence given by equation.( 5.40 ). This was calculated for one million elements of Sn starting with element 1000. The diamonds represent the data on the actual scale, while the other symbols represent the same data on smaller scales (e.g. data shown as diamonds, on the IO2 scale, should be read as the axis value times IO - 2 ) . Chapter 5. Dynamical Systems 91 .Q CO -Q O CO ZS O < I CD O ' tz CD ZS o o o o c • CD Z5 cr CD LL 0.0025 0.0015 0.0005 -0.0005 -0.0015 -0.0025 -0.0035 -0.0045 -0.0055 -0.0065 -0.0075 -0.0085 -0.0095 * a -0.0105 10 0 p. Digit 4 • « o— e- *-o-o0o* T 1 i \ o a- - e o- o-eooe"? ® \ — ° - 0 - 0 e e— -o—e -e -e-eood) V * * #— +- -»|e—?|c » * * _i i ' i t ' e — -o Actual scale a ° 101 scale * * 102 scale A a 103 scale « o 1011 scale 10 10 10 Number of sequence elements counted 10c Figure 5.21: The difference between the frequency of occurrence of the digit 4 in the sequence Sn calculated using the program S E Q U E N C E . C and the actual probability of occurrence given by equation ( 5.40 ). This was calculated for one million elements of Sn starting with element 1000. The diamonds represent the data on the actual scale, while the other symbols represent the same data on smaller scales (e.g. data shown as stars, on the 102 scale, should be read as the axis value times 10~2). Chapter 5. Dynamical Systems 92 Digit 5 CC -O O CO =3 O < • CD O c CD k_ Z3 O o o o >. o c CD cr CD 0.0095 0.0085 0.0075 0.0065 0.0055 0.0045 0.0035 0.0025 0.0015 0.0005 -0.0005 -0.0015 --0.0025 --0.0035 -0.0045 -0.0055 -0.0065 h -0.0075 -0.0085 --0.0095 --0.0105 10 1 j _ i - \ e-o-o-eo A -Actual scale 101 scale 102 scale 10 3 scale 104 scale —o- —e -f*-e- Qoeoo- o-1 / * v \ / \ I \ / /. A / A / N A v. © -e- e e-©oo* A L A / \ / 10° 10' 10& Number of sequence elements counted 10c Figure 5.22: The difference between the frequency of occurrence of the digit 5 in the sequence Sn calculated using the program S E Q U E N C E . C and the actual probability of occurrence given by equation (5.40 ). This was calculated for one million elements of Sn starting with element 1000. The diamonds represent the data on the actual scale, while the other symbols represent the same data on smaller scales (e.g. data shown as stars, on the 102 scale, should be read as the axis value times 10 - 2 ) . Chapter 5. Dynamical Systems 93 Digit 6 CO - Q O CO zs 0 < 1 CD O CD k_ Z5 o o o > » o c CD ZJ cr CD 0.0095 0.0085 0.0075 0.0065 0.0055 0.0045 0.0035 0.0025 0.0015 0.0005 -0.0005 -0.0015 -0.0025 -0.0035 'i 1 \ 1 / ' ' / ' I /I 1/ \ i -1 t \ I \ I \ I \ I o— -- -e Actual scale • 3 - — -a 101 scale * -* 102 scale _ A 103 scale -0 -o 104 scale -* r 1 1 1 ' i r ' 1 i / \ Q Kl X I* 1 ©-©oco e- f e—o- p-e|eoda> i i i I 10' 10° 1 0 \ 10b Number of sequence elements counted 10° Figure 5.23: The difference between the frequency of occurrence of the digit .6 in the sequence Sn calculated using the program S E Q U E N C E . C and the actual probability of occurrence given by equation ( 5.40 ). This was calculated for one million elements of Sn starting with element 1000. The diamonds represent the data on the actual scale, while the other symbols represent the same data on smaller scales (e.g. data shown as stars, on the 102 scale, should be read as the axis value times 10 - 2 ) . Chapter 5. Dynamical Systems f 94 Digit 7 _Q CO - Q O co Z J o < CD O c CD k_ ZJ o o o >^  o cr. CD ZJ cr CD L _ LL 0.0035 0.0025 0.0015 0.0005 -0.0005 -0.0015 -0.0025 -0.0035 -0.0045 -0.0055 10' A- A A-AA \ O r - -e J o / 1 \ -4 e .e -esoe e - -e- e o-ooeea) \ Actual scale 101 scale 10 2 scale 10 3 scale 10° 10' 10° Number of sequence elements counted 10c Figure 5.24: The difference between the frequency of occurrence of the digit 7 in the sequence Sn calculated using the program S E Q U E N C E . C and the actual probability of occurrence given by equation ( 5.40 ). This was calculated for one million elements of Sn starting with element 1000. The diamonds represent the data on the actual, scale, while the other symbols represent the same data on smaller scales (e.g. data shown as stars, on the 102 scale, should be read as the axis value times 10 - 2 ) . Chapter 5. Dynamical Systems 95 Digit 8 CO - Q O CO ZS -+—' o < CD O £= CD i Z3 O O o >< o c CD ZJ cr CD i _ LL 0.0045 0.0035 0.0025 0.0015 0.0005 -0.0005 -0.0015 -0.0025 -0.0035 -0.0045 -0.0055 -0.0065 -0.0075 -0.0085 -0.0095 10' - - ° Actual scale - - ° 10 scale - -* 10 2 scale - 10 3 scale / \ V / 1 I \ 1 \ < /I n / i • e -e-oooeci ' • A / A- - - A 10 10 10 Number of sequence elements counted 10c Figure 5.'25: The difference between the frequency of occurrence of the digit 8 in the sequence Sn calculated using the program S E Q U E N C E . C and the actual probability of occurrence given by equation ( 5.40 ). This was calculated for one million elements of Sn starting with element 1000. The diamonds represent the data on the actual scale, while the other symbols represent the same data on smaller scales (e.g. data shown as stars, on the 102 scale, should be read as the axis value times 10~2). Chapter 5. Dynamical Systems 96 Digit 9 CO - Q O CO ZJ 0 < 1 CD O c CD i— r3 o o o o c CD zs cr CD LL 0.0055 0.0045 0.0035 0.0025 0.0015 0.0005 -0.0005 -0.0015 9 9 10' \/ Jj\ (B q \ G -o Ac tua l s c a l e B - — — n 101 s c a l e - * 102 s c a l e A — A 103 s c a l e /? /1 7 I ri I / I / I / I / I / I ' I <L i I I i i? i> '\ i l b I ' I I ' I I l i \ ^ > - \ 11 i ^ 5,1 I p. / I . I ' w \ • 4 * I A X-^ 0 —e - © e-o e o e e -H I Il I * I I A e — o - e—o-o-eeoe> 10 10 10 Number of sequence elements counted 10c Figure 5.26: The difference between the frequency of occurrence of the digit 9 in the sequence Sn calculated using the program S E Q U E N C E . G and the actual probability of occurrence given by equation ( 5.40 ). This was calculated for one million elements of Sn starting with element 1000. The diamonds represent the data on the actual scale, while the other symbols represent the same data on smaller scales (e.g. data shown as stars, on the 102 scale, should be read as the axis value times 10 - 2 ) . Chapter 5. Dynamical Systems 97 Requiring a flow to be topologically transitive is a stronger condition than ergodicity, which I defined earlier. However, it is equivalent to saying that the flow is ergodic on the whole phase space (i.e. an invariant set must cover the entire phase space), and it is in this stronger sense that I have been and will be using the term ergodic. Def ini t ion 5.5 Let A C IRn be an invariant set under the flow ft(x). The flow ft(x) is said to have sensitive dependence on initial conditions on A if there exists e > 0 such that, for any x £ A and any neighbourhood U of x, there exists y £ U and t > 0 such that | ft(x) - ft(y) |> e. The two terms which I have just defined characterize chaos in dynamical systems, which is defined as follows. • Def in i t ion 5.6 (Chaos) A is said to be chaotic if 1. ft(x) has sensitive dependence on initial conditions on A.. 2. ft{%) is topologically transitive on A. Thus, ergodicity (on the whole phase space) is a necessary condition for chaos. Also, since it is a topological feature, it is well-defined when applied to general relativity. Sen-sitivity to initial conditions, on the other hand, is more difficult to determine. Lyapunov exponents are often used to characterize this (see [33], for example, for a general dis-cussion). They are a measure of the rate at which neighbouring trajectories diverge or converge. Positive Lyapunov exponents indicate sensitivity to initial conditions. One way to define them is the following [25]. Consider two initial conditions po and p'0 which are separated by a small distance lo on the phase space of a dynamical system. Also let pT and p'T, the final states obtained by evolving po and p'0 for a time r , be separated by a distance lT. The principal Lyapunov exponent A of the system is defined as: A = lim lim - l o g ^ . (5.43) Chapter 5. Dynamical Systems 98. If neighbouring trajectories separate exponentially or faster, A will be strictly positive. Otherwise, it will be zero or negative. There is a problem with this definition when it is applied to general relativity because the. latter is invariant under time reparametrizations but clearly the definition of A above is not, due to the ^ factor appearing in it. A trajectory which has a positive Lyapunov exponent with respect to some time parameter may have A = 0 with respect to a different time parameter. We can therefore see that calculating Lyapunov exponents is not a well-defined way of characterizing chaos in a general relativistic system because of sensitivity to reparametrizations of the time parameter. Nevertheless, a number of people have done calculations of this sort, without reaching any conclusive results. For example, such calculations have been done for the Bianchi IX cosmologies in [5], [6], and [15]. So at present we do not know of any good way of characterizing sensitivity to initial conditions in general relativity. Ergodicity on the other hand is well-defined and can be used as a very good indicator of chaos. Although it does not by itself imply chaos, determining whether or not a system is ergodic is an important step towards determining whether the system is chaotic. Also, a non-ergodic system is not chaotic. Chapter 6 Bianch i I : A Case Study A l l of the background material necessary to analyze the phase space dynamics of the Bianchi cosmologies has now been developed. In this chapter I will illustrate how this analysis may be carried out by examining the flow on the phase space for Bianchi I. This is the simplest of the homogeneous cosmological models. Since the structure constants for the Bianchi I algebra are all zero, the resulting spacetime is just flat Minkowski space with a time-dependent scaling of the spatial components of the metric. The resulting dynamics are also very simple, but provide a clear illustration of the kind of analysis which may be applied to the Hamiltonian flows by using concepts presented in chapter 5.' The objective of this analysis is to determine whether or not the phase space flow is ergodic. Ergodicity is an important indicator of chaos ih dynamical systems. It is also a topological effect, and thus it is vitally important to know exactly what the topology of the phase space in question is. To emphasize and illustrate the point I will present in the following section a very well known example: the motion of a free particle on a 2-torus. I will then proceed to construct the phase space and identify the topology of this space for the Bianchi cosmologies. Following this, I will study the phase space dynamics for Bianchi I and determine whether or not it is chaotic. 99 Chapter 6. Bianchi L: A Case Study 100 6.1 The Torus E x a m p l e A standard and very well-known example in dynamical systems is a free particle on a 2-tor.us T 2 . This is a perfect example to demonstrate how a flow can be non-ergodic on one space, but ergodic on. another space obtained from the first one by making certain identifications of points, and thus changing the topology of the space. Consider first the motion of a free particle on 1R2. The Hamiltonian is H = d ± A (6.i) 2m v ; (where as usual p\ and p2 are the components of the momentum and m is the mass). The canonical symplectic on the phase space R 2 x R 2 is LO = dq1 A dpi + dq2 A dp2 . (6.2) The corresponding flow projected onto the configuration space R 2 is FM\q2)=U + -t,q2 + -t) • (6-3) V • m m J This is just a family of straight lines on R 2 , and it is clearly not ergodic. Now let's make some identifications of points on the configuration space, as illustrated in figure 6.1. In effect, I have made q1 Z/a-periodic and q2 Z^-periodic (where La and Lb are the lengths of sides a and b, respectively, in figure 6.1. In other words, there is an equivalence relation ^ such that (g \g 2 ) w (q1 + rLa,q2 + sLb) V integers r,s > 0 . (6.4) The phase space now becomes T*(T2) = T 2 x R 2 . The expressions ( 6.2 ) and ( 6.3 ) for the symplectic structure and configuration space flow remain the same, except that q1 and q2 are now coordinates on the 2-torus T2 rather than on R 2 . Chapter 6. Bianchi I: A Case Study 101 R 2 T > 1 i a b b a V . , J Figure 6.1: A 2-torus is constructed from R 2 by identifying the two sides labeled a and also the two sides labeled b. The coordinate q1 is made L a-periodic and q2 is made Lfe-periodic ( L a and Lb are the lengths of sides a and b, respectively) by defining an equaivalence relation on R 2 . Chapter 6. Bianchi I: A Case Study 102 Let's suppose that Ft is periodic on T2 with period r\ Then, Fr(q\q2) = (q\q2), . (6.5) and so P i T -r and q1 + —r = q1 + rLa ' ' (6.6) m q2 + = q2 + sLb . (6.7) m (6.8) Thus, P2 La s So we are led to the condition that a trajectory if periodic if and only if is a rational number. If this ratio is irrational, then the trajectory is not periodic and will approach any point in T2 arbitrarily closely (this follows from the Poincare recurrence theorem). Thus, the flow in such cases is ergodic. This example illustrates how by modifying the topology of the phase space, the flow on the configuration space has become ergodic. This is an important point, and in the following sections I will show how the phase space flows for the Bianchi cosmologies similarly become considerably more complicated because of some identifications which must be made on K 2 in order to construct the proper configuration space. 6.2 Canonica l Variables The Einstein-Hilbert action for general relativity is I = J y/^jRd^x . (6.9) in units where G = - J ^ . In terms of the coordinates I chose to use in chapter 4, recall that the spacetime line element takes the form ds2- = -dt2 + e~a (e") ( r V , • (6.10) Chapter 6. Bianchi I: A Case Study 103 and so — = e-n+0(i)e-n+/J(2)e-n+j9(3) _ e-an (6.11) The scalar curvature for Bianchi I in these coordinates is given by: Rj = - 6 0 + 12ft2 + 6/?2 + 6/3! (6.12) Now the choice to work in a synchronous gauge with the lapse N(t) set equal to 1 was quite arbirary, and was made for convenience only. I can now decide to pick a different lapse function; this will correspond to making the following change in my choice of coordinates: dt2 -> N2{t)dt2 . (6.13) Thus, d_ It di2 1 d N(t)dt 1 d ( 1 d N(t)dt \N(t)dt) N(t)y/=g = N(t)e~3Q (6.14) (6.15) (6-16) and the Bianchi I Lagrangian density in terms of the lapse N becomes Li = -6e- 3"4f4n)+6 e- 3"f^ 2 + l/9 2+l^ dt \N = -6 d_ 2t 1 N (6.17) + ^ e - 3 " (20 2 + fa + >!) . (6.U The first term in ( 6.17 ) is a total time derivative, which can be discarded because it contributes only a boundary term to ( 6.9 ). It is apparent from this expression that a convenient choice for the lapse is the following: N = 6e -3d (6.19) Chapter 6. Bianchi I: A Case Study 104 The resulting expression for the Bianchi I Lagrangian is then ,Lj = ( - f i 2 + ft + Pi) . . (6.20) The non-vanishing momentum variables conjugate to the variables appearing in the La-grangian (which I called 7r,j in chapter 2) are P + = = 2.V+ . (6.21) P- = djf = W- (6-22) dL Pn = ~( = -2n. (6.23) Oil Note that the scalar curvatures R for the other Bianchi models are of the form Ri + {terms independent of p+, /?_, $7} . (6.24) So it follows that the momentum variables for all of the Bianchi models are given by ( 6.20 ), ( 6.21 ), and ( 6.22 ) with the choice ( 6.18 ) for the lapse. It is important to note that this new choice of lapse will modify the components of the Einstein tensor that were calculated in chapter 4. As one may expect, only the terms which contain derivatives with respect to the parameter t are affected; the spatial components remain unchanged. The modified components of the Einstein tensor are as follows. ' Goo = ^(tf-Pl-fc)-1-^ (6.25) • Gaa' = ^ e 6 a [ 2 n + faa)-d(-n2 + Pl + p2_)] + 3Gaa , (6.26) where — R and 3Gaa are given in tables 4.2 and 4.4, respectively. There is no change in the Gab (a ^ b) and G 0 ; components, given by equations ( 4.90 ) and ( 4.86 ) respectively. Chapter 6. Bianchi I: A Case Study 105 6.3 Phase Space The unconstained configuration space Co is simply R 3 . This is clear because the domain of each of the configuration space variables /5+, /?_, and fl is R . Therefore, the uncon-strained phase space is VQ = T*CQ = R 3 x R 3 . Now the constraints ( 2.39 ) and (. 2.40 ) must be satisfied. If we let co = {(p+,0-,n,p+,p-,m)eT*M\Goo = o} (6.27) and d = {(/?+,./?_, Cl, P + , p_, Pn) £ T*M I Goi = 0} (6.28) denote the set of solutions of the constraints, then the constraint set for the vacuum Einstein system is C = C0nCiC T*M . (6.29) The resulting phase space is exactly this constraint set with all gauge freedom removed: V = Cj ~ . (6.30) The gauge has been completely fixed, but there may yet remain some residual gauge invariance in hij under isometries, as I discussed in chapter 3. Any such isometries must preserve the form of the metric. If I write the line element corresponding to A,j in some basis {xi} as follows cfe2 = X>?(t)d*?, • (6-31) i=l then such an isometry will transform the coordinates as follows: Xi -> oti = fi(xi,x2,x3) .(6.32) where j1(x\, x2, x3) are functions of the spatial variables. Now, 3 dxi = ^2 fi,jdxj • (6.33) Chapter 6. Bianchi! : A Case Study 106 where • ( 6 ' 3 4 ) If the line element in terms of the new coordinates {xi} is given by 3 ds2 = J2AKt)dx2 (6.35) then it follows that ds2 = J2A2(t)(J2h,dxA . • (6.36) i=i \j=i ) In order for the transformation to be an isometry, this expression must have the same form as ( 6.30 ). Thus, all of the cross terms must vanish. Since all the functions are arbitrary, each term in each of the coefficients of the dxidxj (i ^ j) must vanish separately. So for each i (i = 1,2,3), at least two out of the three / , j must vanish. Therefore each /,• can depend on at most one coordinate a;,. Further, the dx2 terms must have non-vanishing coefficients. This implies that for each j, at least one fij ^ 0 (i = 1,2, 3). If for example fi,i — 0, then the coefficient of dx\ vanishes. Thus, each / , must be a function of a different coordinate Xi. Since the coefficients of the dx2 terms in ( 6.35 ) are functions of t only, the /;,/s must all be independent of the coordinates {x{\. Thus each function fi must be a linear function of one of the coordinates xt. The line element is then given by ds2 = J2A2(t)f2dx2 (6.37) which is of the same form as ( 6.30 ), as required. The result is that: •x1-+ fi(xi) = xa (6.38) X2 -+ f2(x2) = xb (6.39) xz -> f3{x3) = xc (6.40) where a, 6, c 6 {1, 2,3} and a ^ b ^  c: Thus the isometries represented by the functions' fi correspond to a permutation of the basis elements Chapter 6. Bianchi I: A Case Study 107 Now that we have identified the remaining gauge freedom in the 3-geometry, we can remove this freedom by identifying all the three-geometries which differ only by a permutation of the basis vectors. In order to carry out this identification, it is convenient to first rewrite the spacetime metric ( 4.6 ) in the following form: ds2 = -N2(t)dt2 + a2{t)(ax)2 + b2{t)(a2)2 + .c2(t){a3)2 , (6.41) where . . . v a(t) = e - 2 0 + 2 W ) (6.42) b{t) = e - 2 Q + 2 p W (6.43) c{t) = e - 2 f i + 2 ^ 3 ) . • . (6.44) It is clear that permuting the elements of the basis is equivalent to permuting the func-tions a(t), b(t), and c(t). Inverting the above expressions yields: Q, = -Iln(afcc) (6.45). "-' = i7!lnG) • <«7> We now examine how f i , /3+, and /?_ are affected when (a, b, c) are acted upon by the group of permutations on three symbols, which I denote TT3. I also denote the elements of 7r3 as follows: — { T a l c , Tbcaj Tcabi Tacbi Tcba-i Tbac\ i (6.48) where the subscript of each group element indicates its effect on {abc}. It is clear that 7r3 is a group: Tabc is the identity element, T b c a and Tcab are inverses and all other elements are their own inverse, and it is closed under the associative group operation of composition. Chapter 6. Bianchi I: A Case Study 108 From ( 6.44 ), it is clear that is unaffected by the action of IT3 since it depends only upon the product abc. Therefore, we can restrict our attention to /3+ and /3_ only. The group action of TT3 on yields the following: Tabc({3+, /?_) Tbca ( A , / ? - ) Tcab T o c 6 (/?+,/?_) T c 6 a (/?+,/?_) Tbac(p+,/3-) A l l of the elements of 7r 3 map (/3+, to points which are equidistant from the origin in the /3+-/3_ plane, with the distance being given by r EE v//5+2 + /?_2 . (6.55) This suggests that the action of each of the elements of 7r3 involves a rotation. We can see this more clearly if we re-express (3+ and•/?_. in terms of two parameters r (given by ( 6.54 ) ) and 9: f3+ = r cos 9 (6.56) = rsinfl . . (6.57) Then, the action of the elements of TT3 may be expressed as follows: Tabc{r,9) = (r,0) (6.58) Tbca(r,e) = , ( ^ + y ) (6-59) ( A H £ - ) (6.49) (-•^ . (/?+ + V3/?_) , | (V3/?+ - /?_)) (6.50) ( - i (/?+ - VzfiJ) , -1 (V3^ +- + /?_)) (6.51) ( - | (/9+ - V3/3_) , | (V3/9+ + /?_)) (6.52) • ( - | (^ + + V3^l), | - ( -V3/9 + + z?:)) (6.53) (6.54) Chapter 6. . Bianchi I: A Case Study 109 (r, 2 7tf3) (r,n) (r,4it/3) r, 5TC/3) Figure 6.2: Identification of points in the ( / ? + , /3_) plane as a result of the group action oi 7 T 3 . The points marked by the same symbol are identified with each other. Tcab (r, 9) Tacb (r, 9) Tcba (r, 0) Tbac (r, 61) = . r,0 + 47T = r,-9 + 2vr • 4TT = (r,-9) . (6.60) (6.61) (6.62) (6.63) We can thus visualize the effect of 7r 3 on (/?+, /?_). Each point lying on a circle centred about the origin with radius given by ( 6.54 ) gets identified with five other points also lying on the same circle, except for {(r, 0), (r, ^ ) , (r, 4^)} and {(r, 7r), (r, | ) , (r, ^ ) } . Only the point (0,0) is mapped back solely to itself under the group action of 7r3, and thus it is the only fixed point. This identification is illustrated in figure 6.3 for some selected points. The points marked by the same symbol are identified with each other. The fundamental domain of TT3 is D = {{r,d)\r>0,0<e<^} . In other words, ^ ( D ) = S1 . (6.64) (6.65) Chapter 6. Bianchi I: A Case Study 110 So, one complete turn around the circle is now completed after just | radians starting from 0 = 0! One way to think of the topology of the resulting space is to imagine the region D to be a wedge-shaped billiard table that extends out to infinity on the open end of the wedge. Any trajectory reaching one of the walls = 0 , B+ > 0 (i.e. 6 = 0 , r > 0) (6.66) and B_ = V3/3+ , P± > 0 (i.e. 6 = | , r > 0) . (6.67) will bounce off of that wall. The resulting space, which I will denote V, is a smooth manifold with boundary. This is clear since the mapping / . : V H 1H+ (where 1H+ is the upper half-plane in R 2) defined by /(r ,0) = (r,30) (6.68) is a diffeomorphism from V onto 1H +. The bounces off of the walls resulting from the identification of points in the (/?+,/?_) plane described above amount to reflections of trajectories about an axis given by ( 6.65 ) or ( 6.66 ). The angle of incidence di (defined as the angle between the tangent line to the trajectory at the point where it intersects a wall and the normal to that wall) equals the angle of reflection otr (defined similarly). This is clear from figure 6.3. Now in the case of Bianchi I, the Hamiltonian constraint is H = Goo = 3e 3 f i (ft2 - p\ - / ? ! ) = 0 , (6.69) and there is no momentum constraint. The momenta are thus constrained to be in the domain: P = R 3 n{ (p + , p_ ,p f i ) | p n 2 - P + 2 -p_ 2 = 0} . (6.70) Chapter 6. Bianchi I: A Case Study 111 Figure 6.3: The bouncing of a trajectory off one of the walls amounts to a reflection off that wall. It is clear from this figure that a,; = aj. . Chapter 6. Bianchi I: A Case Study 112 This is a two dimensional surface W in Euclidean 3-space. It consists of two cones with tips touching at the origin, as illustrated in figure 6.4. The phase space V is therefore V X W , with the natural (canonical) symplectic structure given by LO = dd+ A dp+ + 'dB- A dp_ • (6.71) where (/?+,/?_) € V and (p+,p_) € W. 6.4 Dynamics The equations of motion for Bianchi I are those for the motion of a free point particle: p+=p_=pQ = 0. (6.72) Integration of these equations yields: 8+(t) = P+t + 3°+ (6.73) P-(t) •= P-t + 3°_ (6.74) Q(t) = pQt + ft0 . (6.75) The flow on V is Ft{P+,P^,to,P+,P-,Pti) = (p+t + P+, p-t + 8_, pat + ft, p+,p-, PQ) . (6.76) This flow represents the full dynamics on the constrained phase space Pi = C\ x Pj. Alternatively, we may incorporate the constraints directly into the flow as follows: Fn(8+j^P+,p_) = ( • ; p 2 + =(ft-fto) + /3+, . p ; =(n-no) + /3-, v±^+)2 + (p-)2 ± v W + (p-)2 P + P ~ ^ . (6.77) ±^(p+)2 + (Py ±^/(P+)2 + (p_y Chapter 6. Bianchi I: A Case Study 113 Figure 6.4: The two dimensional momentum space W consists of two cones with tips touching at the origin. Chapter 6. Bianchi I: A Case Study 114 where fl now plays the role of the "time" parameter. This constrained flow represents the same dynamics, but on the unconstrained space. 0 in this case is a very natural choice for the time parameter; the (p+,p-,pn) space looks like a three dimensional Minkowski space (i.e. a 2 + 1 dimensional spacetime). That these two procedures (imposing the constraints on the dynamical equations and then studying the dynamics on the full phase space and imposing the constraints on the phase space and then studying the dynamics on this constrained space) yield the same dynamics is a feature common to all non-quantized dynamical systems. If one is examining a quantized system though, these two procedures will in general yield different results (see [28] for an example and discussion of this point). It is quite clear that the flow (6.76 ) on a flat space IR4 is not ergodic. It is just a family of straight lines, each with fixed p+ and p_ (constant momentum). There are two independent constants of motion, which we take to be p+ and p_ (of course, one could use H or pn instead, but only two will be independent). Therefore, according to the standard classical mechanics terminology, the phase space flow is smoothly integrable. The important question though is whether or not this same flow on the actual phase space V x W is ergodic. In other words, is it possible to construct invariant measurable sets for this flow on the phase space other than the empty set and the whole space itself? A trajectory in (/?+,/?_) space is completely determined by specifying initial values for (/?+,/?_) and p_ /p + , since the equation of any line segment is given by P- = — P+ - — P+° + P-° • (6.78) P+ P+ It is clear that all initial choices of.(p_|_,p_) with the same ratio p_ /p + will yield the same trajectory in configuration space. So we may define an equivalence relation ~ : (p+,P-)^.(p'+,p'L)&P-/p+=p'_/p'+ (6.79) Chapter 6. Bianchi I: A Case Study 115 so that any two elements belonging to the same equivalence class will result in the same trajectories in configuration space. In momentum space W these equivalence classes are straight lines on the surface of the cones, passing through the origin. When projected down onto the plane, they are just radial lines. The behaviour of radial, trajectories in configuration space (i.e. those which pass through the origin) is simple to describe. A radial line 0 = 0o = constant (6.80) will be reflected at the origin to another radial line given by 0 = ^ - 6 Q . (6.81) Each radial line wi l l bounce off a wall exactly once (at the origin), and then go out to infinity. An illustration of this is provided in figure 6.5. Non-radial trajectories are a little more complicated to describe. A l l such trajectories will bounce off the positive (3+ axis (wall 2) at least once. Thus we can pick a generic point P on this axis (away from the origin) and see what happens to all trajectories which bounce off this point. This will provide us with a complete characterization of the configuration space dynamics. Consider the trajectories shown in figure 6.6 as they approach the point P from the right of line 1. For convenience, I will discuss separately trajectories approaching P from each of the following three regions: Region I: Area between lines 1 and 2. These trajectories, will have angles of incidence 60 in the range 0° < 90 < 30°. Region II: Area between lines 2 and 3. These trajectories will have angles of incidence 80 in the range 30° < 90 < 60°. Region III: Area between lines 3 and 4. Chapter 6. Bianchi I: A Case Study 116 Chapter 6. Bianchi I: A Case Study 117 Figure 6.6: A sample of trajectories in configuration, space which bounce off the point P at various angles. Chapter 6. Bianchi I: A Case Study 118 Bounce number W a l l Ang le • 1 1 60° <0_! < 90° 2 ' 2 • 0° < 60 < 30° 3 1 30° <0X< 60° Table 6.1: Bounces for trajectories approaching the point P from region I. 6o is the bounce angle at point P. A l l angles are measured from the normal to the wall. These trajectories will have angles of incidence f90 in the range 60° < 60 < 90' A l l trajectories bouncing off point P will bounce a total of three times. Trajectories in region I will approach wall 1 from infinity, undergo a set of three bounces as specified in table 6.1, and then go out to infinity again. Trajectories in regions 2 and 3 will behave similarly, except that they will approach wall 2 from infinity, bounce 3 times as specified in tables 6.2 and 6.3 below, and then go out to infinity again. The only exception is line 2, which will bounce only twice before going back out to infinity. Lines 1 and 3 will bounce twice at the same point. The slope of any trajectory in configuration space corresponds to the ratio as we can see from ( 6.77 ). For any given trajectory, this will be constant between bounces and will change at each bounce. The different values of the slopes of the line segments making up the trajectories approaching P at an angle 0O are given in table 6.4. Thus the number of distinct values of p-/p+ a given trajectory takes on equals the number of bounces it makes plus one, with the following exceptions: • The radial trajectory 8 = |- takes on only one value of • Trajectories which bounce off wall 2 at an angle 0 = 30° take on three different values Chapter 6. Bianchi I: A Case Study 119 Bounce number W a l l Ang le 1 •2 30° < 0O < 60° 2 1 0° < 0! < 30° 3 2 60° < 02 < 90° Table 6.2: Bounces for trajectories approaching the point P from region II. 90 is the bounce angle at point. P. A l l angles are measured from the normal to the wall. Bounce number W a l l Ang le 1 2 60° < 0O < 90° 2 1 0° < 0t < 30° 3 2 . 30° < 62 < 60° Table 6.3: Bounces for trajectories approaching the point P from region III. 0O is the bounce angle at point P. A l l angles are measured from the normal to the wall. Ang le Slope 1 r Slope 2 Slope 3 Slope, 4 0° < 60 < 30° tan(30° + 60) tan(90° - d0) - tan(90° - d0) . t an(3O°-0 o ) 30° < 0O < 90° tan(90° - 0O) - tan(90° - 0O): tan(30° - 0O) - tan(30° - 60) Table 6.4: The different values of the slopes of the line segments making up the trajec-tories approaching P at an angle 0O. Chapter 6. Bianchi I: A Case Study 120 ofp_/p+. • Trajectories which bounce off wall 2 at an angle 6 = 0° or 6 = 60° take on two different values of p_ /p + . Now that we know how trajectories in configuration space behave, we can attempt to determine whether the flow is ergodic. Since none of the trajectories bounces off the walls in configuration space more than three times, and the corresponding momenta remain bounded in a subset of W, it is clear that there are no dense orbits in the phase space. In other words, if we were to compute the phase space orbit for any set of initial conditions (/?+,•/?_,p+,p_), it would not be dense. If any such orbit were dense, it would indicate an ergodic flow. However, as this is not the case for Bianchi I, it is unlikely that the flow is ergodic. We may make certain that this is the case by constructing an invariant measurable set in the phase space. Each trajectory takes on at most four different values of since it bounces no more than three times. As a specific example, let us begin by choosing the region Sc in configuration space represented by the shaded region in figure 6.7. This is the set of all points which are covered by the full evolution in configuration space of trajectories which approach the point P at an angle in the range 70° < 90 < 80°. The full set of values of momentum Sm corresponding to these trajectories (i.e. such that the ratios p~/p+ of the points (p + , p_) in Sm are the slopes of all the line segments comprising all these trajectories) are shown in figure 6.8. Recall that there is an identification of points in momentum space under the equivalence relation ~ . Now to construct an invariant set from this, let us first figure out what the set W = {FQ (Sc x Sm) | fl G 1R} (6.82) is. Consider the evolution of any point (qo,Po) G V x Sm. The set {Fn(?o,Po) I 0 G R} defines a trajectory in V, and it is clear from the preceding discussion that the slope p~/p+ Chapter 6. Bianchi I: A Case Study 121 Figure 6.7: The shaded region represents the set of points Sc in configuration space covered by the evolution of all the trajectories which approach the point P from an angle in the range 70° < 6Q < 80°. Chapter 6. Bianchi I : A Case Study' 122 p Figure 6.8: The shaded region represents the full set of points Sm in momentum space projected onto the (p+,p_) plane corresponding to the slopes of all the line segments comprising the trajectories which approach the point P from an angle in the range 70° < 90 < 80°. The are the slopes of the indicated lines. Chapter 6. Bianchi I: A Case Study 123 of any of the lines in configuration space which make up this trajectory will correspond to momentum values in SM. In other words, W fl W = SM. The set W is thus an invariant set which is not space-filling in the momentum space. According to definition 5.2, this means that the flow FQ is not ergodic on the phase space V. I have not yet completely identified the set W. In particular, I have not yet indicated which set of points in configuration space are contained in W. As it turns out, W^nV = V. This result.is a consequence of the discussion which follows. Consider first any non-radial trajectory, such as the one represented in figure 6.9 by the solid lines labeled line 1 through line 4. It will make one bounce against wall 2 at an angle of 30° or less; this is clear by examining tables 6.1 - 6.3. In the case depicted in figure 6.9, line 1 makes such a bounce. Let mj be the slope of line 1. Now since we started with the product set SC X 5"m,' any invariant set of non-zero measure containing this trajectory must also contain a line with slope rn\ intersecting line 3 and passing through the origin. This line is represented by line 5 in the figure. Further, all lines parallel to line 3 and intersecting line 5 (represented by dashed lines in the figure) must also be in this set. But this family of lines covers all of V. Thus, if the invariant set contains a non-radial trajectory, it must contain all of V. Now let's consider radial trajectories, such as the one represented by lines 1 and 2 in figure 6.10. These trajectories consist of two lines with slopes m-[ and ?n 2. If the invariant set is of non-zero measure and contains this trajectory, it must also contain all lines of slope m 2 intersecting line 1, such as line 3 in figure 6.10. This line intersects wall 2. According to the discussion in the preceding paragraph concerning non-radial trajectories, it follows that this invariant set must contain all of V. The only exception is the trajectory defined by 9 = | , represented by line 4 in figure (6.4). The sets {(/?+,/?_, P +,p_) | 9 = tan" 1 (/?_//?+) = ^ , (p+,p_) G L } (6.83) Chapter 6. Bianchi I : A Case Study 124 F i gure 6.9: If we construct an invariant set starting with the region'Sc x Sm, then since it contains a non-radial trajectory, such as the one represented by the solid lines 1-4, it will contain the whole configuration space. Chapter 6. Bianchi I: A Case Study 125 m, = tanfc/3- 9) m. = tan@) Figure 6.10: If we construct an invariant set starting with a product set such as the region Sc x Sm, then if the set contains any radial trajectory other than 9 = | (line 4), such as the one represented by lines 1 and 2, it will contain the whole configuration space. where LC = ^} (6.84) are invariant sets in V which do not include all of V, but are of measure zero. It is also possible to construct invariant sets which do not contain all of V. These sets ca be constructed by evolving an initial region of phase space for all time. An example is the following. Consider an initial region intersected with P— tan 9i < — < tan 92 P+ ei < p2+ + pi < e2 (6.85) (6.86) Chapter 6. Bianchi I: A Case Study 126 and tan^i < < tan (^ 2 (6.87) P+ intersected with . . ; " * i < d2+ -r ;j2_ < c2 . • . " (6.88) This is a small wedge shaped region S in phase space! Now we can construct an invariant set Y by taking the evolution of this region for all time fi. It is clear that the resulting set Y - Fu(S) ' . . " (6.89) is invariant and that YD V ^ V. Consider evolving the /?+'/?_ wedge for just one value of momentum. The resulting space will be the union of four strips consisting of one between infinity and wall 1, two between wall 1 and wall 2, and one between wall 2 and infinity, as shown in figure 6.10. Now for each different p~/p+ ratio in the set S, we will get a similar set of four strips. The union of all of these strips is the invariant set. ( 6.85 ) bounds the slope of each of the strips. In particular, the strip with maximum slope tan 0\ and the strip with minimum slope tan 02 bound, the slopes of all of the strips. Finally, any point in a given strip evolves to another point in the same strip. So for each p_ /p + , the four corresponding strips is an invariant set. Therefore, the union of all;these invariant sets for all momenta in the set S is an invariant set in the phase space. To summarize then, the flow FQ is not ergodic on V since there do exist non-trivial invariant measurable sets; V x SM is an example. Following from definition 5.6 then, since the phase space flow for the Bianchi I cosmologies is not ergodic, we may conclude that this flow is also not chaotic. Chapter 6. Bianchi I: A Case Study 127 Figure 6.11: The evolution of the B+8_ wedge, represented by the shaded region, for one particular momentum value yields four strips. This is shown for the maximum and minimum values of the ratio These place bounds on the slopes of the strips corresponding to all the other values of the ratio in the invariant set Y. Chapter 7 Bianchi II : A Brief Look I have just presented to you a complete analysis of the phase space dynamics for the Bianchi type I cosmologies. It turned out that the dynamics in that case are very simple, and thus the analysis could be carried out without encountering any great obstacles. This sets the stage for carrying out a' similar analysis for the remaining Bianchi cosmologies. The analysis will in principle follow along similar lines as for the case of Bianchi I. However, the dynamics will not be nearly as simple and this can lead to significant difficulties when attempting to analyze the dynamics. The equations of motion for the Bianchi II cosmologies are: B+(t) = ^_ln{sech 2 [V^(c - t)]} + fat + fa fl{t) = -^ ln j sech 2 [ \ /6a (c - i ) ]} + 7 i i + 7o P-(t) = — Injsech 2 [Voa(c-i)]} + Plt +.p0 (7.1) (7.2) (7.3) with the following four constraints on the seven constants a, fa, jQ, pQ, fa, 7^ px; a > 0 (7.4) (7-5). (7.6) (7.7) Unfortunately it is not possible to examine the phase space flows in any straightfor-ward way. The main problem one encounters is that there is no clear choice of a time A _ I e 4 ( - 7 o+^o+ \ / 3 p 0 ) 6 -a = , 2 ^ + ^+pf) 0 = _ 7 l + fa + v / 3 / 0 l . 128 Chapter 7. Bianchi II : A Brief Look 129 parameter which will allow us to separate out the configuration space variables from the momentum space ones and we did for Bianchi I (by choosing fi as the time parameter). This is because the Hamiltonian constraint contains a non-vanishing potential term (i.e. the expression for Goo contains a non-vanishing spatial part 3_R, as we can see in table 4.2. As a result, it is not possible to eliminate the parameter t from equations ( 7.1 - 7.3 ) explicitely and rewrite these as functions of fi, for example, (i.e. and /3-(fi) ) as we did for Bianchi I. This makes it very difficult to study the dynamics of these equations in the reduced phase space analytically. It also makes it very difficult to identify what the momentum space is, since we cannot separate it out from the configuration space as we did for Bianchi I. Thus, we do not even know what the topology of this space is. The configuration space though will still have the same topology as for Bianchi I. A numerical approach to solving this problem would not be satisfactory because this would introduce round-off errors, defeating the purpose of this study which is to determine exactly the qualitative features of the phase space flows. Without being able to clearly identify the space of physical degrees of freedom, it becomes very difficult to analyze the dynamics and determine what the phase space portrait looks like. Similar problems are encountered with the equations of motion of all the other remaining Bianchi-cosmologies as well. This brief look at the phase space dynamics of the Bianchi II cosmologies gives us an indication of the hurdles which must be overcome in order to analyze the flows for the other remaining Bianchi cosmologies. However, it also points the way towards how such an analysis may be carried out. The bottom line is that while I have set the framework for analyzing completely the qualitative features of the flows on the reduced phase space for all the Bianchi cosmologies, actually carrying through this analysis is an exceedingly complex task, except in the case of Bianchi I. Chapter 8 Discussion and Concluding Remarks Whether or not cosmological models displaying chaotic behaviour occur generically in the space of solutions of Einstein's equations is a very important question, and one which is still wide open. If such behaviour is generic, then the implications are profound. It would mean that severe limitations exist on the predictive power of the theory of general relativity. Determining whether any particular solution is an accurate model for our own Universe would be very difficult because its evolution would in general depend sensitively on the initial conditions, which we have no way of ascertaining to arbitrary precision. Unfortunately, obtaining a completely general solution of Einstein's equations is a task which is well beyond our means. The best we can hope to do is determine whether chaotic behaviour exists within particular classes of solutions, and from this attempt to draw inferences about the general case. The Bianchi cosmologies have been the focus of attention of a number of researchers for many years. This class of solutions of Einstein's equations has attracted attention because exact solutions can be obtained for these cosmologies. However, the dynamics of these solutions on the space of gravitational degrees of freedom is certainly not well-understood. A number of studies ([5], [6], [13], [15] among them) have attempted to determine numerically whether or not the dynamics of the Bianchi cosmologies (mainly of the Bianchi IX "mixmaster" cosmology) exhibit chaotic behaviour. I have attempted to illustrate the limitations of numerical studies by means of examples of dynamical systems, while at the same time to demonstrate-how indispensable analytical studies are. 130 Chapter 8. Discussion and Concluding Remarks 131 Indeed, an examination of the numerous numerical studies of the dynamics of the Bianchi cosmologies in the light of the results which I have presented in this thesis highlights some very substantial deficiencies and raises some doubts about the validity of these numerical studies. I will elaborate on these concerns in the following paragraphs. First though, I will review some of the more salient features of my analysis of the dynamics of Bianchi I. Two conditions are necessary to establish the presence of chaos: ergodicity on the whole phase space and sensitivity to initial conditions. Studies thus far have focused primarily on the latter condition by the computation of Lyapunov exponents. How-ever, as I discussed earlier, it is well known that these do not provide a good test of chaos for general relativistic systems because of its invariance under time reparametriza-tions. I have instead demonstrated how the first condition, ergodicity, can be used in a well-defined way to test for the presence of chaos in the Bianchi cosmologies. I have rigourously demonstrated that the phase space .flow for Bianchi I is not chaotic. I did this by demonstrating the non-ergodicity of the phase space flow, rather than by the more common approach of showing that the flow does not display a sensitive dependence on initial conditions. I have thus shown that ergodicity does provide a good test for the presence of chaos, and that it is well-defined for general relativistic systems. Further, I have demonstrated that the phase space dynamics of Bianchi I exhibit a complicated behaviour, despite statements made by some that the Bianchi I dynamics are simple ([25], for example). The interesting dynamics are a result of the non-trivial topology of the phase space, which to my knowledge is an aspect of the dynamics which has not yet been addressed in the literature (and which likely explains why the Bianchi I dynamics have been described as being very simple). A strong motivation for attempting to analyze the dynamics of the Bianchi models numerically is that in most cases an analytical study such as the one I carried out for Bianchi I does not seem feasible. I indicated in chapter 7 where the problems arise in Chapter 8. Discussion and Concluding Remarks 132 attempting this. However, one is left with the question of how it is possible to demonstrate the presence of chaos in a system by studying numerically an approximation of the true dynamics of the system. After all, whatever happened to sensitivity to small changes in the initial state, a characteristic of chaotic systems? The aforementioned numerical studies deal with analyzing the dynamics on (j3+, /?_) space for the Bianchi cosmologies. They do not consider the complete phase space dy-namics. While the configuration space dynamics may be interesting to study, it is not clear what the physical relevance is. In order to analyze the dynamical evolution of a cosmology, one must consider the evolution of all the parameters representing the full space of gravitational degrees of freedom - i.e. the phase space. The configuration space represents only one half of the degrees of freedom. There are systems for which it is in-teresting to study the configuration space dynamics. The dynamics-of a free particle on a 2-torus discussed in section 6.1 is one such example. However, it is not clear what the physical relevance of this model is. On the other hand, a number of dynamical systems models are known whose phase space dynamics are physically relevant; the Henon-Heiles potential, which is used to describe motion in an axisymmetric potential'[14], and the geodesic flow on a space with constant negative curvature1 are two such examples. Even if one is.able to argue convincingly that there is a sound physical reason for studying solely the configuration space dynamics rather than the full phase space dy-namics, there is still a serious. problem with the studies carried out thus far. I have shown by construction that the topology of the configuration space is non-trivial for all of the Bianchi cosmologies, and that this non-trivial topology has a very significant ef-fect on the dynamics. However, the above-mentioned studies do not account for these :The geodesic flow on a space with constant negative curvature is equivalent to a constant energy flow of a purely kinetic Hamiltonian system.' According to Hadamard's theorem, this flow is ergodic on constant energy submanifolds of the Hamiltonian phase space. See [3] for details. Chapter 8. Discussion and Concluding Remarks 133 topological effects; in fact, they make no mention of them. They assume that the configu-ration space is just R 2 . So even if any of them is successful in demonstrating the presence of chaos on this Reconfiguration space, it does not follow that the chaotic behaviour will still be present on the true configuration space. Again, this point is clearly brought out by the example of the 2-torus in section 6.1. There remains much work to be done in uncovering the chaotic behaviour which many suspect exists in general relativity, be it a generic feature of the theory or an occurrence only in particular cosmological models. I have shown that ergodicity, which is a necessary condition in the characterization of chaos, is by itself a good test for the presence of chaos in general relativity. I have pointed out some serious deficiencies in much of the recent numerical work that has been done in this area, in addition to demonstrating the limitations of numerical analyses of dynamical systems in a more general context. Further, I have shown that it is essential to know the topology of the phase space in order to obtain the correct phase space dynamics, and the topology is one thing which generally cannot be determined numerically. The implication of this work for future studies into chaotic behaviour in general relativistic systems is that the common approach currently taken in studying such systems is not satisfactory, and that more attention must be devoted to the effects'that topology has on'the dynamics. < Bibl iography [1] Abraham, R.; Marsden, J.E.; Foundations of Mechanics, 2nd ed., Addison-Wesley Publ. Company Inc., 1978 (Sixth Printing, 1987). [2] Arnold, V. I . ; Mathematical Methods of Classical Mechanics, Springer- Verlag, New York, 1989. [3] Arnold, V.I . ; Avez, A. ; Ergodic Problems of Classical Mechanics, W. Benjamin, Inc., New York,. 1968. [4] Barrow, J.D.; Physics Reports 85, No. 1, 1-49, 1982. "Chaotic Behaviour in General Relativity" [5] Berger, B . K . ; Class. Quan. Grav. 7, 203-216, 1990. "Numerical Study of Initially Expanding Mixmaster Universes" [6] Burd, A. ; Buric, N . ; Tavakol, R .K . ; Class. Quantum Grav 8, 123, 1991. "Chaos, entropy and cosmology" [7] Bianchi, L. ; Mem. di Mat. Soc. Ital. Sci. 11, 267, 1897. "Sugli Spazzi a Tre Dimen-sion! che Ammettono un Gruppo Continuo di Movimenti" [8] Collins, C.B. ; Hawking, S.; Astrophys. J. 180, 317, 1973. "Why is the Universe Isotropic?" [9] Ellis, G.F.R.; MacCallum, M . A . H . ; Cornmun. math. Phys., 12, 108, 1969. "A Class of Homogeneous Cosmological Models" [10] Fischer, A . E . ; Marsden, J.E.; "The Initial Value Problem and the Dynamical For- • mulation of General Relativity", pp. 138-211 in: Hawking, S & Israel, W. (eds.); General Relativity : An Einstein Centenary Survey, Cambridge University Press, New York, 1979. " [11] Folland, G.F.; Real Analysis, John Wiley & Sons, Inc., New York, 1984. [12] Gleick, J.; Chaos, Viking, New York, 1987. [13] Halpern, P.; Gen. Rel. Grav., 19, No. 1, 73, 1987. "Chaos in the Long-Term Behavior of Some Bianchi-type VIII Models" [14] Henon, M . ; Heiles, C ; Astr. J., 69, 73, 1964. 134 Bibliography 135 [15] Hobill, D; Bernstein, D.; Welge, M . ; Simkins, D.; Class. Quantum Grav., 8, 1155, 1991. "The Mixmaster cosmology as a dynamical system" [16] Hoenselaers, C. & Dietz, W. (eds.); Solutions of Einstein's Equations: Techniques and Results (Proceedings of the International Seminar on Exact Solutions of Ein-stein's Equations Held in Retzhach, Germany, November 14-18, 1983), Springer-Verlag, 1984. [17] Jantzen, R.T.; A Complete Primer oh the Theory Involved in the Study of Spa-tially Homogeneous Solutions of the Einstein Equations, Junior Paper, Princeton University, Fall 1972. [18] Schmutzer, E . (ed.); Exact Solutions of Einstein's Field Equations, Cambridge Uni-versity Press, New York, 1980. [19] Marsden, J.; Ebin, D.; Fischer, A. ; Proceedings of the Canadian Mathematical Congress, 13, 135, 1972. "Diffeomorphism Groups, Hydrodynamics and Relativity" [20] Marsden, J.; Weinstein, A. ; Reports on Mathematical Physics, No. 5, 121, 1974. "Reduction of Symplectic Manifolds with Symmetry" [21] Misner, C.W.; Thorne, K.S. ; Wheeler, J .A.; Gravitation, W . H . Freeman and Com-pany, New York, 1973. [22] Morrow-Jones, J.; Witt , D . M . ; Phys. Rev. D, 48, No. 6, 1993. "Inflationary Initial Data for Generic Spatial Topology" [23] Moser, A. ; Matzner, R.A. ; Ryan, M.P. ; Annals of Physics, 79, 558, 1973. "Numerical Solutions for Symmetric Bianchi Type IX Universes" [24] Ott, E. ; Chaos in Dynamical Systems, Cambridge University Press, 1993. [25] Pullin, J . , "Time and Chaos in General Relativity" in: SILARG VII, Relativity and Gravitation: Classical and Quantum, D'Olivio, J.C.; Nahmad-Achar, E. ; Rosen-baum, M . ; Ryan, M.P.; Urrutia, L .F . ; Zertuche (eds.); World Scientific, 1991. [26] Reif, F.; Fundamentals of Statistical and Thermal Physics, McGraw-Hill Book Com-pany, 1965. [27] Ryan, M.P.; Shepley, C.S.; Homogeneous Relativistic Cosmologies, Princeton Uni-versity Press, Princeton, New Jersey, 1975. [28] Schleich, K . ; Class. Quantum Gravity 7, 1529, 1990. "Is reduced phase space quan-tisation equivalent to Dirac quantisation?" Bibliography 136 [29] Schutz, B. ; Geometrical Methods of Mathematical Physics, Cambridge University Press, Cambridge, U . K . , 1980. [30] Taub, A . H . ; Annals of Mathematics 53, No. 3, 472, May 1951. "Empty Space-Times Admitting a Three Parameter- Group of Motions" [31] Wald, R . M . ; General Relativity, University of Chicago Press, Chicago, 1984. [32] White, M . ; Scott, D.; Silk, J.; Annu. Rev. Astron. Astrpphys., 32, 319-370, 1994. "Anisotropics in the Cosmic Microwave Background" [33] Wiggins, S.; Introduction to Applied Non-Linear Dynamical Systems and Chaos, Springer-Verlag, 1990. A p p e n d i x A L i s t i ng of Computa t iona l P rog ram P O I N R E C . C The following is a listing of the C program used for the simulation discussed in Chapter 4. /* POINREC.C */ #include<stdio.h> #include<math.h>" #include<string.h> #include<ctype.h> #include<stdlib.h> #include<time.h> #include<stddef.h> /* This version contains a modification to poin.c to add simple interactions between the molecules */ /* This version contains a modification to poinrec.c to average the return time over num_runs separate runs. It w i l l save a l l the results to an output f i l e . */ /* This program i s meant to provide a demonstration of the Poincare Recurrence theorem. This program simulates a simple model of molecules moving around i n a closed box. The user chooses the number of molecules and the dimensions of the box (ie the number of grid points used to represent the box). A l l the molecules are i n i t i a l l y confined to the l e f t half of the box (side 1), and they are "seeded" (ie t h e i r i n i t i a l positions, directions, and speeds are set) either using the pseudo-random number generating function rand() or by input from a user-specified data f i l e . The molecules may move i n one of eight directions, as follows: 1 2 3 \ I / 4 - * - 5 / I \ 6 7 8 At each timestep, each molecule moves a distance i) equal to i t s speed i f i t i s moving in direction 2,4,5, or 7. i i ) equal to i t s speed times sqrt(2) i f i t i s moving in dir e c t i o n 1,3,6, or 8. After each time time step, the program checks whether or not a l l the molecules are i n side one. If they are, the program terminates. The user may also 137 Appendix A. Listing of Computational Program POINREC.C 138 specify a time l i m i t , at which the program w i l l terminate whether or not a l l the molecules are back in side 1. Bugs: The molecules evolve in discrete timesteps. Thus, i t only checks i f a l l the molecules are back in side one after a discrete amount of time. If they were a l l back in side one after say 4.5 timesteps, but after 5 timesteps one molecule jumps back to side 2, then the program w i l l not catch t h i s . The chance of t h i s happening can be reduced by ensuring that the dimensions of the box and speeds of the molecules are not much greater than the number of molecules. This i s not a big problem though because i t i s not necessary to catch the f i r s t occurence of a l l the molecules being back in side.one. If none of the molecules leave side 1 after the f i r s t timestep, then the program w i l l end with a l l the molecules."back" in side 1 after one timestep. Interactions: Molecules w i l l interact as follows: 1) Molecules w i l l interact only i f they are at the same point after a whole number of timesteps. 2) Only 2-particle c o l l i s i o n s are allowed. If more than 2 molecules are at the same point, then they w i l l interact i n pairs. 3) The effect of an interaction between 2 molecules w i l l be: a) to interchange t h e i r directions of motion (this i s e f f e c t i v e l y the same as interchanging t h e i r speeds) i f they c o l l i d e head-on (180 deg) or at right angles (90 deg), or b) to increase or decrease the angle of incidence by 90 degrees, .. depending on whether the angle i s acute or obtuse, resp., otherwise. If more than 2 molecules are present at a point, then they w i l l interact c y c l i c a l l y i n pairs (i.e 1 with 2, 2 with 3,- etc.). */ • struct molecule{ int x; int y; int v; int side; int speed; }•; - ' #define MAX_num_runs 150 int getseed(void); int evolve(struct molecule *, i n t ) ; void evolvel(struct molecule *, i n t , int * void evolve2(struct molecule *, i n t , int * void evolve3(struct molecule *, i n t , i n t * void evolve4(struct molecule *, i n t , int * void evolve5(struct molecule *, i n t , int * void evolve6(struct molecule *, i n t , int * void evolve7(struct molecule *, i n t , int * void evolve8(struct molecule *, i n t , int *,, void bounce(struct molecule *, struct molecule * ) ; Appendix A. Listing of Computational Program POINREC.C 139 main(int argc, char *argv[]) { enum boolean {no, yes}; enum boolean allback, stop; enum boolean seed = yes; /* Random seeding of molecules */ enum boolean setmaxtime = no; /* Set a maximum time? */. unsigned int MAXSPEED; register int i=0; int maxtime =0; int timex =0; /* Counts # of times maxtimes i s exceeded over the MAX_num_runs runs */ register i n t timestep =0; int N, n; char ans; int size = sizeof(struct molecule); FILE *filename; char fname[80]; FILE * o u t f i l e ; char outfname[80]; struct molecule *mol; struct molecule *molp; f l o a t avtime = 0; . register int num.runs = 0; int gseed; register int j , k; int countint =0; /* Counts the t o t a l number of interactions that took place in each run */ struct molecule *molpl; enum boolean match = no; i f (!( (argc==3 II argc==4 || argc==5) && ( (N=atoi(argv[1])) > 0) && ( (n=atoi(argv[2])) > 0) )) { printf("Run the program using the following syntax please:\n\n"); printf("poinrec [N] [n] [filename] [maxtime]\n\n"); printf("where:\n\n"); printf("N specifies the number of molecules in the box . W ) ; p r i n t f ( " n specifies the size of the box (which has dimensions 2n x n).\n"); printf("filename i s optional. It specifies the name of the f i l e from " ) ; p r i n t f ("which i n i t i a l \ndata for the '/,d molecules w i l l be read. ", N); p r i n t f ( " I f none i s specified, then \nmolecules w i l l be seeded randomly.\n"); printf("maxtime i s optional. It specifies a maximum number of timesteps"); p r i n t f ( " after\nwhich execution w i l l terminate. If none i s specified, " ) ; printf("program\nwill terminate when a l l molecules are back i n side 1."); printf("\n\n"); printf("\nEnter a value for [N]: " ) ; scanf ("'/,d",&N); printf("\nEnter a value for [n]: " ) ; scanf ("7,d",&n); getcha ' rO ; printf("\nDo you wish i n i t i a l data to be read from a f i l e " ) ; printf("(<y>es, <n>o, <q>uit)?"); ans = getchar(); getchar(); i f ( tolower(ans) == 'y' ) { Appendix A. Listing of Computational Program POINREC.C 140 printf("\nEnter filename:"); scanf 0"/,s80" ,fname); if ( (filename=fopen(fname, "r")) == MULL ) { printf("Error: cannot open fi l e . Try again."); exit(l); > seed=no; else if ( tolower(ans) == 'q' ) exit(l); printf("\nDo you wish to specify a maximum time? (<y>es, <n>o, <q>uit)?"); ans = getcharQ; getcharO ; if ( tolower(ans) == 'y' ) { printf("\nEnter maxtime: "); scanf ('"/,d" ,&maxtime); setmaxtime=yes; > else if ( tolower(ans) == 'q' ) exit(l); > else-C switch(argc) { case 3: break; case 5: maxtime=atoi(argv[4]); setmaxtime=yes; case 4: if ( strlen(strncpy(fname, argv[3],80)) < 81 ) { if ( (filename=fopen(fname, "r")) == NULL ) { printf("Error: cannot open fi l e . Try again.")'; exit(l); • } • } else { printf("Error: filename too long or invalid. "); printf("Please try again."); exit(l); > . seed=no; break; default: printf("Error in input on command line. Please try again."); exit(O); > > MAXSPEED = (unsigned'int)(4.0*n); mol=(struct molecule *)malloc(N*size); printf("\nEnter output filename:"); scanf C"/,s72" .outfname) ; sprintf (outf name, "'/.s'/.d'/.s'/.d" ,outfname,n, M.",N); if ( (outfile=fopen(outfname, "a")) == NULL ) { printf("Error: cannot open fi l e . Try again."); exit(l); Appendix A. Listing of Computational Program POINREC.C 141 > f p r i n t f (outf i l e , "Results for boxsize n = '/,d, maximum speed = 4*n = '/,d, and\n", n, (unsigned int)(4.0*n) ); f p r i n t f (outf i l e , "number of molecules N = '/,d. \n\n" ,N); gseed=getseed(); srand(gseed); /* Seed random number generator using system clock */ f p r i n t f (outf i l e , "Seed = 7,d \n\n", gseed) ; for (num_runs=l; num_runs <= MAX_num_runs; num_runs++) { timestep=0; molp=mol; /* Check i f user wants molecules randomly seeded, or i f i n i t i a l conditions are to be read from an input f i l e */ i f (seed==no) { i=0; rewind(filename); while (i<N) { ' i f ( fscanf (filename ,"'/,d '/,d '/,d */,d '/.d", &mol[i].x, &mol[i].y, &mol[i] .v, &mol[i].side, &mol[i].speed) == EOF) { p r i n t f ("Insufficient data in f i l e for '/,d molecules. Check.", N); printf("data f i l e . \ n " ) ; e x i t ( l ) ; > else { i f ( !( (mol[i].x < (n+n)) kk (mol[i].x >= 0) ) ) { p r i n t f ("Error in x-position entry of '/,dth molecule in data " , i + l ) ; p r i n t f ( " f i l e . \nPlease try again."); e x i t ( l ) ; . } . -i f ( !( (mol[i].y < n) && (mol[i].y >= 0) ). ) { p r i n t f ("Error i n y-position entry of '/.dth molecule in data " , i + l ) ; p r i n t f ( " f i l e . \nPlease t r y again."); exit(1); > i f (. !( (mol[i].v >= 1) kk (mol[i].v <= 8) ) ) { p r i n t f ("Error in value of v for '/,dth molecule i n data " , i + i ) ; p r i n t f ( " f i l e . \nPlease t r y again."); e x i t ( l ) ; } i f ( !( (moi[i].side ==1) || (moi[i].side == 2) ) ) {. p r i n t f ("Error i n value of side for '/,dth molecule in data " , i + l ) ; p r i n t f ( " f i l e . \nPlease try again."); exit(1); > i f ( mol[i].speed < 0 ) { p r i n t f ("Error: value of speed for '/,dth molecule in data must " , i + l ) ; printf("be non-negative.\nPlease t r y again."); e x i t ( l ) ; > i f ( moi[i].speed > MAXSPEED ) { p r i n t f ("Warning: Speed of */,dth molecule in data exceeds maximum", i + i ) ; Appendix A. Listing of Computational Program POINREC.C 142 p r i n t f (11 speed set\nby variable MAXSPEED. " ) ; > > i++; fclose(filename); > /* Otherwise randomly seed the N molecules to specify the i n i t i a l conditions */ else { i=0; molp=&mol[0]; while (i<M) { (*molp).x = rand()'/,n; /* Produce random number i n the range 0 to n-1 */ (*molp).y = rand()'/,n; (*molp).v = ((rand()7,25)'/,8) + l ; /* Produce random number i n the range 1 to 8 */ (*molp).side = 1; /* A l l molecules start out in side 1 */ (*molp).speed = rand()'/,(MAXSPEED+l); /* Molecule's speed can be in the range 0 to MAXSPEED */ molp++; i++; > > /* Next carry out the evolution of the N molecules, checking at each time step whether a l l the molecules are back in side one and whether maxtime has been exceeded */ stop = no; while( (stop == no) && ( (timestep < maxtime) I I (setmaxtime==no) ) ) { ++timestep; allback=yes; /* Carry out interactions between molecules */ countint=0; molp=mol; for (j=0; j<(N-l); j++) { match=no; molpl=molp+l; . k=j + l.; while ( (match==no) && (k<N) ) { i f ( ( (*molp).x == (*molpl).x ) II ( (*molp).y == (*molpl).y ) ) { bounce(molp, molpl); match=yes; countint++; . > else k++; > molp++; > p r i n t f 0"/,d interactions for run #7,d, timestep = '/,d\n", countint, num.runs, timestep); molp=mol; i=0; Appendix A. Listing of Computational Program POINREC.C 143 while(KN) { if (evolve(molp,n) == 2) allback=no; ++molp; ++i; > if (allback==yes) stop=yes; > if (stop==no) { /* Maximum time exceeded */ f printf (outf ile, "> '/,d\n" , timestep); timex+=l; > if (stop==yes) { f printf (outf ile,'"/,d\n" , timestep); avtime+=timestep; } >. fprintf(outfile, "MAX_num_runs = '/,d\n" ,MAX_num_runs); if (timex < MAX_num_runs) { f printf (outf ile, "7.d & '/.f \n", N, ((float)avtime/(float)(MAX_num_runs-timex))); f printf (outf ile, "Maximum time of '/,d exceeded '/.d times . \n" , maxtime, timex); > else f printf (outf ile, '"/,d & average return time unknown \n" , N) ; fclose(outfile); } /* The following function produces a seed for the pseudo-random number generator randQ using the local time */ int getseed(void) { time_t It; char chr[30]; char tim[6] ; lt=time(NULL); strcpy(chr,ctime(&lt)); sprintf (tim,",/,c7,c'/.c,/1c", chr [18] ,chr [17] ,chr[15] ,chr[14]); return atoi(tim); /* The following functions actually evolve the molecules. */ int evolve(struct molecule *part, int boxsize) { int dleft; dleft = (*part).speed; while(dleft > 0) { switch( (*part).v) { case 1: evolvel(part, boxsize, &dleft); break; case 2: evolve2(part, boxsize, fedleft); Appendix A. Listing of Computational Program POINREC.C 144 break; case 3: evolve3(part, boxsize, ffidleft); break; case 4: evolve4(part., boxsize, fedleft); break; case 5: evolve5(part, boxsize, &dleft); break; case 6: evolve6(part, boxsize, &dleft); break; case 7: evolve7(part, boxsize, &dleft); break; case 8: evolve8(part, boxsize, fedleft); break; default : p r i n t f ( " E r r o r : value of v for p a r t i c l e */,d i s i n v a l i d . " e x i t ( l ) ; > > i f ( (*part).x < boxsize ) { (*part).side = 1; return 1; > else i f ( (*part).x >= boxsize ) { (*part).side =2; return 2; > else { printf("Trouble!!!!"); e x i t ( l ) ; > > void evolvel(struct molecule * p a r t l , int bsize, int * d l e f t p l ) { i f ( ((*partl).x < * d l e f t p l ) || ((bsize-1)-(*partl).y < * d l e f t p l ) ) { /* It w i l l bounce off a wall */ i f ( (*partl).x < ( ( b s i z e - l ) - ( * p a r t l ) . y ) ) { /* Bounces off l e f t wall */ *dleftpl-=(*partl).x; (*partl).y+=(*partl).x; (*partl).x = 0; (*partl).v=3; > else i f ( ( ( b s i z e - l ) - ( * p a r t l ) . y ) < (*partl).x ) { /* Bounces off top wall */ * d l e f t p l - = ( ( b s i z e - l ) - ( * p a r t l ) . y ) ; (*partl).x-=((bsize-l)-(*partl).y); (*parti).y = (bsize-1); (*partl).v=6; } else i f ( ( ( b s i z e - l ) - ( * p a r t l ) . y ) ==.(*partl).x ) { /* Bounces off top l e f t corner */ *dleftpl-=(*partl).x; (*partl).x = 0; (*partl).y = (bsize-1); Appendix A. Listing of Computational Program POINREC.C O p a r t l ) .v=8; } > else i f ( ( O p a r t l ) . x >= * d l e f t p l ) kk ( ( ( b s i z e - l ) - ( * p a r t l ) . y ) >= * d l e f t p l ) ) { /* Doesn't h i t a wall or just makes i t to a wall */ Opartl).x-=*dleftpl; O p a r t l ) .y+=*dleftpl; *dleftpl=0; /* If i t just h i t s a wall, direction w i l l change: */ i f ( ( O p a r t l ) . x== 0) kk ( O p a r t l ) . y == (bsize-1)) ) O p a r t l ) .v = 8; else i f ( ( O p a r t l ) . x == 0) kk ( O p a r t l ) . y < (bsize-1)) ) O p a r t l ) .v = 3; • else i f ( ( O p a r t l ) . x > 0) kk .(Opartl).y == (bsize-1)) ) O p a r t l ) .v = 6 ; } y void evolve2(struct molecule *part2, int bsize, int *dleftp2) { i f ( ((bsize-l)-Opart2) .y) < *dleftp2 ) { /* Bounces off top wall */ *dleftp2-=((bsize-1)-(*part2).y); Opart2).y=(bsize-l); (*part2).v=7; > else i f ( ((bsize-l)-(*part2).y) >= *dleftp2 ) { /* Doesn't reach wall or just makes i t to top wall */ Opart2).y+=*dleftp2; *dleftp2=0; i f ( Opart2).y == (bsize-1) ) Opart2) .v=7; > }. void evolve3(struct molecule *part3, int bsize, int *dleftp3) i f ( ((bsize+bsize-l-(*part3).x) < *dleftp3) II . ((bsize-l)-Opart3) .y < *dleftp3) ) { /* It w i l l bounce off a wall */ i f ( (bsize+bsize-1-Opart3).x) < (bsize-l-(*part3).y) ) { /* Bounces off right wall */ *dleftp3-=(bsize+bsize-l-(*part3).x); Opart3).y+=(bsize+bsize-l-(*part3).x); (*part3).x=(bsize+bsize-l); (*part3).v=l; > else i f ( (bsize-1-Opart3).y) < (bsize+bsize-1-Opart3).x) /* Bounces off top wall */ *dleftp3-=(bsize-l-Opart3) .y) ; Opart3).x+=(bsize-l-(*part3).y); Opart3) .y=(bsize-l); (*part3).v=8; Appendix A. Listing of Computational Program POINREC.C 146 > else i f ( (bsize-1-(*part3).y) == (bsize+bsize-1-(*part3).x) ) { /* Bounces off top right corner */ *dleftp3-=(bsize-l-(*part3).y); (*part3).x=(bsize+bsize-l); (*part3).y=(bsize-l); (*part3).v=6; } > else i f ( ( (bsize+bsize-l-(*part3).x) >= *dleftp3 ) && ( (bsize-1-(*part3).y) >= *dleftp3 ) ) { /* Doesn't h i t a wall or just makes i t to a wall */ (*part3).x+=*dleftp3; (*part3).y+=*dleftp3; *dleftp3=0; /* If i t just h i t s a wall, direction w i l l change: */ i f ( ((*part3).x == (bsize+bsize-1)) kk ((*part3).y == (bsize-1)) ) (*part3).v=6; else i f ( ((*part3).x == (bsize+bsize-1)) && ((*part3).y < (bsize-1)) ) (*part3).v=l; else i f ( ((*part3).x < (bsize+bsize-1)) kk ((*part3).y == (bsize-1)) ) (*part3).v=8; > > void evolve4(struct molecule *part4, int bsize, int *dleftp4) { i f ( (*part4);X < *dleftp4 ) { /* Bounces off l e f t wall */ *dleftp4-=(*part4).x; (*part4).x=0; (*part4).v=5; > else i f ( (*part4).x >= *dleftp4 ) {. /* Doesn't reach wall or just makes i t to l e f t wall */ (*part4).x-=*dleftp4; *dleftp4=0; i f ( (*part4).x == 0 ) (*part4).v=5; } > void evolve5(struct molecule *part5, int bsize, int *dleftp5) { i f ( (bsize+bsize-1)-(*part5).x < *dleftp5 ) { •/* Bounces, off right wall */ *dleftp5-=(bsize+bsize-l-(*part5).x); (*part5).x=(bsize+bsize-l); (*part5).v=4; > • else i f ( (bsize+bsize-l)-(*part5).x >= *dleftp5 ) .{ /* Doesn't reach wall or just makes i t to righ t wall */ (*part5).x+=*dleftp5; *dleftp5=0; Appendix A. Listing of Computational Program POINREC.C i f (, (*part5).x == (bsize+bsize-1) ) (*part5). v=4; > } • • . • • • • void evolve6(struct molecule *part6, int bsize, int *dleftp6) { • ' i f ( ((*part6).x < *dleftp6) II ((*part6).y < *dleftp6) ) { /* It w i l l bounce off a wall */ i f ( (*part6).x < (*part6).y ) .{ /* Bounces off l e f t wall */ *dleftp6-=(*part6).x; (*part6).y-=(*part6).x; (*part6).x=0; (*part6).v=8; } else i f ( (*part6).y < (*part6).x ) { /* Bounces off bottom wall */ *dleftp6-=(*part6).y; (*part6).x-=(*part6).y; (*part6).y=0; (*part6).v=l; > else i f ( (*part6).y == (*part6).x ) { /* Bounces off bottom l e f t corner */ *dleftp6-=(*part6).y; (*part6).x=0; (*part6).y=0; (*part6).v=3; • } . ' > else i f ( ((*part6).x >= *dleftp6) && ((*part6).y >= *dleftp6) ) { /* Doesn't h i t a wall or just makes i t to a wall */ (*part6>.x-=*dleftp6; (*part6).y-=*dleftp6; *dleftp6=0; /* I f . i t just h i t s a wall, direction w i l l change: */ i f ( ((*part6).x == 0) && ((*part6).y == 0) ) (*part6).v=3; else i f ( ((*part6).x == 0) kk ((*part6).y > 0) ) (*part6).v=8; else i f ( ((*part6).x > 0) kk ((*part6).y ==0) ) (*part6).v=l; } > void evolve7(struct molecule *part7, int bsize, int *dleftp7) i f ( (*part7).y < *dleftp7 ) { /* Bounces off bottom wall */ *dleftp7-=(*part7).y; (*part7).y=0; (*part7).v=2; > Appendix A. Listing of Computational Program POINREC.C 148 else i f ( (*part7).y >= *dleftp7 ) { /* Doesn't reach wall or just makes i t to bottom wall */ Opart7).y-=*dleftp7; *dleftp7=0; i f ( (*part7).y == 0 ) . (*part7).v=2; > void evolve8(struct molecule *part8, int bsize, int *dleftp8) { i f ( ((bsize+bsize-l-(*part8).x) < *dleftp8) || ((*part8).y < *dleftp8) ) { /* It w i l l bounce off a wall */ i f ( (bsize+bsize-l-(*part8).x) < (*part8).y ) { /* Bounces off right wall */ *dlef.tp8-=(bsize+bsize-l-(*part8) .x) ; (*part8).y-=(bsize+bsize-l-(*part8).x); (*part8).x=(bsize+bsize-l); (*part8).v=6; > else i f ( (*part8).y < (bsize+bsize-1-(*part8).x) ) { /* Bounces off bottom wall */ *dleftp8-=(*part8).y; (*part8).x+=(*part8).y; (*part8).y=0; (*part8).v=3; > else i f ( (*part8).y == (bsize+bsize-l-(*part8) .x) ) {. /* Bounces off bottom right corner */ *dleftp8-=(*part8).y; (*part8).x=(bsize+bsize-l); (*part8).y=0; (*part8).v=l; > > else i f ( ((bsize+bsize-l-(*part8).x) >= *dleftp8) kk ((*part8).y >= *dleftp8) ) { " /* Doesn't h i t a wall or just makes i t to a wall */ (*part8).x+=*dleftp8; (*part8).y-=*dleftp8; *dleftp8=0; /* If i t just h i t s a wall, direction w i l l change: */ i f ( ((*part8).x == (bsize+bsize-1)) kk ((*part8).y == 0) ) (*part8).v=l; else i f ( ((*part8). .x == (bsize+bsize-1)) kk ((*part8).y > 0) ) (*part8).v=6; else i f ( ((*part8).x < (bsize+bsize-1)) kk ((*part8).y == 0) ) (*part8).v=3; > } /* The following function takes care of the interactions (or bounces). */ void bounce(struct molecule *partA, struct molecule *partB) int tempv =0; Appendix A. Listing of Computational Program POINREC.C int i = 0; enum boolean {no, yes}; enum boolean found = yes; enum boolean toggle = yes; struct molecule *tempp; while( (toggle==yes) kk (i<2) ) { toggle=no; if ( OpartA) .v==l ) { switch( (*partB).v ) { case 2: OpartA). v=6; OpartB) .v=5; break; case 4: OpartA). v=3; OpartB) .v=7; break; case 5: OpartA) .v=6; OpartB) .v=7; break; case 7: OpartA) .v=3; (*partB).v=5; break; default: found=no; . } } else if ( OpartA) . v==3 ) { switch( OpartB).v ) { case 2: OpartA) • v=8; OpartB).v=4; break; case 4: (*partA).v=8; OpartB) .v=7; break; case 5: OpartA) .v=l; OpartB) .v=7; break; case 7: OpartA).v=l; OpartB) .v=4; break; default: found=no; } } else if ( OpartA) .v==6 ) { switch( (*partB).v ) { case 2: OpartA) .v=8; (*partB).v=5; break; case 4: (*partA).v=8; (*partB).v=2; break; case 5: (*partA).v=l; OpartB) .v=2; break; case 7: OpartA) • v=l; Appendix A. Listing of Computational Program POINREC.C Opar tB) .v=5; break; de fau l t : found=no; > > e lse i f ( Opar tA) .v==8 ) { switch( O p a r t B ) . v ) { case 2: Opar tA) . v=6; OpartB).v=4; break; case 4: O p a r t A ) . v=3; Opar tB) .v=2; break; case 5: Opar tA) .v=6; O p a r t B ) . v=2; break; case 7: Opar tA) .v=3; OpartB).v=4; break; de fau l t : found=no; } } e lse { toggle=yes; tempp=partA; partA=partB; partB=tempp; i++; > } • i f ( (found==no) | | (i==2) ) { tempv=(*partA).v; (*par tA) .v = O p a r t B ) . v ; O p a r t B ) . v = tempv; > > A p p e n d i x B L i s t i ng of Computa t iona l P rog ram S E Q U E N C E . C The following is a listing of the C program used for calculating the elements of the sequence discussed in Chapter 4. /* SEQUENCE.C */ #include<stdio.h> #include<math.h> #include<string.h> #include<ctype.h> #include<stdlib.h> /* Program to count the occurence of each number (1-9) i n the sequence consisting of the f i r s t d i g i t of 2~n */ void output(int s t a r t , int *occur, double *act, FILE *rfp[] , FILE * f p [ ] ) ; void i n i t ( i n t *array, int array_size, int.value); main(int argc, char *argv[]) { FILE * f i l e s [ 9 ] ; FILE * r f i l e s [ 9 ] ; register int i , 1; enum boolean {no, yes}; int found; int c [10]; double n l , j ; double actual [10]; char mode[1]; int s t a r t , stop, stepsize; int ans; char namein[77], name[80]; i f (!( (argc==6) kk ( (start=atoi(argv[1])) > 0) kk ( (stop=atoi(argv[2])) > 0) kk ( (stepsize=atoi(argv[3])) > 0) kk ( strlen(strncpy(namein,argv[4],77)) < 78 ) kk ( ((mode[0]=*argv[5])=='a') || ((mode[0] =*argv[5] )=='w>) ) ).) { printf("Run the program using the following syntax please:\n\n"); printf("sequence [start] [stop] [stepsize] [filename] [mode]\n\n"); printf("where:\n\n"); p r i n t f ( " s t a r t specifies the sequence element at which sequence w i l l " ) ; printf("begin execution.\n"); 151 Appendix B. Listing of Computational Program'SEQUENCE.C 152 printf("stop specifies the sequence element after which sequence w i l l . " ) ; printf("terminate execution.\n"); pr i n t f ( " s t e p s i z e specifies the intervals at which output w i l l be " ) ; pr i n t f ( " w r i t t e n \n to the data f i l e s , s t a r t i n g at start and " ) ; printf("ending at stop.\n"); printf("mode must be either 'w' or 'a'. 'w' specifies that sequence " ) ; printf("should open \n new data f i l e s while 'a' specifies that " ) ; printf("data should be appended \n to already existing data " ) ; p r i n t f ( " f i l e s ( i f sequence i s started with \n mode 'a 5 and the " ) ; printf("appropriate data f i l e s do not exist i n the \n " ) ; prin t f ( " f i l e s y s t e m , an error w i l l occur and sequence w i l l terminate \n"); p r i n t f ( " execution).\n"); printf("filename specifies the names of the output f i l e s , and may " ) ; printf("not be more \n than 77 characters long. Output " ) ; p r i n t f ( " w i l l be saved to the f i l e s \n. filenamel, " ) ; printf("filename9 and filenamel.r, rfilename9_r.\n\n"); /* Put in error checking f o r the following user inputs */ . printf("\nEnter a value for [ s t a r t ] : " ) ; scanf ("'/,d",&start); i printf("\nEnter a value for [stop]: '"); scanf 0"/.d" ,&stop); printf("\nEnter a value f o r [stepsize]: " ) ; scanf 0"/.d" ,&stepsize); printf("\nEnter a filename: " ) ; scanf C"/,s77",namein); printf("\nEnter a filemode: " ) ; scanf ("'/.si'.',mode) ; /* mode[0]=getchar(); */ > /* Open the 18 output f i l e s , checking for error upon opening each one */ for (i=0; i<9;i++) { strncpy(name,namein,77); sprintf (name, '"/.s'/d" ,name,i+l); i f ( ( f i l e s [i]=f open (name, mode)) == NULL )'•{ pr i n t f ( " E r r o r : cannot open f i l e . Try again."); e x i t ( l ) ; } strcat(name,"_r"); i f ( (rfiles[i]=fopen(name.mode)) == NULL ) { printf("\nError: cannot open f i l e . Try again.\n"); e x i t ( l ) ; } > /* Any previous data? */ i f (start > 1 H p r i n t f ("You have specified start = '/,d. Do you wish to enter ", s t a r t ) ; p rintf("data f o r \nthe frequency of occurence of each d i g i t up to and " ) ; p r i n t f ("including the \n'/,dth element of the sequence. If you do " , s t a r t - l ) ; printf.("not have any previous data, \nthe f i r s t '/.d elements of " , s t a r t - l ) ; p r i n t f ( " t h e sequence w i l l be ignored.\n"); printf("Do you have any previous data to enter (<y>es, <n>o, <q>uit)? " ) ; ans=getchar(); while ( (tolower( ans)!—'y') && (tolower(ans)!='n') && (tolower(ans)! ='q') )•[ Appendix B. Listing of Computational Program SEQUENCE. C printf("\nPlease answer <y>, <n>, or <q>: " ) ; ans=getchar(); > i f ( tolower(ans) == 'q' ) exit(O); i f ( tolower(ans) == 'n' ) { init(c,10,0); c[0]=start; . } i f ( tolower(ans) == 'y' ) { printf("\nEnter your data: How many elements of the sequence were ") printf("counted? " ) ; scanf ('7,d",&c[0]) ; printf("Enter the frequency of occurence for each d i g i t . \ n " ) ; for ( i = .1; i < 10; ++i) { pr i n t f ("\n 7.d: ", i ) ; scanf ("'/.d",&c[i]); > start=0; /* Check to make sure inputted data i s v a l i d */ i f ( c[0] != ( c[l]+c[2]+c[3]+cC4]+c[5]+c[6]+c[7]+c[8]+c[9] ) ) { printf("\nThere i s a problem with your data. Please check i t " ) ; printf("and t r y again.\n"); exit(O); > > > • else { /* i n i t i a l i z e array to zero */ init(c,10,0); start=0; } /* Prob. of occurence of each d i g i t i n sequence i s P(k) = loglO[(k+l)/k] actual[0]=0.0; for ( i =1; i < 10; ++i) actual[i]=logi0( (i+1.0)/(double)(i) ); /* Here's the part that actually computes the elements of the sequence * while (c[0] < stop) { . " ' . ' . for (1=0; 1 < stepsize; ++1){ found=no; nl=c[0]*logl0(2.0); nl=n l - f l o o r ( n l ) ; for (i=0; i<10 kk found==no; ++i){ j=(double)(i); i f ( (nl < logl0(j+1.0)) kk (nl >= logl0(j))){ ++c[i]; found=yes; } . } • ++c[0]; > /* This l i n e makes stepsize go up by a factor of 10 every decade */ /* i f ( c[0]/stepsize == 10.0 ) stepsize=stepsize*10; Appendix B. Listing of Computational Program SEQUENCE. C 154 */ output(start, c, actual,- r f i l e s , f i l e s ) ; . • } > •' ' /* function to i n i t i a l i z e an integer array of size array_size to value */ void i n i t ( i n t *array, i n t array_size, int value) { register int i ; for ( i = 0; i ,.< array_size; ++r) array[i]=value; > /* function to ''output', result's to a set of output f i l e s •*/ void output'(int' s t a r t , int *occur, double *act, FILE * r f p [ ] , FILE *fp[]) register int i ; • double prob; for (i=0; i<9; i++){ f p r i n t f ( r f p [ i ] ,"7,d y.9d 7,2.16f \n" , (occur[0]-start) , occur[i+1] , ( prob=(double)(occur[i+1])/(double)((occur[0]-start)) ) ); f p r i n t f (fp[i] , "7,d 7.2.16f \n" , (occur [0]-start), act [i+1] -prob); } 

Cite

Citation Scheme:

        

Citations by CSL (citeproc-js)

Usage Statistics

Share

Embed

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

Comment

Related Items