High Order WENO Scheme for Computational Cosmology by Ishani Roy B.S., University of Colorado; Boulder, CO, 2003 M.S., University of Massachusetts; Amherst, MA, 2006 M.S., Brown University; Providence, RI, 2007 A dissertation submitted in partial fulfillment of the requirements for the degree of Doctor of Philosophy in The Division of Applied Mathematics at Brown University PROVIDENCE, RHODE ISLAND May 2010 c Copyright 2010 by Ishani Roy This dissertation by Ishani Roy is accepted in its present form by The Division of Applied Mathematics as satisfying the dissertation requirement for the degree of Doctor of Philosophy. Date Chi-Wang Shu, Ph.D., Advisor Recommended to the Graduate Council Date Johnny Guzm´an, Ph.D., Reader Date Jan S. Hesthaven, Ph.D., Reader Approved by the Graduate Council Date Sheila Bonde, Dean of the Graduate School iii Vitae Personal Information Born: 15th October 1980, Kolkata, India Address: Division of Applied Mathematics Phone: (413) 687-2354 Brown University, Providence, RI E-mail: iroy@dam.brown.edu Education Brown University • PhD in Applied Mathematics May 2010 expected Thesis - High Order WENO Method in Computational Cosmology Thesis Adviser - Professor Chi-Wang Shu • Masters in Applied Mathematics May 2007 University of Massachusetts at Amherst Masters in Applied Mathematics May 2006 University of Colorado at Boulder Bachelors in Applied Mathematics December 2003 Research Interests • High order numerical methods for conservation laws, especially discontinuous Galerkin method, weighted ENO method and spectral methods. • Application of WENO method in Computational Cosmology • Numerical computation of nonlinear φ4 field theory Research Experience Brown University, Division of Applied Mathematics 2007 - Present WENO Methods in Computational Cosmology iv Los Alamos National Laboratory, Center for Nonlinear Studies Summer 05,06 University of Massachusetts at Amherst 2005 - 2006 Numerical Solution of φ4 Field Theory University of Colorado at Boulder, Physics Fall 2003 Analysis of Rare Decay(s) of B Meson University of Colorado at Boulder, Applied Mathematics 2002 - 2003 Multiresolution Analysis and Wavelets Publications/Preprints 1. I. Roy, C-W. Shu, L-Z. Fang, Resonant Scattering and Lyα Damping Wing of Neutral Hydrogen Halos, to appear in The Astrophysical Journal 2. I. Roy, W. Xu, J-M. Qiu, C-W. Shu and L-Z. Fang, Wouthuysen-Field coupling in 21 cm region around high redshift sources, The Astrophysical Journal 703: 1992-2003 (2009) 3. I. Roy, W. Xu, J-M. Qiu, C-W. Shu and L-Z. Fang, Time Evolution of Wouthuysen- Field coupling, The Astrophysical Journal 694: 1121-1130 (2009) 4. I. Roy, J-M. Qiu, C-W. Shu and L-Z. Fang, A WENO algorithm for radia- tive transfer with resonant scattering: the time scale of the Wouthuysen-Field Coupling, New Astronomy 14: 513-520 (2009) 5. I. Roy, S. Dmitriev, P. G. Kevrekidis, A. Saxena, Comparative Study of Dif- ferent Discretizations of the φ4 Model, Physical Review E 76: 026601-026615 (2007) Visiting Researcher Indian Institute of Science, Bangalore Summer 2009 Department of Aerospace Engineering University Pierre Marie Curie, Paris Summer 2008 Laboratoire Jacques-Louis Lions v Presentations Brown University - Laboratoire Jacques-Louis Lions, joint Seminar February 2010 Metamaterials: Applications, Analysis and Modeling workshop January 2010 Institute for Pure and Applied Mathematics, UCLA Joint Mathematics Meetings, San Francisco January 2010 AMS-SIAM Special Session on Mathematics of Computation International Conference on Advances in Scientific Computing, December 2009 Brown University, Providence SIAM Annual Meeting, Denver July 2009 Mini-symposium on Advanced Numerical Methods for Kinetic Equations University of New Mexico, Los Alamos July 2006 Los Alamos Student Symposium Los Alamos National Laboratory July 2005 and 2006 Center for Nonlinear Studies Colloquium University of Colorado at Boulder, Department of Applied Mathematics July 2002 4th Annual VIGRE Symposium Teaching Experience Brown University 2007 - 2008 Graduate Teaching Assistant University of Massachusetts at Amherst 2004 - 2006 Graduate Teaching Assistant University of Colorado at Boulder • Science, Technology, Engineering and Mathematics, Teacher Preparation Fall 2003 • NSF supported Teaching Assistant in Calculus for Engineers Upward Bound Program, Teaching Assistant Summer 2003 Taught pre-calculus and calculus to gifted minority high school students vi Other Activities Challenges of Mathematical Modeling and Computation in Finance Brown University, Providence September 2009 • Received International Colloquium Grant - Office of International Affairs • Organized two day colloquium with speakers from leading financial institutes and academia Awards and Honors Graduate Student Fellowship, Brown University 2006 - 2007 Undergraduate Research Opportunity Program, Scholarship Summer 2003 J Ranald Fox Scholarship, University of Colorado 2002 - 2003 Md. Zanborie Scholarship, University of Colorado 2001 - 2002 Interests Indian classical dance, Chhandika Institute of Kathak Dance, Boston vii Acknowledgements This is after all a dissertation in Mathematics. It is only natural that I start with acknowledging the four very fun, and in some ways mathematical, years I have spent at the Division of Applied Mathematics at Brown University. First and foremost I would thank my adviser Professor Chi-Wang Shu for his help and support throughout the last four years. I was very fortunate when he accepted my tremulous approach in asking him to be my adviser with a hearty laugh. From then on I have learned a lot from him along the way, not just from his scientific knowledge, clear explanations but also from his unwavering patience towards all my questions and concerns. I would like to acknowledge the contributions of the members of my committee, Professor Jan Hesthaven and Professor Johnny Guzman who have given me useful suggestions which has made my thesis better in many ways. The division has always felt like an extended family mostly due to all our Pro- fessors. It would not be complete if I do not acknowledge the influence of Professor David Gottlieb. I was one of the fortunate few to be his student, to have a chance to listen to his lectures as well as numerous interesting stories about mathematics celebrities and also probably be the only person in the division to receive a chocolate truffle from him on the day of the dreaded preliminary examinations. I am glad to know his daughter Professor Sigal Gottlieb whose ever smiling face and encouraging words have always motivated me to take on the daunting task of becoming a female academician. I would like to acknowledge the help and support of our division chair viii Professor Paul Dupuis who took time to help two graduate students in their naive notions of organizing an international colloquium, who had very little idea of orga- nizing let alone know how to apply for a grant. Professor Dupuis’ advice and the division’s support made it possible and also successful in many ways. On a similar note, I would acknowledge the thoughts and advice of Professor Jerome Stein who very kindly took it on his hand to guide us along the way to organize this colloquium. I am glad to have received the support of Professor Hui Wang, Professor Constantine Dafermos and Professor Wendall Fleming. I am also thankful to Professor Govind Menon, who has always had time to talk with me about my future plans. I am glad to know Professor Jennifer Ryan whose sound advice, in matters personal and pro- fessional, I treasure in every step. Finally I would like to thank the ever helpful Jean Radican, Laura Leddy, Roselyn Winterbottom, Stephanie Han and Janice D’Amico, who all make the division so friendly and special. I can safely say that I was not one who always knew that I would one day pursue a doctoral degree in mathematics. Mathematics slowly sort of grew on me and for that I would have to thank my professors at UMass Amherst, most importantly Professor Panayotis Kevrekidis. I was able to work as a summer intern in the Los Alamos National Laboratory due to his reference and advice. The project I worked on there lead me to my decision to continue in a research career and to apply for doctoral studies. I attribute a big part of my interest in numerical analysis and pursuit of higher degree in the subject to Professor Kevrekidis mentoring. This would be a good opportunity for me to thank my supervisor Dr. Avadh Saxena at Los Alamos National Laboratory, who has not only taught me a lot of science but also taught me how to explain science and the secrets of a good science presentation. In the beginning, Astrophysics seemed like an insurmountable subject when I was initiated to it as my dissertation topic, shrouded with mysterious concepts like black ix holes and cosmic energy. I am grateful to my collaborator Professor Li-Zhi Fang for making those mysteries more understandable as time went on. I would also thank my other collaborator and friend Professor Jing-Mei Qiu who has always answered every question with tremendous patience, however trivial. I would like to thank my mother Jaya Roy who has been my hero since I was a little girl. I have always been fascinated with my mother’s unending energy, and this is what motivates me take on hundred tasks at one time. She would have been proud of me even if I had been a dancer or a mountaineer and not pursued a PhD in mathematics. I am very thankful towards my parents who have so much faith in me. When I think back, it is because of my father Prasun Roy without whom I wouldn’t have even traveled from beautiful warm Calcutta to freezing Colorado “in the pursuit of mathematics”. Thank you for saying I have fulfilled your own dreams. I would like to acknowledge the influence of my sister Biyas Roy, she is the smartest woman I know and I am lucky that I have had the chance to look up to her all my life. I want to thank my most favorite mathematician of all, Saugata Bandyopadhyay, who has counseled me in many things from career choices to the silliest things. I also cherish the guidance of my uncle Professor Atanu Sinha and the occasional arguments that we have in almost every subject. Where else will I get a chance to thank the few most precious family members who are also the some of the strongest and kindest people I have known - my very loving grandmother Pratibha Rani Dey and my aunts Maitri Mitra and Manjusree Dutta. I would like to thank my other set of parents, Nilima Chitnis and Ravindra Chitnis for being so proud of me at every little thing I do. Lastly and very fondly, my grandfather Sripad Vidhwans for giving me one of the best compliments I have ever received. I am glad that during the four long years of my PhD studies I have been able to share the joys (and agonies) with many of my close friends. I can’t thank enough x my childhood friend, from our kindergarten days in Kolkata, Sayak Mukherjee for patiently listening to all my silly troubles. Also, Sourish Basu and Sohini Dasgupta whose words and encouragements have always been very spacial to me. I would like to thank my peers in the division at Brown especially Akil Narayan and An- dreas Kl¨ ockner who helped me with smallest real analysis questions and the biggest “my-computer-just-crashed” ones. I am fortunate to have been friends with Kenny Chowdhary and Chia Ying Lee. We suffered together through the dreaded analysis and probability in our first year and survived, couldn’t have asked for better friends then them. I am thankful to my group members Yinhua Xia, Xiangxiong Zhang, Sirui Tan and especially Yanlai Chen who has very methodically looked over my work and helped me at every instance. I want to thank my former office mate and friend Wei Wang for many interesting scientific and non scientific discussions including tu- torials in Mandarin and my present office mate Xinghui Zhong with whom I have continued this friendly tradition and who has also helped me numerous times during my defense and thesis preparation. I am thankful to the talented Rusha Sopariwala for preparing the beautiful plots which make my boring math talks a little bit more bearable. Last but definitely not the least, I thank my valuable friendship with Catherine Benincasa and Rosalyn Benincasa for many wonderful moments I have spent with them. No, I didn’t forget him :). I am not sure how to thank him, so here’s to my best friend, Shailesh Chitnis, (who also happens to be my husband). Thank you for putting up with me for the last ten years. xi Abstract of “ High Order WENO Scheme for Computational Cosmology ” by Ishani Roy, Ph.D., Brown University, May 2010 This doctoral dissertation is concerned with the formulation and application of a high order accurate numerical algorithm suitable in solving complex multi dimensional equations and the application of this algorithm to a problem in Astrophysics. The algorithm is designed with the aim of resolving solutions of partial differential equations with sharp fronts propagating with time. This high order accurate class of numerical technique is called a Weighted Essentially Non Oscillatory (WENO) method and is well suited for shock capturing in solving conservation laws. The numerical approximation method, in the algorithm, is coupled with high order time marching as well as integration techniques designed to reduce computational cost. This numerical algorithm is used in several applications in computational cosmol- ogy to help understand questions about certain physical phenomena which occurred during the formation and evolution of first generation stars. The thesis is divided broadly in terms of the algorithm and its application to the different galactic processes. The first chapter deals with the astrophysical problem and offers an introduction to the numerical algorithm. In chapter 2 we outline the mathematical model and the various functions and parameters associated with the model. We also give a brief description of the relevant physical phenomena and the conservation laws associated with them. In chapter 3, we give a detailed description of the higher order algorithm and its formulation. We also highlight the special techniques incorporated in the algorithm in order to make it more suitable for handling cases which are computationally intensive. In the later chapters, 4-7, we explore in detail the physical processes and the different applications of our numerical scheme. We calculate different results such as the time scale of a temperature coupling mechanism, radiation and intensity changes etc. Different tests are also performed to illustrate the stability and accuracy of this algorithm. Chapter 8 is concerned with ongoing work and the future research that will be carried out in this topic. xiii Contents Vitae iv Acknowledgments viii 1 Introduction 1 2 Radiative transfer of photons in early universe 6 2.1 Radiative transfer with resonant scattering in frequency space . . . . 7 2.1.1 Re-scaling the equations . . . . . . . . . . . . . . . . . . . . . 9 2.2 Radiative transfer equation in the 21 cm region . . . . . . . . . . . . 11 2.2.1 Eddington approximation . . . . . . . . . . . . . . . . . . . . 13 3 The high order accurate algorithm 16 3.1 High order WENO based scheme for the RT equation in frequency space 17 3.1.1 Computational domain and computational mesh . . . . . . . . 17 3.1.2 Algorithm of the spatial derivative . . . . . . . . . . . . . . . 18 3.1.3 Time evolution using TVD Runge-Kutta . . . . . . . . . . . . 20 3.2 Numerical integration: Higher order scheme and an O(N ) algorithm . 20 3.2.1 High order numerical integration . . . . . . . . . . . . . . . . 21 3.2.2 O(N ) algorithm for numerical integration . . . . . . . . . . . 21 3.3 Algorithm for the RT equation in a spherical halo . . . . . . . . . . . 25 3.3.1 The WENO algorithm: approximations to the spatial derivatives 25 3.3.2 Implementation of the boundary condition . . . . . . . . . . . 27 3.3.3 Adaptive mesh procedure for non-uniform grid . . . . . . . . . 27 3.4 Conservation and normalization techniques . . . . . . . . . . . . . . . 28 3.4.1 Test with the conservation of the photon number . . . . . . . 29 3.4.2 Normalization of the re-distribution term . . . . . . . . . . . . 29 4 Radiative transfer in the frequency space 31 4.1 Radiative transfer (RT) equation in the frequency space . . . . . . . . 32 4.1.1 RT equation in static background . . . . . . . . . . . . . . . . 32 xiii 4.1.2 Expanding background . . . . . . . . . . . . . . . . . . . . . 34 4.1.3 The effect of recoil of the atom . . . . . . . . . . . . . . . . . 36 4.2 The Wouthuysen-Field Coupling . . . . . . . . . . . . . . . . . . . . . 38 4.3 Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 42 5 Time scale of Wouthuysen Field coupling and intensity calculation 43 5.1 Wouthuysen-Field coupling in a static background . . . . . . . . . . . 45 5.2 Wouthuysen-Field coupling in an expanding background . . . . . . . 48 5.2.1 Width of the local Boltzmann distribution . . . . . . . . . . . 48 5.2.2 Photon source and W-F coupling . . . . . . . . . . . . . . . . 52 5.2.3 Intensity . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 54 5.2.4 Time scales of the frequency shift . . . . . . . . . . . . . . . . 54 5.2.5 Bounce back mechanism . . . . . . . . . . . . . . . . . . . . . 57 5.2.6 Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 58 6 Radiative transfer in the 21 cm region 60 6.1 Basic properties of the 21 cm region . . . . . . . . . . . . . . . . . . . 61 6.1.1 Problems with the W-F coupling . . . . . . . . . . . . . . . . 63 6.2 Radiative transfer (RT) of Lyα photons in the 21 cm regions . . . . . 65 6.2.1 Diffusion in the physical space . . . . . . . . . . . . . . . . . . 66 6.2.2 Resonant photon restoration and W-F coupling onset . . . . . 68 6.2.3 Effect of injected photons . . . . . . . . . . . . . . . . . . . . 72 6.2.4 Effect of the Hubble redshift . . . . . . . . . . . . . . . . . . . 72 6.3 Discussion and conclusions . . . . . . . . . . . . . . . . . . . . . . . . 73 7 Radiation from neutral Hydrogen and flash sources 77 7.1 Solutions of time-independent sources . . . . . . . . . . . . . . . . . . 78 7.1.1 Optically thick halos . . . . . . . . . . . . . . . . . . . . . . . 78 7.1.2 Radiative transfer equation in spherical halo . . . . . . . . . . 80 7.2 Solutions of time-independent sources . . . . . . . . . . . . . . . . . . 82 7.2.1 Time scale of local thermalization . . . . . . . . . . . . . . . . 82 7.2.2 Restoring ν0 photons . . . . . . . . . . . . . . . . . . . . . . . 83 7.2.3 Two peaks in the flux profile . . . . . . . . . . . . . . . . . . . 84 7.2.4 Absorption spectrum . . . . . . . . . . . . . . . . . . . . . . . 86 7.2.5 The estimation of the HI column density . . . . . . . . . . . . 87 7.3 Modeling of flash sources . . . . . . . . . . . . . . . . . . . . . . . . . 89 7.3.1 Frequency-dependence of Lyα transfer . . . . . . . . . . . . . 89 7.3.2 The evolution of time duration of a flash . . . . . . . . . . . . 89 7.4 Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 92 8 Ongoing work 93 9 Conclusion 96 xiv List of Figures 1.1 Big Bang timeline . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2 1.2 Halo of the first star . . . . . . . . . . . . . . . . . . . . . . . . . . . 3 4.1 Static solutions (γ = 0) for pure Doppler redistribution (a = 0) of equation(2.12), in which C = 1 and J(x, 0) = 0. The analytical solutions are shown by dashed lines, while the numerical results are shown by solid lines. . . . . . . . . . . . . . . . . . . . . . . . . . . . 34 4.2 Static solutions (γ = 0) for pure Doppler redistribution (a = 0) of equation(2.12), in which C = 0 and J(x, 0) = π −1/2 exp(−x2 ). The analytical solutions are shown by dashed line, while the numerical results are shown by solid lines. . . . . . . . . . . . . . . . . . . . . . 35 4.3 Solutions of J(x, τ ) of equation(4.4) with initial condition J(x, 0) = 0. The parameter γ is taken to be 10−1 (left panel) and 10−3 (right). The analytical solutions equation(4.5) are shown by dashed lines . . . . . 36 4.4 Solutions of J(x, τ ) of equation (2.12) with initial condition J(x, 0) = 0. The parameter γ is taken to be 10−1 (left) and 10−3 (right). . . . . 36 4.5 Solutions of J(x, τ ) of equation(2.12) with redistribution function equation (2.5), and γ = 0. . . . . . . . . . . . . . . . . . . . . . . . . 37 4.6 F (τ ) as function of τ for solutions shown in figure 4.5 (left) and figure 4.1 (right) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 38 4.7 Solutions of J(x, τ ) of equation (2.12) with γ = 10−5 (left panel) and 10−7 (right panel). . . . . . . . . . . . . . . . . . . . . . . . . . . . . 40 4.8 F (τ ) as function of τ for solutions of a = 0 and γ = 10−5 (left) and 10−7 (right). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 41 5.1 WENO numerical solutions of equation (2.12) with γ = 0 and b = 0.03. The source is taken to be S = φ(x) and initial condition J(x, 0) = 0. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 45 5.2 WENO numerical solutions of equation (2.12) with γ = 0 and b = 0.03 (solid), 0.015 (dashed), 0.0079 (dot-dot-dashed) and 0.0025 (dot- dashed). The source is taken to be S = φ(x) and initial condition J(x, 0) = 0. The bottom panel is a blow-up of the dashed square of the top panel. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 47 5.3 B(τ ) vs. τ for solutions in figure 5.2 with b = 0.03 (solid), 0.015 (dashed), 0.0079 (dot-dot-dashed) and 0.0025 (dot-dashed). . . . . . . 48 xv 5.4 WENO numerical solutions of equation (2.12) with b = 0.03 and γ = 10−3 (two top panels) and γ = 10−5 (two bottom panels). The source is taken to be S = φ(x) and initial condition J(x, 0) = 0. The right panels are blow-up of the dashed square of the corresponded left panel 50 5.5 WENO numerical solutions of equation (2.12) with γ = 10−3 and b = 0.0079 (left), 0.015 (middle), 0.03 (right). The source is taken to be S = φ(x − 10) and initial condition J(x, 0) = 0. The bottom panels are the blow-up of the dashed square of the top panels. . . . . 51 5.6 WENO numerical solutions of equation (2.12) with γ = 1 (left) and γ = 10−1 (right). Parameter b = 0.0079 (dot-dot-dashed), 0.015 (dashed), 0.03 (solid). The source is S = φ(x − 10) and initial con- dition J(x, 0) = 0. The bottom panels are the blow-up of the dashed square of the top panels. . . . . . . . . . . . . . . . . . . . . . . . . . 52 5.7 WENO numerical solutions of equation .(2.12) with γ = 10−3 and b = 0.0079 (dot-dot-dashed), 0.015 (dashed), 0.03 (solid). The source is given by equation (5.11). The right panel is a blow-up of the dashed square of left panel. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 53 5.8 WENO numerical solutions of equation (2.12) with γ = 10−2 (left) and 10−4 (right) with source S(x) = φ(x − 10) and b = 0.0079 (dot- dot-dashed), 0.015 (dashed), 0.03 (solid). The bottom panels are the zoom in of the dashed square of top panels. . . . . . . . . . . . . . . 55 5.9 Solutions J(τ, x) of equation (2.12) with the re-distribution function equation (2.7). Parameters b = 0, and γ = 10−6 . . . . . . . . . . . . . 56 5.10 x¯ vs. γτ for the solution shown in figure 5.9 (solid) and the solution of equation (5.15) (dashed). . . . . . . . . . . . . . . . . . . . . . . . 57 6.1 The solution f (η, r) of equation (6.9) and (6.10). The straight line is f = S0 = 106 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 67 6.2 The flux f (η, r, x) of the solutions of equation (6.7) and (6.8) at r = 50 (left), 100 (middle) and 500 (right). The frequency profile of the √ 2 source is φg (x) = (1/ π)e−x . . . . . . . . . . . . . . . . . . . . . . 68 6.3 The flux f (η, r, x) of the solutions of equation (6.7) and (6.8) at r = 50 (left), 100 (middle) p and−2x 500 (right). The frequency profile of the source 2 is φs (x) = (1/ π/2)e . . . . . . . . . . . . . . . . . . . . . . . . 69 6.4 The mean intensity j(η, r, x) of equation (6.7) and (6.8) at r = 50 (left), 100 (middle) and 500 (right). . . . . . . . . . . . . . . . . . . . 70 6.5 The mean intensity j(η, r, x) (left) and flux f (η, r, x) (right) of the solution of equation (6.7) and (6.8) at r = 500. The recoil parameter is taken to be b = 0.3. . . . . . . . . . . . . . . . . . . . . . . . . . . 71 6.6 The mean intensity f (η, r, x) and flux j(η, r, x) of equation (6.7) and (6.8) at r = 50 (top), 100 (middle) and 500 (bottom), b = 0.3. The frequency profile of the source is φs (x) = φg (x − 3) . . . . . . . . . . 75 6.7 The solutions f (η, r, x) and j(η, r, x) of equation (6.7) and (6.8) at r = 50 (top), 100 (middle) and 500 (bottom). γ = 10−3 . The parameters are the same as those in figures 6.2 and 6.4, except that γ = 10−3 . . . 76 xvi 7.1 Top two panels are the solutions of j(η, r, x) (left) and f (η, r, x) (right) of equations (7.6) and (7.7) at r = R = 102 and η = 200, 300, 500 with the boundary condition equation (7.11). Bottom two panels are the solutions without the boundary condition equation (7.11). The source is taken to be (7.13 )with S0 = 1 . The parameters are taken to be a = 10−3 and γ = 0−5 . . . . . . . . . . . . . . . . . . . . . . . . 82 7.2 Solutions of j(η, r, x) (left) and f (η, r, x) (right) of (7.8) and (7.9) at r = 102 and η = 200, 300, 500. The source and all parameters are the same as in figure 7.1. . . . . . . . . . . . . . . . . . . . . . . . . . . 83 7.3 Flux f (η, r, x) at r = 102 at radius r = 102 and time η = p 200, 300 and 2 500. The frequency profile of the source is φs (x) = (1/ π/2)e−2x . Other parameters are the same as in figure 7.1 . . . . . . . . . . . . . 84 7.4 The positions of the peaks |x+ | as a function of (ar)1/3 for the solutions of equations (7.6) and (7.7) with a = 10−2 . Other parameters are the same as in figure 7.1. The dashed line is for x± = ±(ar)1/3 . . . . . . . 85 7.5 Solutions of j(η, r, x) (left) and f (η, r, x) (right) of equations (7.6) and (7.7) at r = 102 and η = 200, 300 and 500. The source is given by 7.14. The parameters are a = 10−3 and γ = 10−5 . . . . . . . . . . . . 87 7.6 Red damping wing (a) the DLA model τ0 = 4 × 103 (dotdot dashed line), 6 × 103 (dot dashed line) and 104 (dashed line); (b) the solution f (x) of (7.15) for τ0 = 104 (solid line). For both (a) and (b), the parameter a = 10−2 . . . . . . . . . . . . . . . . . . . . . . . . . . . . 88 7.7 The profiles of j(η, r, x) and f (η, r, x) with time-dependent source [eq.(7.14)]. Top panel: j(η, r, x) (left) and f (η, r, x) (right) of η0 = 100 at r = 102 . Bottom panel: η0 = 500 at r = 103 . The parameters a = 10−2 and γ = 10−5 . . . . . . . . . . . . . . . . . . . . . . . . . . . 90 7.8 The time-dependence of f (η, r, x) at x = 0. Top panels: f (η, r, 0) at r = 102 (left) and r = 103 (right) for flash source with η0 = 50. Bottom panel: f (η, r, 0) at r = 103 of flash source with η0 = 100 (left) and η0 = 500 (right). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 91 xvii Chapter One Introduction 2 Figure 1.1: Big Bang timeline During the early years of the formation of our galaxy, around 14 billion years ago, there was a transient phase when the first stars were born. To broadly understand the galactic events which occurred during that time, astrophysicists look at the time scale of collision and absorption between different particles. One such event is called the Wouthuysen-Field coupling and is determined by the kinetics of photons undergoing scattering inside a spherical halo of first generation stars during the formation of a galaxy. An accurate computation of the time scale of this process is essential in order to determine certain information about this coupling which relates spin and kinetic temperature of the atoms at this particular epoch in time. The coupling process is related to an emission and absorption process of paticles at a 21 cm wavelength. It is of interest since the lifetimes of the first stars are short and the ionized and heated halos around the first luminous objects are strongly time-dependent. Thus it is theorized and shown that the time scale of the evolution of the expected 21 cm emission/absorption regions can be as small as 105 years. Also, the intensity of the 3 21-cm emission line depends on the density of the neutral atomic hydrogen along the line of sight. Atomic hydrogen along the line of sight contributes to the energy received. This again determines the distance to each clump of hydrogen gas that is detected, enabling physicists to look at a three-dimensional picture of the galaxy. Figure 1.2: Halo of the first star In our work, we also refer to the Λ-CDM or the Lambda Cold Dark Matter model, a simple model which explains certain theories of the Big Bang cosmology. In general most observed physical phenomena are validated by this model. It takes into account the accelerated expansion of the universe and explains the existence of large scale structures of galaxy clusters and other distributed light elements. To illustrate some of these phenomena, we consider a scattering problem inside the halo of a first star. The collisions between the photons are characterized by absorption and redistribution mechanisms that reach their equilibrium to give us the time scale of the event in question. This problem is modeled by a coupled system of Boltzmann type integral differential equations in a non-homogeneous medium and the solution profile involves a sharp front propagating at a given speed. In problems such as this a Weighted Essentially Non-oscillatory (WENO) algorithm is usefull is resolving the solution with high accuracy and limiting oscillations at the same time. However, WENO algorithm was first designed to find solutions of conservation laws. Hyperbolic conservation laws are a class of partial differential equations (PDE)s 4 arising naturally in the physical world in several fields ranging from fluid dynamics to astrophysics. The problems are characterized by the property of appearance of shock waves, across which the solution of the PDE becomes discontinuous. In these areas, a wide range of scientific and engineering problems require efficient computation and high resolution. Different classes of numerical algorithms have been designed such that the numerical solution converges to the unique entropy solution of the conservation law efficiently. These techniques comprise of essentially non-oscillatory methods in the finite difference or finite volume framework, discontinuous Galerkin methods and other finite element schemes and spectral methods. In this work we focus on a numerical algorithm which is based on the first kind: a type of higher order method called Essentially Non-oscillatory (ENO) method. The accuracy criterion in this algorithm is enforced such that when applied to piecewise smooth function, the algorithm gives higher order accuracy in smooth regions and limits oscillations at discontinuities. In our problem, a weighted ENO algorithm coupled with a high order total variation diminishing Runge Kutta method for time evolution is designed to solve a one dimensional scattering equation, with the results published in [27]. It is an extremely stable and robust algorithm to solve such a problem involving propagation of a sharp front with highest resolution. The computation of the time scale of the coupling of the kinetic and spin temperature of the photons required extensive CPU time. Certain advanced techniques, such as an O(N ) algorithm for grouping of terms in the summation at every step of integration of the scattering term in the equation reduce overall computational cost from O(N 2 ) without changing mathematically the algorithm and its accuracy. The algorithm and the results of the physical problem have been published in [29]. A more generalized multi dimensional scattering event inside the spherical halo of a star is solved in the paper in [30]. 5 One of the major challenges of the numerical scheme is that the time of simulation involved in a galactic event is at the scale of millions of years and the time step, in order to guarantee stability, is considerably smaller. Moreover with the spatial domains being large and a refined mesh being a main requirement, an adaptive mesh algorithm has been implemented to successfully solve the problem in realistic time. The adaptive mesh model is a part of the recent work which is submitted to [28]. The motivation for designing this algorithm comes from solving large systems in astrophysics, but the technique is suitable for any problem arising in computational fluid dynamics. Chapter Two Radiative transfer of photons in early universe 7 In this chapter we outline the partial differential equation (PDE) governing the ra- diative transfer of photons which undergo resonant scattering. This chapter also describes the initial and boundary conditions used in the problem. We begin by describing the transfer equation defined in frequency space and continue with de- scribing the four dimensional equation modeling the scattering of photons in the halo of a star. 2.1 Radiative transfer with resonant scattering in frequency space Considering a spatially homogeneous and isotropically expanding infinite medium consisting of neutral hydrogen with temperature T , the kinetics of photons in the frequency space is described by the radiative transfer equation with resonant scat- tering ([17], [31]). ∂J(x, t) cH ∂J(x, t) + 2HJ(x, t) − = ∂t Z v T ∂x −kcφ(x)J(x, t) + kc R(x, x′ )J(x′ , t)dx′ + C(t)φ(x) (2.1) ˙ where J is the specific intensity in terms of the photon number, H(t) = R(t)/R(t) is the Hubble parameter, R(t) is the cosmic factor, vT = (2kB T /m)1/2 is the thermal velocity of hydrogen atom, the dimensionless frequency x is related to the frequency ν and the resonant frequency ν0 by x = (ν − ν0 )/∆νD , and ∆νD = ν0 vT /c is the Doppler broadening. The parameter k = χ/∆νD , and the intensity of the resonant absorption χ is given by χ = πe2 n1 f12 /me c, where n1 is the number density of neutral hydrogen HI at ground state, and f12 = 0.416 is the oscillation strength. The cross 8 section at the line center is πe2 σ0 = f12 (∆νD )−1 . (2.2) me c In equation (2.1), C(t) is the source of photons with the frequency distribution φ(x), which is taken as the Voigt function of the frequency profile with the center at ν0 , i.e. ∞ 2 a e−y Z φ(x) = dy (2.3) π 3/2 −∞ (x − y)2 + a2 Here a is the ratio of the natural to the Doppler broadening. We have a = A21 /(8π∆νD ), and A21 = 6.25 × 108 Hz is the Einstein spontaneous emission coefficient. φ(x) is Z normalized with φ(x′ )dx′ = 1. When a → 0, we have pure Doppler broadening as 1 2 φD (x) = √ e−x . (2.4) π The redistribution function R(x, x′ ) gives the probability of a photon absorbed at the frequency x′ , and re-emitted at the frequency x. It depends on the details of the scattering ([15], [18], [20]). If we consider coherent scattering without recoil, it is R(x, x′ ) = (2.5) ∞      1 xmin + u xmax − u Z −u2 −1 −1 e tan − tan du π 3/2 |x−x′ |/2 a a where xmin = min(x, x′ ) and xmax = max(x, x′ ). In the case of a = 0, i.e. considering only the Doppler broadening, equation (2.5) becomes 1 R(x, x′ ) = erfc[max(|x|, |x′ |)]. (2.6) 2 9 Considering the recoil of atoms, the redistribution function for the Doppler broad- ening is ([12],[2]) 1 ′ 2 R(x, x′ ) = e2bx +b erfc[max(|x + b|, |x′ + b|)], (2.7) 2 where the parameter b = hν0 /mvT c = 2.5 × 10−4 (104 /T )1/2 . It is in the range of 3 × 10−2 - 3 × 10−4 , if the temperature T is in the range of 1 K - 104 K. The redistribution function is normalized as Z R(x, x′ )dx = φ(x′ ). (2.8) Therefore, we have Z Z Z ′ kc φ(x)J (x, t)dx = kc R(x, x′ )J ′ (x′ , t)dxdx′ . (2.9) It means that the total number of photons absorbed given by the term kcφ(x)J(x, t) of equation (2.1) is equal to the total number of scattered photons. Therefore, with equation (2.1), the number of photons is conserved. 2.1.1 Re-scaling the equations We use a new time variable τ defined as τ = cn1 σ0 t, which is in units of the mean free flight time of photons at resonant frequency. The number density of neutral hydrogen atoms n1 = fHI nH , with fHI being the fraction of neutral hydrogen. For the concordance ΛCDM model (briefly described in the introductory chapter), we have nH = 1.88 × 10−7 (Ωb h2 /0.02)(1 + z)3 cm−3 . The factor 0.75 is from hydrogen abundance. Therefore, 10  1/2  3   −1 T 10 0.022 t= 0.054τ fHI yrs (2.10) 104 1+z Ωb h2 We re-scale the equation (2.1) by the following new variables J ′ (x, t) = [R(t)/R0 ]2 J, C ′ (t) = [R(t)/R0 ]2 C. (2.11) Thus, equation (2.1) becomes ∂J ′ (x, τ ) = −φ(x)J ′ (x, τ ) ∂τ ∂J ′ Z + R(x, x′ )J ′ (x′ , τ )dx′ + γ + C ′ (τ )φ(x). (2.12) ∂x We will use J for J ′ and C for C ′ in the equations below. The parameter γ is the so-called Sobolev parameter, γ = (H/vT k) = (8πH/3A12 λ3 n1 ) = (Hme ν0 /πe2 n1 f12 ) (2.13) where λ is the wavelength for Lyα transition. γ is simply related to the Gunn- Peterson optical depth τGP by 1/2  3/2 Ωb h2   −1 5 −1 0.25 1+z γ = τGP = 4.9 × 10 h fHI . (2.14) ΩM 0.022 10 The redshift evolution of fHI is dependent on the reionization models. Before reion- ization fHI ≃ 1; after reionization fHI ≃ 10−5 in average. Therefore, the parameter γ has to be in the range from 1 (z ≤ 7) to 10−7 (z ≥ 10). For a static background, γ = 0. Most of the previous work however has been done on the Fokker-Planck approximation of the equation (2.12), ∂J ′ (x, t)   1 ∂ Z φ(x) = −φ(x)J (x, t) + R(x, x′ )J ′ (x′ , t)dx′ ′ (2.15) 2 ∂x ∂x 11 Which makes equation (2.12),   ∂J(x, t) 1 ∂ ∂J(x, t) ∂J = φ(x) +γ + S ′ (x, t) (2.16) ∂t 2 ∂x ∂x ∂x The initial condition can be taken as zero, ie. when the initial radiative field J(x, τ ) = 0 or the initial radiative field can be 2 J(x, 0) = π −1/2 e−x . (2.17) 2.2 Radiative transfer equation in the 21 cm re- gion The mathematical equations in the last section models only the case of spatial homo- geneity and isotropy. It is equivalent to assume that the Lyα photons are uniformly distributed in the entire 21 cm signal region. This assumption is not trivial, since all the Lyα photons in the 21 cm regions come from first stars, either from direct emission of Lyα photons, or from the Hubble redshifted photons. On the other hand, the 21 cm regions are highly opaque for Lyα. The W-F coupling may not be uni- formly available in the 21 cm signal region, if the time scale of the transfer of Lyα in the physical space is longer than that of the evolution of the 21 cm region. The radiative transfer of Lyα photons in the spherical halo is described by the equation of the specific intensity I(η, r, x, µ) as ∂I ∂I (1 − µ2 ) ∂I ∂I +µ + −γ = ∂η ∂r r ∂µ ∂x Z −φ(x; a)I + R(x, x′ ; a)I(η, r, x′ , µ′ )dx′ dµ′ + S, (2.18) 12 Here we use dimensionless time η defined as η = cnHI σ0 (∆νD )−1 t and dimensionless coordinate r defined as r = nHI σ0 (∆νD )−1 rp , with rp being the physical radial coor- dinate. That is, η and r are, respectively, in the units of mean free flight-time and mean free path of photon ν0 . For a signal propagated in the radial direction with the speed of light, we have r = η + const. In equation (2.18) µ = cos θ is the direction relative to the radial vector r. Similar to the transfer equation in the frequency space, for this equation the re-distribution function R(x, x′ ; a) gives the probability of a photon absorbed at the frequency x′ , and re-emitted at the frequency x. If we consider coherent scattering without recoil, the re-distribution function with the Voigt profile is given by equation (2.5). The re-distribution function of equation (2.5) is normalized as Z ∞ 2 R(x, x′ )dx′ = φ(x, 0) = π −1/2 e−x (2.19) −∞ With this normalization, the total number of photons is conserved in the evolution described by equation (2.18). That is, the destruction process of Lyα photons, such as the two-photon process ([35], [25]), is ignored in equation (2.18). If we take into account the recoil of the atom then the re-distribution can be written as equation (2.5). In equation (2.18), the term with the parameter γ is due to the expansion of the universe. If nH is equal to the mean of the number density of cosmic hydrogen, we −1 have γ = τGP , and τGP is the Gunn-Peterson optical depth. Since the Gunn-Peterson optical depth is of the order of 106 at high redshift ([30]), we take γ = 10−5 − 10−6 in most of our calculations although higher values are also used for different test cases. 13 2.2.1 Eddington approximation When the optical depth is large, we can take the Eddington approximation as I(η, r, x, µ) ≃ J(η, r, x) + 3µF (η, r, x) (2.20) where +1 1 Z J(η, r, x) = I(η, r, x, µ)dµ 2 −1 is the angularly averaged specific intensity and +1 1 Z F (η, r, x) = µI(η, r, x, µ)dµ 2 −1 is the flux. Now we define two new variable as, j = r2 J f = r2 F Equation (2.18) yields the set of equations of j and f as ∂j ∂f ∂j Z + = −φ(x; a)j + R(x, x′ ; a)jdx′ + γ + r2 S, (2.21) ∂η ∂r ∂x ∂f 1 ∂j 2 j + − = −φ(x; a)f. (2.22) ∂η 3 ∂r 3 r The mean intensity j(η, r, x) describes the x photons trapped in the halo by the resonant scattering, while the flux f (η, r, x) describes the photons in transit. 14 Without resonant scattering, equations (2.21) and (2.22) give ∂j ∂f ∂j + = −φ(x; a)j + γ + r2 S, (2.23) ∂η ∂r ∂x ∂f 1 ∂j 2 j + − = −φ(x; a)f. (2.24) ∂η 3 ∂r 3 r For photons with frequency shifts away from ν0 , the halo would be optically thin, and therefore, the Eddington approximation will no longer be valid when the x is large. However we are mainly interested in the profile around x = 0, hence the Eddington approximation is available. The source term S in the equations (2.21) and (2.22) can be described by a boundary condition of j and f at r = r0 . We can take r0 = 0, as r0 is not important if the optical depth of the halo is large. Thus, we have j(η, 0, x) = 0, f (η, 0, x) = S0 φs (x), (2.25) where S0 , and φs (x) are, respectively, the intensity and normalized frequency profile of the sources. Since equation (2.21)-(2.24) are linear, the intensity S0 can be taken as any constant. That is, the solution f (x) of S0 = S is equal to Sf1 (x), where f1 (x) is the solution of S = 1. On the other hand, equations (2.21) and (2.23) are not linear with respect to the function φs (x), i.e. the solution f (x) of φs (x) is not equal to φs (x)f1 (x), where f1 is the solution of φs (x) = 1. In the range outside the halo, r > R, no photons propagate in the direction µ < 0. Therefore, the boundary condition at r = R given by Z −1 µJ(η, R, x, µ)dµ = 0 0 15 is then ([36]) j(η, R, x) = 2f (η, R, x). (2.26) There is no photon in the field before t = 0. Therefore, the initial condition is j(0, r, x) = f (0, r, x) = 0. (2.27) We solve equations (2.20) and (2.21) or equations (2.23) and (2.24) with time- independent sources S and the results are discussed in chapters 6 and 7. Chapter Three The high order accurate algorithm 17 In this chapter we describe the formulation and computation of the high order ac- curate algorithm to solve the Boltzmann type partial differential equation described in the previous section. 3.1 High order WENO based scheme for the RT equation in frequency space This section describes the WENO based higher order method for the solution of equation (2.12) and the numerical treatment of the time evolution. 3.1.1 Computational domain and computational mesh The computational domain in the case of static background is x ∈ [−6, 6]. The initial condition is J(x, 0) = 0. The boundary condition is J(x, τ ) = 0, at |x| = 6. (3.1) In the case of the expanding background, i.e.γ 6= 0, the computational domain is bigger than x ∈ [−6, 6] depending on the value of the Sobolev parameter γ. The domain (xlef t , xright ) is chosen such that for the particular value of γ we have J(xlef t , τ ) ≈ 0, J(xright , τ ) ≈ 0. (3.2) For example, the domain is taken to be x ∈ [−100, 6] for the case of γ = 10−3 . Also for different values of γ the solutions reach saturation at different time. For 18 example, we see in our numerical results that the solution J(ν, t) of equation (2.12) reaches saturation at time τ = 100 for γ = 10−1 , whereas it reaches saturation at time τ = 104 for γ = 10−3 . The computational domain is discretized into a uniform mesh in the x direction, xi = xlef t + i∆x, i = 0, 1, 2, · · · , Nx , where ∆x = (xright − xlef t )/Nx , is the mesh size. We also denote Jin = J(xi , τ n ), the approximate solution values at (xi , τ n ). 3.1.2 Algorithm of the spatial derivative ∂J To calculate ∂x , we use the fifth order WENO method ([21]). That is, ∂J(xi , τ n ) 1 ˆ ˆ j−1/2 ) ≈ (hj+1/2 − h (3.3) ∂x ∆x ˆ j+1/2 is obtained by the procedure given below. We use where the numerical flux h the upwind flux in the fifth order WENO approximation because the wind direction is fixed (negative). First we denote hi = J(xi , τ n ), i = −2, −1, · · · , Nx + 3 (3.4) where n is fixed. The numerical flux from the WENO procedure is obtained by ˆ i+1/2 = ω1 h h ˆ (1) + ω2 h ˆ (2) + ω3 h ˆ (3) , (3.5) i+1/2 i+1/2 i+1/2 19 ˆ (m) are the three third order fluxes on three different stencils given by where h i+1/2 ˆ (1) 1 5 1 h i+1/2 = − hi−1 + hi + hi+1 , 6 6 3 ˆ (2) 1 5 1 h i+1/2 = hi + hi+1 − hi+2 , 3 6 6 ˆ (3) 11 7 1 h i+1/2 = hi+1 − hi+2 + hi+3 , 6 6 3 and the nonlinear weights ωm are given by, ω ˇm γl ωm = , ω ˇl = , (3.6) 3 X (ǫ + βl )2 ω ˇl l=1 where ǫ is a parameter to avoid the denominator to become zero and is taken as ǫ = 10−8 . The linear weights γl are given by 3 3 1 γ1 = , γ2 = , γ3 = , (3.7) 10 5 10 and the smoothness indicators βl are given by 13 1 β1 = (hi−1 − 2hi + hi+1 )2 + (hi−1 − 4hi + 3hi+1 )2 , 12 4 13 1 β2 = (hi − 2hi+1 + hi+2 )2 + (hi − hi+2 )2 , 12 4 13 1 β3 = (hi+1 − 2hi+2 + hi+3 ) + (3hi+1 − 4hi+2 + hi+3 )2 . 2 12 4 20 3.1.3 Time evolution using TVD Runge-Kutta To evolve in time, we use the third-order TVD Runge-Kutta time discretization ([34]). For systems of ODEs uτ = L(u), the third order Runge-Kutta method is u(1) = un + ∆τ L(un , τ n ), 3 n 1 (1) u(2) = u + (u + ∆τ L(u(1) , τ n + ∆τ )), 4 4 1 n 2 (2) 1 u(3) = u + (u + ∆τ L(u(2) , τ n + ∆τ )). 3 3 2 TVD property is the Total Variation Diminishing property of a numerical solution u which says that X T V (u) = |uj+1 − uj | (3.8) j does not increase ie., T V (un+1 ) ≤ T V (un ) (3.9) This is one of the highly applied class of high order TVD time discretization methods which is useful for solving conservation laws with stable spatial discretizations. 3.2 Numerical integration: Higher order scheme and an O(N ) algorithm In this section we describe the algorithm to numerically integrate the scattering term as well as the source terms in equations (2.3) and (2.7) etc. We need a high order accurate numerical integration technique. We also device a grouping algorithm which reduces the number of computations of the integration routine. Since our time 21 evolution of the solution profile is long, in terms of millions of years, and for stability the time step is usually small (in the order of a few years) the numerical algorithm is computationally intensive. To be able to compute the solution in a realistic time we incorporate the following grouping method. 3.2.1 High order numerical integration The integration of the resonance scattering term is calculated by a fifth order quadra- ture formula ([32]) Z xright Nx X f (x)dx = ∆x wj f (xj ) + O(∆x5 ), (3.10) xlef t j=0 where the weights are defined as, 251 299 211 739 w0 = , w1 = , w2 = , w3 = , 720 240 240 720 739 211 299 251 wNx −3 = , wNx −2 = , wNx −1 = , wN x = , 720 240 240 720 and wj = 1 otherwise. 3.2.2 O(N ) algorithm for numerical integration Z We need to numerically integrate R(x, x′ )J(x′ , t)dx′ , denoted as 1 Z ′ 2 I(x) = e2bx +b erfc(max(|x + b|, |x′ + b|))J(x′ , t)dx′ , (3.11) 2 22 with R(x, x′ ) as in equation (2.7). To evaluate Im = I(xm ), ∀m = −Nl , · · · , Nr , we apply the rectangular rule, which is spectrally accurate for smooth functions vanishing at boundaries, xright 1 Z ′ 2 Im = erfc(max(|xm + b|, |x′ + b|))e2bx +b J(x′ , t)dx′ (3.12) 2 xlef t Nr 1 X 2 ≈ ∆x erfc(max(|xm + b|, |xi + b|))e2bxi +b J(xi , t). (3.13) 2 i=−N l Notice that this numerical integration is very costly and uses most part of the CPU time. For the equation with recoil and the redistribution function with a = 0, we have used a grouping of the numerical integration operations at different x locations so that the computational cost can be reduced to order O(N ) rather than O(N 2 ), where N is the number of grid points in x, without changing mathematically the algorithm and its accuracy. The O(N ) algorithm for numerical integration is given below which further highlights the speed and accuracy of the numerical algorithm proposed in this paper. The proposed scheme with order O(N ) computational effort is the following. Let b 2b Nb = f loor( ∆x ) and N2b = f loor( ∆x ). The integration algorithm is designed for two cases: m ≥ −Nb and m < −Nb . 23 In the case of m ≥ −Nb or equivalently xm + b ≥ 0: −m−N X2b −1 1 2 Im = ∆x( erfc(|xi + b|)e2bxi +b J(xi , t) 2 i=−Nl m X 2 +erfc(|xm + b|) e2bxi +b J(xi , t) i=−m−N2b Nr X 2 + erfc(|xi + b|)e2bxi +b J(xi , t)) (3.14) i=m+1 . 1 = ∆x(I1,m + erfc(|xm + b|)I2,m + I3,m ) (3.15) 2 1. Evaluate I1,−Nb , I2,−Nb and I3,−Nb respectively as Nb −N X2b −1 2 I1,−Nb = erfc(|xi + b|)e2bxi +b J(xi , t), (3.16) i=−Nl −Nb X 2 I2,−Nb = e2bxi +b J(xi , t), (3.17) i=Nb −N2b Nr X 2 I3,−Nb = erfc(|xi + b|)e2bxi +b J(xi , t), (3.18) i=−Nb +1 which leads to O(N ) cost. 2. DO m = −Nb + 1, Nr Evaluate I2,m , I3,m respectively by 2 I1,m = I1,m−1 − erfc(|x−m−N2b + b|)e2bx−m−N2b +b J(x−m−N2b , t) (3.19) 2 2 I2,m = I2,m−1 + e2bxm +b J(xm , t) + e2bx−m−N2b +b J(x−m−N2b , t) (3.20) 2 I3,m = I3,m−1 − erfc(|xm + b|)e2bxm +b J(xm , t) (3.21) ENDDO To be consistent with the indices, if Nl − N2b < Nr then, we will set I1,m = 0, for m = Nl − N2b , Nr . The algorithm leads to O(1) cost per m, therefore O(N ) 24 computation overall. In the case of m < −Nb , or equivalently xm + b < 0: m−1 1 X 2 Im = ∆x( erfc(|xi + b|)e2bxi +b J(xi , t) 2 i=−N l −m−N X2b −1 2 +erfc(|xm + b|) e2bxi +b J(xi , t) i=m Nr X 2 + erfc(|xi + b|)e2bxi +b J(xi , t)) (3.22) i=−m−N2b 1 = ∆x(I1,m + erfc(|xm + b|)I2,m + I3,m ) (3.23) 2 1. Evaluate I1,−Nb −1 , I2,−Nb −1 and I3,−Nb −1 as −N X b −2 2 I1,−Nb −1 = erfc(|xi + b|)e2bxi +b J(xi , t), (3.24) i=−Nl NbX −N2b 2 I2,−Nb −1 = e2bxi +b J(xi , t), (3.25) i=−Nb −1 Nr X 2 I3,−Nb −1 = erfc(|xi + b|)e2bxi +b J(xi , t), (3.26) i=Nb +1−N2b which leads to O(N ) cost. 2. DO m = −Nb − 2, −Nl Evaluate I1,m , I2,m , I3,m respectively by 2 I1,m = I1,m+1 − erfc(|xm + b|)e2bxm +b J(xm , t) (3.27) 2 2 I2,m = I2,m+1 + e2bxm +b J(xm , t) + e2bx−m−N2b −1 +b J(x−m−N2b −1 , t) (3.28) 2 I3,m = I3,m+1 − erfc(|x−m−N2b −1 + b|)e2bx−m−N2b −1 +b J(x−m−N2b −1 , t) (3.29) 25 ENDDO To be consistent with the indices, if Nl − N2b > Nr , we will set I3,m = 0, for m = −N2b − Nr − 1, −Nl . Again, the algorithm leads to O(1) cost per m, therefore O(N ) computation overall. Unfortunately, this technique does not work for the case a 6= 0, hence the CPU cost for the case with a 6= 0 is much larger. The function in the Voigt function, in equation (2.3) is such that no grouping of terms can be done. In this case the computation is significantly more costly. 3.3 Algorithm for the RT equation in a spherical halo In this section we outline the numerical method for solving the radiative transfer (RT) equation in the three dimensional spherical halo of the first star (2.18). After the application of the Eddington approximation (2.20) we get the the set of scalar equations (2.21) and (2.22). The dimensionless time variable here is represented by η. 3.3.1 The WENO algorithm: approximations to the spatial derivatives The spatial derivative terms in equations (2.21) and (2.22) are approximated by a fifth order finite difference WENO scheme. 26 ∂j We first give the WENO reconstruction procedure in approximating ∂x , ∂j(η n , ri , xj ) 1 ˆ ˆ j−1/2 ), ≈ (hj+1/2 − h (3.30) ∂x ∆x ˆ j+1/2 is obtained by the fifth with fixed η = η n and r = ri . The numerical flux h order WENO approximation in an upwind fashion, because the wind direction is fixed (negative). Denote hj = j(η n , ri , xj ), j = −2, −1, · · · , N + 3 (3.31) with fixed n and i. The numerical flux from the WENO procedure is obtained by a similar procedure as given in the previous section. To approximate the r- derivatives in the system of equations (2.21) and (2.22), we need to perform the WENO procedure based on a characteristic decomposition. We write the left hand side of equations (2.21) and (2.22) as ut + Aur (3.32) where u = (j, f )T and    0 1  A=  1 3 0 is a constant matrix. To perform the characteristic decomposition, we first compute the eigenvalues, the right eigenvectors, and the left eigenvectors of A and denote them by, Λ, R and R−1 . We then project u to the local characteristic fields v with v = R−1 u. Now ut + Aur of the original system is decoupled as two independent equations as vt + Λvr . We approximate the derivative vr component by component, each with the correct upwind direction, with the WENO reconstruction procedure ∂j similar to the procedure described above for ∂x . In the end, we transform vr back 27 to the physical space by ur = Rvr . We refer the readers to Cockburn et al. 1998 [[9]] for more implementation details. 3.3.2 Implementation of the boundary condition The source term given in equation (2.21) is implemented as a boundary condition on f (η, r = r0 , x). f (η, r = r0 , x) = s0 φs (x) For the intensity j, a reflective boundary condition is used at r = r0 . At the boundary of r = rmax , x = xlef t and x = xright , we use zero boundary conditions for both j and f , because of the way we choose rmax , xlef t and xright . 3.3.3 Adaptive mesh procedure for non-uniform grid A fifth order conservative finite difference WENO scheme can only be applied to a uniform grid or a smoothly varying grid [33]. A smooth transformation, ξ = ξ(r) (3.33) gives us a uniform grid in a new variable ξ. In this case ξ is sufficiently smooth, i.e., it has as many derivatives as the order of accuracy of the scheme. Therefore the left hand side of the (2.21) and (2.22) as ut + Aur (3.34) 28 where u = (j, f )T and    0 1  A=  1 3 0 is transformed to ut + Aξr uξ (3.35) and the WENO derivative approximation is now applied to uξ . We use a gaussian function for ξ(r) which clusters points tracking the left moving front. The velocity of the moving front is calculated aprroximately by the Sobolev parameter γ. The process of creating the non-uniform mesh is to start with a partic- ular r-domain and number of points in the mesh. Thereafter we increase the domain size and points at particular intervals of time, in our case at η = 1000, 2000 . . .. To calculate the solution at the points in the domian before the increment we use a high order interpolation technique. Our algorithm has been implemented for short time calculation as in under investigation for longer time simulations. 3.4 Conservation and normalization techniques We have seen in the earlier chapter that the number of photons is conserved in the process of scattering. We can use the equation (2.9) to test the accuracy of our algorithm. 29 3.4.1 Test with the conservation of the photon number From equation (2.21) we have ∂ ∂ Z Z jdx + f dx = 0 (3.36) ∂η ∂r R Therefore, for time-independent solution we have f (r, x)dx = const. It yields Z Z f (r, x)dx = f (0, x)dx = S0 (3.37) That is, the flux at saturated states is r-independent. This is the conservation of the number of photons. 3.4.2 Normalization of the re-distribution term R(x, x′ )jdx. The rectangular rule is R We apply the rectangular rule to evaluate known to have spectral accuracy for smooth integrated function R(x, x′ )j(x′ ) with compact support. In equations (2.21) and (2.22), it is known that Z R(x, x′ )dx′ = φ(x), (3.38) which is crucial for photon conservation. However, X R(x, xj )∆x = φ(x) j 30 may not be true in general. To numerically preserve the photon conservation, we first compute φ(x) ratio(x, ∆x) = P , (3.39) j R(x, xj )∆x then we normalize the collision term by approximating Z R(x, x′ )j(t, r, x′ )dx′ (3.40) with X ratio(x, ∆x) R(x, xj )j(η, r, xj )∆x (3.41) j with fixed η and r. We see that this quadrature formula is exact when j(t, r, x) = 1. By using this technique the quadrature formula (3.41) exactly approximates the integration term (3.40). Chapter Four Radiative transfer in the frequency space 32 In this chapter we describe the numerical solution of the integral differential equa- tion (2.1) which describes the radiative transfer of photons in the frequency space with resonant scattering of Lyα photons by hydrogen gas. The scattering event is characterized by 21 cm absorption and emission of photons from the halo of the first generation stars. It is believed that detecting the red shifted 21 cm signals the from early universe is one of the new frontiers in observational cosmology. The study of the scattering phenomenon is directly related to that of Wouthuysen-Field (W-F) coupling ([5],[23],[8]) which is the coupling of the spin temperature of Lyα photons with the kinetic temperature of Hydrogen gas. 21 cm absorption and emission de- pend on several factors including the fraction of neutral Hydrogen and Lyα photons that are available for the W-F coupling. Hence, the goal is to show that the coupling is possible only if the time scale of the onset of the coupling is less that the forma- tion and evolution of 21 emission and absorption shells. Most of the previous work ([7], [16], [13], [8]) on this subject has been done on a time independent solution using a Focker-Planck approximation of equation (2.1). In this chapter we study the time dependent equation (2.1) using a higher order WENO algorithm coupled with a higher order TVD-Runge Kutta method. 4.1 Radiative transfer (RT) equation in the fre- quency space 4.1.1 RT equation in static background WENO schemes have been widely used in solving Boltzmann equations ([4]). In this section we first present the test for the numerical solver against known analytic 33 solutions and then use it to solve for the time scale of the W-F coupling. We first test the WENO solver with two analytical solutions of equation (2.12) with static background, H = γ = 0 (2.13) and Doppler broadening a = 0 ([11]). The first analytical solution is for the initial radiative field J(x, τ ) = 0 and the constant source C(τ ) = 1. It is 2 J(x, τ ) = π −1/2 [1 − exp(−τ e−x )] (4.1) Z ∞ 2 2 2 + ew [1 − (1 + τ e−w ) exp(−τ e−w )]erf(w)dw. x The second analytical solution of equation(2.12) is also for H = γ = 0, a = 0, but the source C = 0, while the initial radiative field is 2 J(ν, 0) = π −1/2 e−ν . (4.2) The solution is Z ∞ −1/2 −x2 −x2 2 2 J(x, τ ) = π e exp(−τ e )+τ e−w exp(−τ e−w )erf(w)dw. (4.3) x The analytical solutions (4.1) and (4.3) are shown in figures 4.1 and 4.2 respectively. We also plot the numerical solutions given by our algorithm in the same figures. The numerical results show very small deviations from the analytical solutions. A common feature of figures 4.1 and 4.2 is that the original Doppler peak at the center of the frequency profile gradually becomes a flat plateau. The width of the plateau increases with time. This is because the resonant scattering makes a non-uniform distribution in the frequency space a uniform one. It is similar to the physical space when diffusion leads to an evolution from non-uniform distribution to a uniform one. The height of the plateau of figure 4.1 is increasing with time, while 34 τ = 10000 3 10 τ = 1000 2 10 τ = 100 101 τ = 10 J(x) 0 10 τ=1 -1 10 τ = 0.1 10-2 τ = 0.01 -3 10 -4 -2 0 2 4 6 x Figure 4.1: Static solutions (γ = 0) for pure Doppler redistribution (a = 0) of equation(2.12), in which C = 1 and J(x, 0) = 0. The analytical solutions are shown by dashed lines, while the numerical results are shown by solid lines. in figure 4.2 it is decreasing, because for the case C = 1, the number of photons increases, while for the case of C = 0, the total number of photons is conserved. 4.1.2 Expanding background In this section we look at the profile of the radiation field in an expanding background ie, when γ 6= 0 We also perform a test to show the stability and accuracy of our algorithm by considering expanding background. Without absorption and scattering, equation (2.12) becomes ∂J ∂J =γ + C(τ )φ(x). (4.4) ∂τ ∂x If C is τ -independent, the analytical solution of equation (4.4) is given in ([31]) as, J(x, τ ) = J(x + γτ, 0) + Cγ −1 [Φ(x) − Φ(x + γτ )] (4.5) 35 τ = 0.01,0.1 τ=1 τ = 10 τ = 100 τ = 1000 τ = 10000 10-1 J(x) -2 10 -3 10 -4 -2 0 2 4 x Figure 4.2: Static solutions (γ = 0) for pure Doppler redistribution (a = 0) of equation(2.12), in which C = 0 and J(x, 0) = π −1/2 exp(−x2 ). The analytical solutions are shown by dashed line, while the numerical results are shown by solid lines. where Z ∞ Φ(x) = φ(x′ )dx′ . (4.6) x The physical interpretation of the solution of equation (4.5) is that the redshift of photons in frequency space is described by −x = γτ. (4.7) We take √ 2 φ(x) = (1/ π)e−x , (4.8) and assume that the initial field to be J(x, 0) = 0. Figure 4.3 shows that the numerical results obtained by using the WENO algorithm agree well with the an- alytical solution of equation(4.5). Figure 4.4 gives the solutions of equation (2.12) with parameter γ = 10−1 and 10−3 . In these cases the analytical solutions are not available. We can see from figures 4.3 and 4.4 that the intensity J stops increasing and approaches a saturated value when τ is large. This happens when the number 36 3 10 101 τ = 10000 τ = 1000 τ = 100000 τ = 1000 τ = 10 τ = 100 102 τ = 100 1 10 0 10 0 J(x) 10 J(x) 10-1 -2 10 -1 10 -3 10 10-4 -100 -80 -60 -40 -20 0 -100 -50 0 x x Figure 4.3: Solutions of J(x, τ ) of equation(4.4) with initial condition J(x, 0) = 0. The parameter γ is taken to be 10−1 (left panel) and 10−3 (right). The analytical solutions equation(4.5) are shown by dashed lines . 103 101 τ = 1000 τ = 40000 τ = 20000 τ = 10000 τ = 1000 102 τ = 100 τ = 10 τ = 100 101 τ = 10 0 10 100 J(x) J(x) -1 10 -2 -1 10 10 -3 10 -4 10 -100 -80 -60 -40 -20 0 -40 -30 -20 -10 0 x x Figure 4.4: Solutions of J(x, τ ) of equation (2.12) with initial condition J(x, 0) = 0. The param- eter γ is taken to be 10−1 (left) and 10−3 (right). of photons produced from the sources is equal to that of the red shifted photons. 4.1.3 The effect of recoil of the atom The final test is achieved by using the redistribution function with recoil equation (2.7). Figure 4.5 gives the solutions of equation (2.12) with parameter b = 0.03, and the initial condition and source term are taken to be the same as figure 4.1. The 37 τ = 10000 3 10 τ = 1000 102 τ = 100 1 10 τ = 10 J(x) 0 10 τ=1 10-1 τ = 0.1 -2 τ = 0.01 10 -4 -2 0 2 4 x Figure 4.5: Solutions of J(x, τ ) of equation(2.12) with redistribution function equation (2.5), and γ = 0. results of figure 4.5 are similar to figure 4.1, but the center flat plateau of figure 4.1 now is replaced by a sloping plateau. The latter is a local Boltzmann distribution around the resonant scattering. Previous work in [12] has shown that, once the solution J(x, τ ) of equation (2.1) with b = 0 has a flat plateau in the central region, the effect of recoil is to make the flat plateau to approach Boltzmann-like distribution, i.e. if the solution with no recoil shows J(x, τ ) ≃ J(0, τ ), within |x| < xf , the solution of equation (2.12) with redistribution equation(4.4) will be J(0, τ )e−2bx = J(0, τ )e−h(ν−ν0 )/kT , |x| < xf . (4.9) One can test this property with F (τ ) defined as F (τ ) = log J(0, τ ) − log(1, τ ) (4.10) If the local Boltzmann distribution within |x| < 1 is realized at time τ , F (τ ) should be equal to 2b [equation (4.9)]. In figure 4.6, we present the relation of F (τ ) vs. 38 τ . At η = 0, the Gaussian sources [equation (4.3)] yields F (0) = 1. Then, F (τ ) approaches to 2b at τ > 102 . These results show that the WENO solver is reliable and effective to study the time-dependence of J described by the resonant scattering equation. Figure 4.6 also shows the relation F (τ ) vs. τ for solutions b = 0.03 (left panel) and b = 0 (right panel). The two curves actually are similar. When the right curve approaches to F (τ ) → 0 at time τ , the left curve, at the same time, approaches to F (τ ) → 2b. Therefore, it is reasonable to calculate the time-scale of the formation of local Boltzmann distribution by the formation of the flat plateau. 0.8 0.8 0.6 0.6 F(τ) F(τ) 0.4 0.4 0.2 0.2 -2 -1 0 1 2 3 4 0 -2 -1 0 1 2 3 4 10 10 10 10 10 10 10 10 10 10 10 10 10 10 τ τ Figure 4.6: F (τ ) as function of τ for solutions shown in figure 4.5 (left) and figure 4.1 (right) 4.2 The Wouthuysen-Field Coupling In this section we introduce the idea of the Wouthuysen-Field coupling and using the time evolution of the solution profile of equation (2.12), we answer some questions ababoutt the physical properties associated with the coupling mechanism. We wish to estimate the time scale, τWF , of the onset of the W-F coupling. As 39 mentioned in the earlier section, in the first phase, J(x, τ ) keeps the initial profile φ(x) of the photons from source C. In the second phase, the profile is no longer the initial one. A flat plateau (without recoil) or local Boltzmann distribution (with recoil) form the central part (|x| ≃ 0) due to the resonant scattering. The width of the flat plateau increases with the time τ . Finally, in the third phase, the injection of Lyα photons from the source C is balanced by the redshift, the height of the flat plateau or the local Boltzmann of J(x, τ ) will stop increasing, and will reach a saturated value. The width of the flat or sloping plateau on the red side continuously increases. As mentioned in the previous chapter, most previous works on the effect of the W-F coupling of 21 cm signal are based on the static solution of the Fokker-Planck approximation of equation (2.1) ([7], [16], [13], [8]). They have solved the approxi- mated equation (2.16). Obviously, the solution of the scattering equation corresponds to the third phase. Thus, the time scale of the W-F coupling is of order τIII . There- fore, this approximation would be reasonable if the time scale of the formation and evolution of the 21 cm signal regions is larger than τIII . From figures 4.4 and 4.6, one can see that the time scale of the onset of the third phase, τIII , has to be of the order of 102 , and 104 for γ = 10−1 , and 10−3 respectively, i.e. τIII is roughly equal to ≀N of Gunn-Peterson optical depth. On the other hand, we can see from equation (2.12) that γ = 10−1 , and 10−3 correspond to fHI = 10−4 , and 10−2 for the source at redshift 1 + z = 10. Thus, for all cases, equation (2.9) yields that the time scale tIII is of the order of 1 Myr. The static solution of the Fokker-Planck equation would not be valuable for the 21 cm problem, if the time scale of the evolution of the 21 cm region is equal to or less than 1 Myr. Actually the time-independent solution of equation (2.12) would not be neces- 40 sary for the 21 cm signal. The relative occupation of the two hyperfine-structure components of the hydrogen ground state depends only upon the shape of the spec- trum near the Ly-α frequency, regardless whether the solution is time-independent, or saturated. What we need for the W-F coupling is only the frequency distri- bution to show a local Boltzmann-like distribution J(x, τ ) ∝ exp[−h(ν − ν0 )/kT ] around ν0 , where T is the kinetic temperature of hydrogen gas. The formation of the Boltzmann-like distribution is irrelevant for the time scale τIII , but depends on the onset of the second phase of the J(x, τ ) evolution. Therefore, τWF should be estimated by the time scale of the onset of the second phase, τII , not by τIII . According to the property shown by figure 4.6, τII can be found by the formation of a flat plateau around the resonant frequency. For small γ, τII is much less than τIII . To show this point, we calculate the solution of equation (2.12) with γ = 10−5 and 10−7 on static background. The results are given in figure 4.7. It shows that a small flat plateau has already formed at τ = 100. On the other hand, τIII will be as large as τ ≃ 106 for γ = 10−5 . τ = 10000 τ = 10000 103 103 τ = 1000 τ = 1000 2 2 10 10 τ = 100 τ = 100 1 1 10 10 τ = 10 τ = 10 J(x) J(x) 100 τ=1 100 τ=1 -1 -1 10 τ = 0.1 10 τ = 0.1 -2 -2 10 τ = 0.01 10 τ = 0.01 -3 -3 10 10 -6 -4 -2 0 2 4 6 -6 -4 -2 0 2 4 6 x x Figure 4.7: Solutions of J(x, τ ) of equation (2.12) with γ = 10−5 (left panel) and 10−7 (right panel). In figure 4.2, we show the function F (τ ) for the parameters (a = 0, γ = 10−5 ), and (a = 0, γ = 10−7 ). It is interesting to observe that the two sets of F (τ ) are almost identical. That is, the time scale τWF is actually independent of the parameter 41 0.8 0.8 0.6 0.6 F(τ) F(τ) 0.4 0.4 0.2 0.2 0 -2 -1 0 1 2 3 4 0 -2 -1 0 1 2 3 4 10 10 10 10 10 10 10 10 10 10 10 10 10 10 τ τ Figure 4.8: F (τ ) as function of τ for solutions of a = 0 and γ = 10−5 (left) and 10−7 (right). γ. For both cases of the parameter γ, τWF is equal to the ≀N of 102 . In other words, a flat plateau with proper width can form with a ON of 102 of the mean free flight time of the resonant photons, regardless of the expansion of the universe. The reason for the γ-independence of τII or τWF can be directly seen from equation (2.12). This equation describes two kinetic processes of approaching a statistically steady state. The first is the resonant scattering, which leads to the formation of a flat plateau or a local Boltzmann distribution around the resonant frequency. The second is the redshift of photons, which leads to a statistical equilibrium between the injected photons and redshifted photons. The steady state due to resonant scattering corresponds to the second phase, while the steady state due to redshift corresponds to the third phase. Similar to various statistical equilibrium or steady states maintained by collision or scattering, the steady state of resonant scattering can be realized via a few tens or hundreds of the resonant scattering. On the other hand, γ is a ratio between the time scales of the resonant scattering and the expansion of the universe. Therefore, when γ is much less than 1, the time scale of the formation of the flat plateau and local Boltzmann distribution are much shorter than the expansion of the universe. This point can also be seen with figure 4.6, which shows that for γ = 10−5 and 10−7 the evolutions of J(x, τ ) actually are the same until τ = 104 . In short, 42 the W-F coupling will take place once a statistical equilibrium localized around the resonant frequency is realized. 4.3 Discussion In this chapter we solve the radiative transfer equation in one space dimension and also introduce the W-F coupling and show the γ-independence of the coupling pro- cess. In the next chapter we answer the questions of the dependence of this physical process on other parameters such as recoil of atoms. Chapter Five Time scale of Wouthuysen Field coupling and intensity calculation 44 In this chapter we present the study of the Wouthuysen-Field coupling in the early universe with numerical solutions of the integro-differential equation describing the kinetics of photons undergoing resonant scattering. We focus on the time evolution of the Wouthuysen-Field (W-F) coupling in relation to the 21 cm emission and absorption at the epoch of reionization. We are looking at the equation (2.12) with the same redistribution and source functions as in the previous chapter. In order to calculate the width of the Boltzmann distribution we explain a few concepts here. The physical meaning of the terms on the right hand side of equation (2.12) is clear. The first term is the absorption at frequency x, the second term is the re- emission of photons with frequency x by scattering, and the third term describes the redshift of photons. The time scale of a photon moving ∆x in the frequency space is equal to ∆τ = γ −1 ∆x. (5.1) This actually is due to the Hubble expansion. Considering equation (2.9), equation (2.12) gives d Z Z Jdx = Sdx. (5.2) dτ Z This equation shows that the total number of photons Jdx is dependent only on the sources, regardless of the parameter b of the resonant scattering. Since numerical errors accumulated over a long time evolution could be significant, equation (5.2) is also useful to check the reliability of a numerical code. 45 5.1 Wouthuysen-Field coupling in a static back- ground To study the effect of atomic recoil, we first solve the time evolution of J(x, τ ) in a static background, i.e. γ = 0. A typical time-dependent result is shown in figure 5.1, in which b = 0.03. The solution of figure 5.1 is actually similar to that shown in figure 4.1. Figure 4.1 shows that a flat plateau around x = 0 is to be formed at τ > 100, while figure 5.1 shows that J(x, τ ) evolves into a Boltzmann distribution around x = 0 as J(x, τ ) ≃ J(0, τ )e−2bx = J(0, τ )e−h(ν−ν0 )/kT , |x| ≤ xl . (5.3) This is the so-called “Planck shape in a region around the initial frequency” [38]. The expression (5.3) has also been found by Field in [12]. We will call this feature to be a local Boltzmann distribution. The width xl of the local Boltzmann distribution is numerically defined by the frequency range |x| ≤ xl , in which the slope ln J(x, τ )/dx deviating from 2b is small (see below). τ = 10000 3 10 τ = 1000 102 τ = 100 1 10 τ = 10 J(x) 0 10 τ=1 10-1 τ = 0.1 -2 τ = 0.01 10 -4 -2 0 2 4 x Figure 5.1: WENO numerical solutions of equation (2.12) with γ = 0 and b = 0.03. The source is taken to be S = φ(x) and initial condition J(x, 0) = 0. 46 Figure 5.2 plots the time-evolution of J(x, τ ) with different parameter b, and a blow-up figure at time τ = 104 . The zoom-in figure shows clearly that the integral R J(x, τ )dx is b-independent. This is consistent with the photon number conservation equation (2.18). It shows again that the WENO algorithm is robust. We can also see from the right panel of figure 5.2 that all the curves of ln J(x, τ ) vs. x at τ = 104 and in the range −2 < x < 2 can be approximated by a straight line. That is, the width xl of the local Boltzmann distribution shown in figure 5.2 is equal to about 5.1, and it is approximately b-independent. The formation and evolution of the local Boltzmann distribution can be quantitatively described by B(τ ) defined as B(τ ) = 2b/[ln J(0, τ ) − ln J(1, τ )], (5.4) where ln J(0, τ )−ln J(1, τ ) is the slope of the straight line ln J(x, τ ) vs. x for |x| ≤ 1. This function is the reciprocal of the function F (τ ) defined in the previous chapter, equation (4.10) except for the constant 2b in the numerator. It is defined in this way to better demonstrate the effect of recoil of the atom on the solution profile. For a Gaussian source, 2 √ S(x) = φ(x) = e−x / π (5.5) we have B(0) = 2b, and B(τ ) approaches 2b at large τ . Figure 5.3 presents the numerical relation of B(τ ) vs. τ . The slopes [log J(0, τ ) − log J(1, τ )] at τ = 105 are, respectively, 0.0601 for b = 0.030, 0.0303 for b = 0.015, 0.0159 for b = 0.0079, and 0.0051 for b = 0.0025. That is, within the frequency range |x| ≤ xl and xl = 1, the relative deviation of the slope d ln J(x, τ )/dx from 2b is no larger than 2%. Thus, τ = 105 can be considered as the time scale of forming a local Boltzmann distribution within |x| < xl = 1. For a small width xl < 1, this time scale is as small as ≃ 104 . Therefore, the time scale of the onset of W-F coupling with the width xl equal to about Doppler broadening is 104 -105 . 47 τ = 10000 3 10 τ = 1000 102 τ = 100 1 10 τ = 10 J(x) 0 10 τ=1 10-1 τ = 0.1 -2 τ = 0.01 10 -4 -2 0 2 4 x 2400 2200 2000 τ = 10000 1800 1600 J(x) 1400 1200 1000 800 -2 -1 0 1 2 x Figure 5.2: WENO numerical solutions of equation (2.12) with γ = 0 and b = 0.03 (solid), 0.015 (dashed), 0.0079 (dot-dot-dashed) and 0.0025 (dot-dashed). The source is taken to be S = φ(x) and initial condition J(x, 0) = 0. The bottom panel is a blow-up of the dashed square of the top panel. We can relate the width xl to the mean number of scattering, Nc , needed to form the local Boltzmann distribution. Although the redistribution function equation (2.5) is b-dependent, the probability of x photons undergoing a resonant scattering per unit time is φ(x), which is b-independent. Thus, at a given time τ , the mean number Nc of resonant scattering of photons within |x| ≤ xl is approximately, xl 1 Z Nc ≃ τ φ(x)dx. (5.6) xl 0 48 1.2 1 0.8 0.6 B(τ) 0.4 0.2 0 -0.2 -2 -1 0 1 2 3 4 5 10 10 10 10 10 10 10 10 τ Figure 5.3: B(τ ) vs. τ for solutions in figure 5.2 with b = 0.03 (solid), 0.015 (dashed), 0.0079 (dot-dot-dashed) and 0.0025 (dot-dashed). Equation (5.6) gives the “finite but large number of scattering” for realizing a lo- cal Boltzmann distribution within x ≤ xl (τ ) ([38]). Therefore, the approximate b-independence of xl (figures 5.2 and 5.3) would imply the b-independence of Nc . 5.2 Wouthuysen-Field coupling in an expanding background 5.2.1 Width of the local Boltzmann distribution Considering an expanding background, i.e. γ 6= 0, we solve equation (2.12) by the WENO algorithm. Figure 5.4 plots solutions with the same source S = φ(x) and parameter b = 0.03 as in figure 5.1, but with γ = 10−3 and 10−5 . Similar to figure 5.1, a local Boltzmann distribution has formed when τ ≥ 103 for both γ = 10−3 and 10−5 . The section of the spectrum near x = 0 becomes τ -independent when τ ≥ 104 for γ = 10−3 , and τ ≥ 106 for γ = 10−5 . We refer to this τ -independence as saturation of 49 the profile J(x, τ ) around the resonant frequency. In the saturated state, the number of photons red shifted from ν > ν0 to the local Boltzmann distribution area ≃ ν0 due to Hubble expansion is equal to the number of photons leaving from ν0 to the red wing. Therefore, we see from figure 5.4 that once J reaches the saturation state, the boundary on the red wing (x < 0) of J is moving left (red). On the other hand, the boundary on the blue wing (x > 0) is almost time-independent. Unlike in figure 5.1, the width xl does not always increase with time. For γ = 10−3 the width ceases to increase when τ > 103 , and for γ = 10−5 , it is stopped at τ > 105 . One can find the mean scattering number Nc in the similar way as equation (2.21). When γ 6= 0, the time duration of photons staying in the frequency space from x to x + ∆x is roughly ∆x/γ [equation (5.1)]. On the other hand, the mean probability 1 xl Z of |x| ≤ xl photons being scattered in a unit τ is φ(x)dx. The larger the xl , xl 0 the less the probability. Thus, all photons within |x| < xl on averagely undergo Nc scattering given by xl 1 Z Nc ≃ φ(x)dx. (5.7) γ 0 From figure 5.4, the maximum width for γ = 10−3 is estimated as xl = 1.9, corre- sponding to Nc ≃ 0.5 × 103 . While for γ = 10−5 , maximum width is xl = 2.8, and Nc ≃ 5.0 × 105 . Once the width xl stops to increase, all quantities in equation (2.22), γ, xl and φ(x), are τ -independent. Thus, Nc should also be τ -independent. The τ -independence of the width xl is also shown in figure 5.5, in which we still use γ = 10−3 , and J(x, 0) = 0 initially. However, the source is taken to be S(x) = φ(x − 10) (5.8) That is, the source photons have frequency x = 10, or ν = ν0 +10∆νD . The resonant 50 1500 3 10 τ = 10000 1000 τ = 10000 τ = 1000 102 J(x) 500 J(x) 1 10 τ = 1000 0 10 -10 -5 0 5 -4 -3 -2 -1 0 1 2 3 4 x x 5 τ = 1000,000 10 τ = 100,000 τ = 1000,000 5 10 3 τ = 10000 10 τ = 1000 1 τ = 100 10 J(x) J(x) 10-1 τ = 100,000 4 10 -3 10 -15 -10 -5 0 5 -4 -2 0 2 x x Figure 5.4: WENO numerical solutions of equation (2.12) with b = 0.03 and γ = 10−3 (two top panels) and γ = 10−5 (two bottom panels). The source is taken to be S = φ(x) and initial condition J(x, 0) = 0. The right panels are blow-up of the dashed square of the corresponded left panel scattering at x = 0 (ν = ν0 ) will occur when these photons have red shifted from ν to ν = ν0 , which takes a time of about τ = 10/γ ≃ 104 . Figure 5.5 shows that the whole distribution of J(x, τ ) evolves significantly with time, but the width of the local Boltzmann distribution around x = 0 is maintained at xl ≃ 2 from τ = 104 to 4×104 . We also find from our numerical calculations that when τ > 4 × 104 , the intensity of the photon flux around x = 0 remains constant, or is in a saturated state. Similar to figure 5.2, figure 5.5 also shows that the width of the local Boltzmann distribution is approximately b-independent. From equation (5.7), one can also expect that the width will be smaller for larger γ. A local Boltzmann distribution can form only if 51 3 τ = 40000 3 τ = 40000 3 τ = 40000 10 τ = 10000 10 τ = 10000 10 τ = 10000 τ = 1000 τ = 1000 τ = 1000 102 102 102 J(x) J(x) J(x) 1 1 1 10 10 10 0 0 0 10 10 10 -10 -5 0 5 10 -10 -5 0 5 10 -10 -5 0 5 10 x x x 1500 1500 1500 τ = 40000 τ = 40000 τ = 40000 1000 1000 1000 τ = 10000 τ = 10000 τ = 10000 J(x) J(x) J(x) 500 500 500 -4 -3 -2 -1 0 1 2 3 4 -4 -3 -2 -1 0 1 2 3 4 -4 -3 -2 -1 0 1 2 3 4 x x x Figure 5.5: WENO numerical solutions of equation (2.12) with γ = 10−3 and b = 0.0079 (left), 0.015 (middle), 0.03 (right). The source is taken to be S = φ(x − 10) and initial condition J(x, 0) = 0. The bottom panels are the blow-up of the dashed square of the top panels. γ −1 is large enough. This property is shown with figure 5.6, in which we use the same photon source S(x) and parameter b as in figure 5.5, but we take larger γ. Figure 5.6 presents the results of γ = 1 and 10−1 . We see from figure 5.6 that in the case of γ = 1, there is no local Boltzmann distribution at any time. The resonant scattering leads only to a valley around x = 0. It is because the strongest scattering is at x ≃ 0, which moves photons with frequency x ≃ 0 to x 6= 0. However, the redshift allows the photon quickly leave from x ≃ 0. They do not undergo sufficient number of scatterings to form a local statistical equilibrium distribution. For γ = 10−1 , it seems to show a local Boltzmann distribution, but its width is very small for all time. 52 τ = 100 τ = 1000 1 10 τ = 10 8 τ = 100 0.95 6 4 J(x) J(x) 0.9 0.85 2 -5 0 5 10 -5 0 5 10 x x τ = 100 τ = 1000 1 10 8 0.99 6 J(x) J(x) τ = 10 4 0.98 -2 -1 0 1 2 -2 -1 0 1 2 x x Figure 5.6: WENO numerical solutions of equation (2.12) with γ = 1 (left) and γ = 10−1 (right). Parameter b = 0.0079 (dot-dot-dashed), 0.015 (dashed), 0.03 (solid). The source is S = φ(x − 10) and initial condition J(x, 0) = 0. The bottom panels are the blow-up of the dashed square of the top panels. 5.2.2 Photon source and W-F coupling The τ - and b-independence of the shape and width of the local Boltzmann distribu- tion yield the important conclusion that for given parameters γ and b, the forma- tion and evolution of the local Boltzmann distribution is independent of the photon sources S(x, τ ). This is because equation (2.12) is linear in J. Any source S(x, τ ) can always be considered as a superposition of many monochromatic sources around frequency x = xi , or X S(x, τ ) = Si φ(x − xi , τ ), (5.9) i where Si is the intensity of photon source φ(x − xi , τ ) with frequency xi . J(x, τ ) can P be decomposed into J(x, τ ) = i Ji (x, τ ), where Ji (x, τ ) is the solution of equation (2.12) with the source i. Thus, if the formation of the local Boltzmann distribution 53 around x = 0 is independent of τ and b, the superposition X J(x, τ ) = Ji (x, τ ) (5.10) i should also show the same local Boltzmann distribution around x = 0. Although the overall amplitude does depend on the source, the shape in the neighborhood of the resonant frequency does not. As an example, figure 5.7 presents a solution with the same parameters as in figure 5.5, but the source is with continuous spectrum given by   (10/x)2 , 10 < x < 15,  S(x, τ ) = (5.11)  0  otherwise Photons with frequency x = 10 will arrive earlier at x = 0 with higher intensity, while those with frequency x = 15 will arrive at x = 0 later with lower intensity. The flux J(x, τ ) of figure 5.7 has very different shape from figure 5.5, while the local Boltzmann distribution at −2 < x < 2 of figure 5.7 is exactly the same as that in figure 5.5. τ = 10000 600 3 10 τ = 1000 τ = 10000 500 102 400 1 10 J(x) J(x) 0 10 300 -1 10 -5 0 5 10 15 -3 -2 -1 0 1 2 3 x x Figure 5.7: WENO numerical solutions of equation .(2.12) with γ = 10−3 and b = 0.0079 (dot- dot-dashed), 0.015 (dashed), 0.03 (solid). The source is given by equation (5.11). The right panel is a blow-up of the dashed square of left panel. 54 5.2.3 Intensity From the results in the previous sections, we see that the intensity of the photon flux J at the local Boltzmann distribution is strongly dependent on γ and τ . Figures 5.4 and 5.5 show that at an early time J is smaller than its saturation state. For the solution of figure 5.5, the flux in the frequency range of the local Boltzmann distribution is saturated at about τ = 4 × 104 with a saturated flux J ≃ 104 , while the intensity at τ = 104 is significantly lower than that of the saturated state. Figure 5.8 is the same as figure 5.5, but taking γ = 10−2 and 10−4 . The lower panels of figure 5.8 shows once again that the time-evolution of the intensity is approximately b-independent. Figure 5.8 shows also that for γ = 10−2 , the photon flux approaches the saturated state with intensity of J ≃ 102 at τ = 104 . While for γ = 10−4 , the saturated state has not yet been approached even with an intensity of J ≃ 104 , and τ = 105 . Generally, the smaller the γ, the larger the saturated intensity and the longer the τ needed to approach its saturated state. 5.2.4 Time scales of the frequency shift On the right hand side of equation (2.12), the first term is the absorption at the resonant frequency x, the second term is the re-emission of photons with frequency x by scattering, and the third term describes the Hubble redshift of the photons. We first solve equation (2.12) by dropping the terms of absorption and re- emission. The equation is ∂J ∂J =γ + S(x, τ ). (5.12) ∂τ ∂x 55 τ = 10000 τ = 100000 4 2 10 10 τ = 10000 τ = 1000 103 1 10 2 10 J(x) J(x) 0 101 10 0 10 10-1 10-1 -5 0 5 10 -5 0 5 10 x x 120 7000 τ = 10000 100 6000 τ = 100000 5000 80 4000 τ = 1000 60 J(x) J(x) 3000 40 2000 -3 -2 -1 0 1 2 3 -4 -3 -2 -1 0 1 2 x x Figure 5.8: WENO numerical solutions of equation (2.12) with γ = 10−2 (left) and 10−4 (right) with source S(x) = φ(x − 10) and b = 0.0079 (dot-dot-dashed), 0.015 (dashed), 0.03 (solid). The bottom panels are the zoom in of the dashed square of top panels. Assume the source is S = Cφs (x), where φs (x) is the normalized frequency profile of the source photons. The analytic solution of equation (5.12) is ([31]), Z x+γτ −1 J(x, τ ) = J(x + γτ, 0) + Cγ φs (x′ )dx′ (5.13) x where J(x, 0) is the initial flux. Define the mean frequency by R xJ(τ, x)dx x¯(τ ) ≡ R . (5.14) J(τ, x)dx If the initial flux is J(0, x) = 0, one can show that for any profile φs (x) one has, x¯(τ ) = −γτ /2. (5.15) As expected, the speed of the redshift d¯ x/dτ is a constant. For the Hubble expansion, 56 a frequency shift x needs a time τ ≃ γ −1 x. Thus, in order to have the frequency shift x ≥ 3, the time scale is τ ≥ 6 × 106 , corresponding to t ≥ 4 × 104 years. This scale seems to be short enough compared with the lifetime of first stars. However, Hubble redshift is less effective compared with the resonant scattering. This point can be seen in figure 5.9, which gives the solution of equation (2.12) with the re-distribution function equation (2.7). The source is taken to be √ 2 S(x) = φg (x) = (1/ π)e−x (5.16) We also take the parameter γ to be 10−6 , and ignore recoil, i.e. b = 0. η = 10 7 6 10 η = 10 6 5 10 η = 10 5 104 η = 10 4 J(η,x) 3 10 η = 10 3 2 10 η = 10 2 1 10 0 10 -15 -10 -5 0 5 x Figure 5.9: Solutions J(τ, x) of equation (2.12) with the re-distribution function equation (2.7). Parameters b = 0, and γ = 10−6 . Figure 5.9 shows that the diffusion in the frequency space leads to a flat plateau with width |x| ≤ 3 when τ is as small as τ ≃ 104 , or t ≃ 102 years. This time scale is much less than that of the Hubble redshift. Therefore, the major mechanism for producing photons with frequency |x| ≥ 3 is given by the resonant scattering. 57 5.2.5 Bounce back mechanism From figure 5.9, we can see that the profile of photons are almost symmetric with respect to x = 0 (or to ν = ν0 ) until about the time τ = 106 . The redshift effect can only be seen from the curve of τ ≥ 106 . That is, the resonant scattering impedes cosmic redshift. The impediment is due to the “bounce back” of resonant scattering. Regardless of whether the frequency of the absorbed photons is larger or smaller than ν0 , the mean frequency of the re-emitted photons is always ν0 . Therefore, the resonant scattering will bring back some red shifted photons to the frequency ν0 . The net effect of resonant scattering (absorption and re-emission) is, on average, to bounce red shifted photons back to the resonant frequency ν0 , and to restore the symmetry with respect to x = 0. This bounce back mechanism is the key of restoring ν0 photons and the W-F coupling Figure 5.10: x ¯ vs. γτ for the solution shown in figure 5.9 (solid) and the solution of equation (5.15) (dashed). To illustrate the bounce back effect, we calculate the mean frequency x¯(τ ) for the solutions shown in figure 5.9. The result is plotted in figure 5.10. The solution for pure cosmic redshift equation (5.15) is also shown in figure 5.10 as the dashed 58 curve. We can see from figure 5.10 that the speed of photon frequency shift d¯ x/dτ when considering resonant scattering is in general less than in the case of pure cosmic redshift. When the time τ is large, many photons have been red shifted out of the frequency range of φg (x). In this case, the “bounce back” effect ceases, and redshift speed is recovered to d¯ x/dτ = −γ/2. From figure 5.10 one observes once again that the Hubble redshift can produce frequency shift from x = 0 to |x| ≃ 3 only when τ ≥ 107 , or t ≃ 105 years. This time scale may not be short enough to match the 21 cm region with a short lifetime. When the time scale of resonant scattering is less than the time scale of Hubble expansion, the frequency shift of the Hubble expansion slows down significantly due to the resonant scattering. 5.2.6 Discussion In this chapter we have concentrated on the formation of the local Boltzmann distri- bution in the photon frequency spectrum around the resonant frequency ν0 within width νl , i.e. |ν − ν0 | ≤ νl . The distribution takes the form, e−(ν−ν0 )/kT (5.17) We show that a local Boltzmann distribution will be formed if photons with fre- quency ∼ ν0 have undergone a O(104 ) scattering, which corresponds to the order of 103 yrs for neutral hydrogen density of the concordance ΛCDM model. The time evolution of the shape and width of the local Boltzmann distribution does not depend much on the details of atomic recoil, photon sources, or initial conditions. However, the intensity of the photon flux at the local Boltzmann distribution shows substan- tial time-dependence. The time scale of approaching the saturated intensity can be as long as 105 -106 yrs for typical parameters of the ΛCDM model. The intensity 59 of the local Boltzmann distribution at time less than 105 yrs is significantly lower than that of the saturation state. Therefore, it may not be reasonable to assume that the deviation of the spin temperature of 21 cm energy states from cosmic back- ground temperature is mainly due to the W-F coupling if the first stars or their emission/absorption regions evolved with a time scale equal to or less than Myrs. Chapter Six Radiative transfer in the 21 cm region 61 The 21 cm emission and absorption from gaseous halos around the first generation of stars substantially depend on the Wouthuysen-Field (W-F) coupling, which relates the spin temperature to the kinetic temperature of hydrogen gas via the resonant scattering between Lyα photons and neutral hydrogen. Therefore, the existence of Lyα photons in the 21 cm region is essential for this process. Although the center object generally is a strong source of Lyα photons, the transfer of Lyα photons in the 21 cm region is very inefficient, as the optical depth, as explained in equation (2.14), of Lyα photons is very large. Consequently, the Lyα photons from the source may not be able to efficiently transfer to the entire 21 cm region to provide the W-F coupling understanding. The W-F coupling is crucial to estimate the redshifted 21 cm signal from the halos of first generation of stars, because the deviation of the spin temperature from the temperature of cosmic microwave background (CMB) Tcmb is considered to be mainly caused by the W-F coupling ([13]). This problem is also particularly important considering that the lifetime of first stars generally is short. We investigate this problem with the numerical solution of the integrodifferential equation, which describes the kinetics of Lyα resonant photons in both physical and frequency spaces. 6.1 Basic properties of the 21 cm region The property of the ionized and heated regions around an individual luminous object is dependent on the luminosity (or mass), the spectrum of UV photon emission, and the time-evolution of the center object ([5], [23]). In the following work, We will not consider a specific model, but discuss only the common features. These halos generally consist of three spheres. The inner region of the halo is the highly ionized Str¨omgren sphere, or the HII region, in which the fraction of neutral hydrogen 62 fHI = nHI /nH is no more than 10−5 , where nHI and nH are, respectively, the number densities of neutral hydrogen HI and total hydrogen. The temperature of the HII region is about 104 K. The physical radius of this sphere is ranges from a few to a few tens of kpc (kilo parsec). The second region is the 21 cm emission shell, which is just outside the HII region. The physical size of this shell is similar to that of the HII sphere. The temperature of hydrogen gas in the emission shell is in the range 102 < T < 104 K due to the heating of UV photons and the fraction fHI is in the range of 0.1 to 1. The third region is the 21 cm absorption shell, which is outside the 21 cm emission region. The physical size is on the order of 100 kpc. The temperature of this region is lower than Tcmb = TCMB (1 + z), where TCMB is the temperature of CMB today and z is the measure of redshift. The time scale of the formation of the halos is about 106 years. The lifetime of the halos is about the same as the lifetime of the first stars. At the epoch of redshift z ≃ 20, the reionization region consists of isolated patches around the first sources. Most Lyα photons in the 21 cm regions are expected to come from the central source and the subsequent re-emission processes. If one can estimate the center object as a Lyα emitter, the emission of Lyα photons in number per unit time would be about dNLyα /dt = 1053 s−1 . The recombination of HII and electron in the Str¨omgren sphere is also a source of Lyα photons. If the physical radius of the Str¨omgren sphere is ∼ 10 kpc, the emission intensity is about dNLyα /dt = 1049 s−1 . Using the parameters of the concordance ΛCDM model, the optical depths of photons with Lyα resonant frequency ν0 in the 21 cm region are −1/2  3  Ωb h2    6 T 1+z R τ (ν0 ) = nHI σ0 R = 4.9 × 10 fHI , (6.1) 104 20 0.022 10kpc where σ0 is the cross section of the resonant scattering at the frequency ν0 , R is the physical size of the considered sphere, and τ is the distance R in the units of the 63 mean free path of ν0 photons at redshift z. Equation (6.1) shows that the optical depth of the 21 cm signal regions with fHI ≥ 0.1 is τ (ν0 ) ≥ 106 . It is also useful to define a dimensionless time η = cnHI σ0 t as  1/2  3   T 20 0.022 t = η/cnHI σ0 = 6.7 × −1 10−3 fHI η yrs. (6.2) 104 1+z Ωb h2 This is the time in the units of mean free flight-time. For a time scale t ≃ 105 yrs, the scale of η is ≃ 107 . 6.1.1 Problems with the W-F coupling There are, at least, two open questions regarding radiative transfer problems with the W-F coupling in the 21 cm regions: A.) How to provide enough Lyα photons in such an opaque medium? B.) Can the radiative transfer in the 21 cm region keep the W-F coupling working well? If photons always maintain the frequency ν0 , their spatial diffusion can be de- scribed as a random walk process ([6]). The mean number of scattering required for the diffusion over a range R is on the order of τ 2 . Thus, the time for the diffusion over the 21 cm range is η ≃ τ 2 , corresponding to the time t = τ 2 /cnHI σ0 ≥ 1012 years. This diffusion mechanism, obviously, is useless for the 21 cm regions. The time scale of the diffusion would be substantially reduced if the frequency of the Lyα photons can have a small shift from ν0 to ν = ν0 ± ∆ν, because the optical depth τ at the frequency ν0 ± ∆ν is significantly less than τ (ν0 ). One can estimate 64 the frequency shift ∆ν by the condition of optical depth τ (ν0 ± ∆ν) ≃ 1, which is φ(x, a) τ (ν0 ) ≤1 (6.3) φ(0, a) where φ(x, a) is the Voigt function for the resonant line ν0 profile as ([19]) ∞ 2 a e−y Z φ(x, a) = dy (6.4) π 3/2 −∞ (x − y)2 + a2 where the dimensionless variable x is defined by x = (ν −ν0 )/∆νD and ∆νD = ν0 vT /c p is the Doppler broadening of hydrogen gas with thermal velocity vT = kb T /2mH . Therefore, x measures the frequency deviation ∆ν = |ν − ν0 | in the units of the Doppler broadening. The parameter a in equation (6.4) is the ratio of the natural to the Doppler broadening. For Lyα line, a = 2.35 × 10−4 (T /104 )−1/2 . In terms of x, the solution of equation (6.3) is x ≥ 3. An effective mechanism of the frequency shift is provided by the diffusion of the photon distribution in the frequency space. Considering this mechanism, the number of scattering required for diffusion over a physical distance R is no longer equal to τ 2 , but roughly on the order of τ ([25], [14]). Since τ ∝ R, the time scale of the spatial diffusion over size R is comparable to R/c. Problem A would then be solved. Problem B still remains. If photons with the frequency ν < ν0 −∆ν or ν > ν0 +∆ν take a faster spatial diffusion than the ν0 photons, how can we maintain the W-F coupling? Without ν0 photons, one cannot have a local Boltzmann distribution around ν0 . Therefore, the photons with the frequency ν < ν0 − ∆ν or ν > ν0 + ∆ν should be brought back to the frequency ν0 . We need to study whether the frequency space diffusion mechanism can restore ν0 photons from photons with the frequency ν 6= ν0 . 65 Cosmic expansion also leads to a deviation of photon frequency from ν0 to a lower one ν0 − ∆ν, which speeds up the spatial transfer. However, we also need a mechanism to restore Lyα photons from Hubble redshifted photons. Therefore, in terms of the W-F coupling in the 21 cm region, we need to study the kinetics of Lyα photons in both physical and frequency spaces. 6.2 Radiative transfer (RT) of Lyα photons in the 21 cm regions Considering a photon source located at the central region of a uniformly distributed expanding medium, we can use the RT equation of the specific intensity I(η, r, x, µ) of the form, ∂I ∂I (1 − µ2 ) ∂I ∂I +µ + −γ = ∂η ∂r r ∂µ ∂x Z −φ(x)I + R(x, x′ )I(η, r, x′ , µ′ )dx′ dµ′ + S, (6.5) where µ = cos θ is the direction relative to the radius vector r. It is the same equation as the RT equation in a spherical halo (2.18). When the optical depth is large, we use the same Eddington approximation, I(η, r, x, µ) ≃ J(η, r, x) + 3µF (η, r, x) (6.6) 1 R +1 where J(η, r, x) = 2 −1 I(η, r, x, µ)dµ is the angularly averaged specific intensity 1 R +1 and F (η, r, x) = 2 −1 µI(η, r, x, µ)dµ is the flux. Defining j = r2 J and f = r2 F , 66 equation (2.19) yields the equations for j and f as ∂j ∂f ∂j Z + = −φ(x)j + R(x, x′ )jdx′ + γ + r2 S, (6.7) ∂η ∂r ∂x ∂f 1 ∂j 2 j + − = −φ(x)f. (6.8) ∂η 3 ∂r 3 r 6.2.1 Diffusion in the physical space We first solve equations (6.7) and (6.8) by dropping all terms related to the transfer in frequency. It has been shown that the source term can be replaced by a boundary condition of f and j at r = r0 ([26]). The equations then become ∂j ∂f + = 0, (6.9) ∂η ∂r ∂f 1 ∂j 2 j + − = −f. (6.10) ∂η 3 ∂r 3 r We take the boundary condition at r0 = 1 to be j(η, 1) = 3f (η, 1) = 3S0 . (6.11) Since equations (6.9) and (6.10) are linear, the parameter S0 is not important if we are only interested in the shape of j and f as function of η and r. In our calculation we use S0 = 106 in most cases. The initial condition is taken to be j(0, r) = f (0, r) = 0. (6.12) The solution of f (η, r) is presented in figure 6.1. The dashed line is f = S0 = 106 , which is the time-independently exact solution of equation (6.9). It is a horizontal straight line because the flux f = r2 F is r-independent and satisfies the conservation 67 of the photon number. Figure 6.1 shows that the spatial size r of f (η, r) roughly √ satisfies r ∝ η. Therefore, without resonant scattering, the diffusion basically is a random walk process as discussed in section 6.1. 6 10 105 η = 10 4 4 10 3 10 f(η,r) η = 10 3 102 η = 10 2 1 10 0 10 50 100 150 200 250 r Figure 6.1: The solution f (η, r) of equation (6.9) and (6.10). The straight line is f = S0 = 106 . We now turn to the solution of equations (6.7) and (6.8) with resonant scattering and Hubble redshift. The relevant parameters are given by b = 0 and γ = 10−5 . We use the boundary condition at r = 0 j(η, 0, x) = 0, f (η, 0, x) = S0 φs (x) (6.13) where the frequency profile is taken to be the Gaussian profile, i.e. φs (x) = φg (x). The parameter S0 is still taken to be 106 . The initial condition is similar to equation (6.12), i.e. j(0, r, x) = f (0, r, x) = 0. The solutions of f (η, r, x) are plotted in figure 6.2. All the solutions in figure 6.2 show two remarkable peaks at x ≃ ±(2 − 3) for all radius r. The amplitude of f (η, r, x) at the peaks is higher than at x = 0 by a factor of 10 to 102 . It shows that the flux is dominated by photons with frequency x ≃ ±(2 − 3). That is, the spatial transfer is carried out by photons of x ≃ ±(2 − 3). 68 The amplitude of the flux at the saturated peaks is basically r-independent, and R the saturated values of f (η, r, x)dx are r-independent. This is consistent with the conservation of the photon number. 6 6 10 10 5 5 10 η = 400 10 η = 800 η = 300 104 104 η = 600 η = 100 3 3 10 10 η = 200 2 2 10 10 f(η,r,x) f(η,r,x) 1 1 10 10 100 100 10-1 10-1 -2 -2 10 10 -3 -3 10 10 -4 -2 0 2 4 -4 -2 0 2 4 x x 6 10 5 10 104 3 η = 2000 10 η = 1500 2 10 f(η,r,x) η = 900 1 10 0 10 10-1 -2 10 -3 10 -4 -2 0 2 4 x Figure 6.2: The flux f (η, r, x) of the solutions of equation (6.7) and (6.8) at r = 50 (left), 100 √ 2 (middle) and 500 (right). The frequency profile of the source is φg (x) = (1/ π)e−x . More importantly, the amplitude of the peaks approaches its saturation at about η ≃ 200 for r = 50, η ≃ 400 for r = 100, and η ≃ 2, 000 for r = 500. That is, the time η needed for the spatial transfer over size r is roughly proportional to r. This is very different from the random walk relation η ∝ r2 , or the results shown in figure 6.1. 6.2.2 Resonant photon restoration and W-F coupling onset As mentioned in the previous sections the flux f (η, r, x) is lacking ν0 (or x = 0) photons even when f approaches its saturation. This is, obviously, not good for the 69 6 6 10 10 5 5 10 η = 400 10 η = 800 η = 300 104 104 η = 600 η = 100 3 3 10 10 η = 200 2 2 10 10 f(η,r,x) f(η,r,x) 1 1 10 10 100 100 10-1 10-1 -2 -2 10 10 -3 -3 10 10 -4 -2 0 2 4 -4 -2 0 2 4 x x 6 10 5 10 104 3 η = 2000 10 η = 1500 2 f(η,r,x) 10 η = 900 1 10 100 10-1 -2 10 -3 10 -4 -2 0 2 4 x Figure 6.3: The flux f (η, r, x) of the solutions of equation (6.7) and (6.8)pat r = 50 (left), 100 2 (middle) and 500 (right). The frequency profile of the source is φs (x) = (1/ π/2)e−2x . W-F coupling. However, we found that the situation may not be so if we consider the solution of the mean specific intensity j. Figure 6.4 presents the solution j(η, r, x) of equations (6.7) and (6.8) at r = 50, 100 and 500. The parameters b = 0 and γ = 10−5 are the same as the solutions in figure 6.2. The results are given in figure 6.4. When the time η is small, j(η, r, x) is similar to f (η, r, x), with a valley around x = 0 and two peaks at x ≃ ±(2 − 3). However, different from the flux f (η, r, x), the amplitude of j(η, r, x) around x = 0 is quickly increasing. In the saturated state it is about the same as the peaks i.e., although the flux always has a valley around x = 0, the ν0 photons are quickly restored in the mean intensity. The time scale for approaching its saturated state is proportional to r, not r2 . Therefore, the restoration of resonant photons is due to the resonant scattering “bounce back”, that pushes photons with frequency x ≃ ±(2 − 3) back to x ≃ 0. The shape of 70 η = 400 106 106 η = 800 η = 300 η = 600 10 5 105 j(η,r,x) j(η,r,x) η = 100 η = 200 104 104 3 3 10 10 -4 -2 0 2 4 -4 -2 0 2 4 x x η = 2000 5 10 η = 1500 104 η = 900 j(η,r,x) 3 10 2 10 101 -6 -4 -2 0 2 4 6 x Figure 6.4: The mean intensity j(η, r, x) of equation (6.7) and (6.8) at r = 50 (left), 100 (middle) and 500 (right). j(η, r, x) around x = 0 is a flat plateau, which is similar to that in figure 5.9. As have been shown in [27], the flat plateau of j(η, r, x) at b = 0 tends towards the local Boltzmann distribution if recoil is considered. We can expect that the flat plateau of figure 6.4 will also show a local Boltzmann distribution if b 6= 0. We calculate the solutions j(η, r, x) and f (η, r, x) of equation (2.22) and (2.23) at r = 500 with the same parameters as in figure 6.4, but with b = 0.3. We use a large b, because it is easier to see the slope of the local Boltzmann distribution. The results are given in figure 6.5. Figure 6.5 clearly shows a local Boltzmann distribution within the range |x| ≤ 2 as j(η, r, x) ≃ j(η, r, 0)e−2bx = j(η, r, 0)e−h(ν−ν0 )/kB T . (6.14) We find the slope to be ln j(η, 0) − ln j(η, 1) = 0.59, which is well consistent with 2b = 0.6. Figure 6.5 shows that at r = 500, the onset of the W-F coupling can 71 occur as early as η = 900, but the amplitude of the local Boltzmann distribution at that time is lower by a factor of 102 than its saturated value. The amplitude of the local Boltzmann distribution is substantially increasing with time. Therefore, the “bounce back” mechanism keeps the W-F coupling to work with a timescale η larger than a few hundreds, i.e. a few hundred collisions. This result is the same as that in [27]. 106 η = 2000 5 10 5 10 η = 1500 104 η = 2000 104 3 η = 900 10 10 3 η = 1500 f(η,r,x) j(η,r,x) 102 η = 900 102 1 10 1 10 100 100 -1 10 -1 10 10-2 10-2 -5 0 5 -5 0 5 x x Figure 6.5: The mean intensity j(η, r, x) (left) and flux f (η, r, x) (right) of the solution of equation (6.7) and (6.8) at r = 500. The recoil parameter is taken to be b = 0.3. The solution of the flux f (η, r, x) in figure 6.5 is not very different from that in figure 6.2. The only difference between the two figures is that the former is asymmetric with respect to x = 0, i.e. the peak at x < 0 is stronger than that of x > 0, while the latter is symmetric. This is due to the recoil of b 6= 0 inducing more photons to move to x < 0. Neither a flat plateau nor a local Boltzmann distribution is shown in the flux f (η, r, x). There is always a valley around ν0 even when f (η, r, x) is in its saturated state. It once again indicates that the flux is dominated by photons of x ≃ ±(2 − 3). 72 6.2.3 Effect of injected photons Hubble redshift is another mechanism to produce photons with frequency ∼ ν0 , which also have the problems of the spatial transfer and the W-F coupling. Since the 21 cm region basically is optical thin for photons with x > 3, we first model the cosmic redshift by a photon source with the profile φs (x) = φg (x − 3). The solutions of f (η, r, x) and j(η, r, x) with the same parameters as in figure 6.5 are shown in figure 6.6. 6.2.4 Effect of the Hubble redshift Although the cosmic expansion is considered in all the above-mentioned solutions, the effect of cosmic expansion seems to be negligible. This is because the optical depth is large and the parameter γ is very small. The number of resonant scattering within the Hubble time is very large. The cosmic expansion is too small in one free flight time and the “bounce back” is dominant. When fHI is smaller, the optical depth of the ν0 photons is smaller, and γ is larger, such that the effect of the Hubble redshift would appear. Figure 6.7 presents the solution f (η, r, x) and j(η, r, x) of equations (6.7) and (6.8) with the same parameters as those in figures 6.2 and 6.4, except that the parameter γ = 10−3 is large. Both j(η, r, x) and f (η, r, x) show the same features as in figures 6.2 and 6.4. The Hubble redshift makes the profile asymmetric with respect to x = 0. The red wing is stronger than the blue wing. j(η, r, x) still shows a flat plateau in the range −2 < x < 2. Hence, the W-F coupling work in case of γ = 10−3 , corresponding to fHI ≃ 10−3 . 73 6.3 Discussion and conclusions The kinetics of Lyα resonant photons in the HI media with high optical depth τ can basically be described as diffusion in both the physical space and the frequency space. If Lyα photons do not join the diffusion in the frequency space, the transfer of Lyα photons in the physical space is very inefficient, as the number of scattering needed for escape is proportional to τ 2 . The resonant scattering of Lyα photons and neutral hydrogen makes the diffusion processes in the physical space coupled to the diffusion processes in the frequency space. First, the diffusion in the frequency space provides a shortcut for the diffusion in the physical space. It makes the mean number of scattering for escape to be approximately proportional to τ . Second, the bounce back of resonant scattering provides a mechanism of quickly restoring ν0 photons from x ≃ ±(2 − 3) photons. Finally the W-F coupling is realized simultaneously with the restoration of the x = 0 photons. The mechanism of “escape via shortcut” plus “bounce back” is mainly carried out by the photons with frequency x ≃ ±(2 − 3). In a 21 cm emission region of physical size R ≃ 10 kpc and fHI > 0.1, the optical depth of x ≃ ±(2 − 3) photons is still larger than 1. Therefore, it is reasonable to use Eddington approximation. On the other hand, the optical depth of x ≃ ±(2 − 3) photons is much less than that of the x = 0 photons. The x ≃ ±(2 − 3) photons can transfer and enter the 21 cm emission region in a time scale less than 105 years. Therefore, the mechanism of “escape via shortcut” plus “bounce back” is able to timely support the W-F coupling of the 21 cm emission shell with Lyα photons from the center objects. The time dependence of the W-F coupling would make the 21 cm signals weaker than the predication given by models which do not consider this time-evolution. 74 Especially at the early stage of the formation of the 21 cm signal regions, the intensity of the local Boltzmann distribution is still very low, and therefore, one cannot assume that the spin temperature of 21 cm is locked to the kinetic temperature of gas. It may yield a low brightness temperature of the 21 cm signals. Although the mechanism of “escape via shortcut” plus “bounce back” helps Lyα diffusion, it does not mean that this mechanism will reduce the Gunn-Peterson opti- cal depth of the Lyα photons. On the contrary, the resonant scattering will lead to a slight increase of the optical depth of the Lyα in the 21 cm region, as the resonant scattering impedes the cosmic redshift (Figure 5.10). Consequently, there should be no observable redshifted optical signal with the frequency ν0 /(1 + z) to be spatially correlated with the (1 + z) redshifted 21 cm signal. The evolution of photons described by equations (6.7) and (6.8) conserves photon numbers. The number of Lyα photons is basically conserved if one can ignore the Lyα photon destruction processes, such as the two-photon process ([25]). Thus, subsequent evolution of the Lyα photons in the 21 cm region is to diffuse to a large sphere around the first stars. At the same time, the Lyα photons will be redshifted. When the redshift is large enough, their Gunn-Peterson optical depth will be small, and finally these photons will escape from the halo ([39]). The escaping sphere should be larger than the size of the 21 cm region. Therefore, redshifted Lyα optical signal with low surface brightness may come from a big halo around the 21 cm emission region. 75 5 10 η = 400 5 10 η = 300 η = 100 4 10 f(η,r,x) f(η,r,x) 4 10 η = 400 103 η = 300 η = 100 102 103 -4 -2 0 2 4 -2 0 2 x x 5 5 10 10 η = 800 η = 600 η = 200 4 4 10 10 f(η,r,x) j(η,r,x) η = 800 103 10 3 η = 600 η = 200 102 2 10 -4 -2 0 2 4 -4 -2 0 2 x x 105 10 5 η = 2000 η = 1500 4 η = 900 10 4 10 f(η,r,x) j(η,r,x) 3 10 3 10 η = 2000 η = 1500 2 10 η = 900 2 10 101 -4 -2 0 2 4 6 -4 -2 0 2 x x Figure 6.6: The mean intensity f (η, r, x) and flux j(η, r, x) of equation (6.7) and (6.8) at r = 50 (top), 100 (middle) and 500 (bottom), b = 0.3. The frequency profile of the source is φs (x) = φg (x − 3) 76 6 7 10 10 10 6 η = 400 5 10 η = 400 η = 300 η = 300 10 5 η = 100 4 10 η = 100 f(η,r,x) j(η,r,x) 4 10 103 3 10 2 10 2 10 1 101 10 -4 -2 0 2 4 -4 -2 0 2 4 x x 7 10 105 106 η = 800 η = 600 η = 800 5 10 104 η = 600 η = 200 f(η,r,x) j(η,r,x) 104 η = 200 103 3 10 102 102 1 1 10 10 -4 -2 0 2 4 -4 -2 0 2 4 x x 7 6 10 10 106 105 5 η = 2000 4 10 10 η = 2000 η = 1500 3 104 10 f(η,r,x) j(η,r,x) η = 1500 3 2 10 10 η = 900 η = 900 1 102 10 1 0 10 10 100 10-1 -1 -2 10 10 -5 0 5 -5 0 5 x x Figure 6.7: The solutions f (η, r, x) and j(η, r, x) of equation (6.7) and (6.8) at r = 50 (top), 100 (middle) and 500 (bottom). γ = 10−3 . The parameters are the same as those in figures 6.2 and 6.4, except that γ = 10−3 . Chapter Seven Radiation from neutral Hydrogen and flash sources 78 In this chapter we look at the Lyα emission from high red shifted sources such as galaxies, quasars and gamma ray bursts. The initial task is to investigate the flux of Lyα photon emergent from an optical halo containing a central light source. The goal is to show that even though the mean intensity of radiation is highly time dependent the frequency distribution of photon in the halo approaches a locally thermalized state around the resonant frequency. In the later part of this chapter we study the radiative transfer when the light emitted from the central source is a flash. We show that the nature of the transfer is a diffusion process, although it is not a purely Brownian diffusion. This leads to the trapping and storage of photons thermalized around the Lyα frequency for a long time after the cease of the central source emission. 7.1 Solutions of time-independent sources 7.1.1 Optically thick halos The property of the halo around individual luminous object depends on the lumi- nosity, the spectrum of UV photon emission, and the time-evolution of the center object. For luminous objects at high redshift, the halos generally consist of three spheres ([5], [23]). The most inner region is the highly ionized Str¨omgren sphere, or the HII region, which is optically thin for Lyα photons. The second region, which is just outside the HII region, is optically thick for Lyα photons. The temperature of the baryon gas is about 104 K, which is due to the heating of the UV photons. The third region is outside of the heated region. It is un-heated, and therefore, the temperature of the baryon gas can be as low as 102 K. 79 In this context, we study the transfer of Lyα photons in a radius R spherical halo of neutral hydrogen with temperature T in the range of 102 to 104 K. Assuming the uniformly distributed HI gas has number density nHI , the optical depth over a light path dl is dτ (ν) = nHI dl. σ(ν) is the cross section of the resonant scattering of Lyα photons by hydrogen, given as σ(x) = σ0 (∆νD )−1 φ(x, a) (7.1) where x is the dimensionless frequency defined by x ≡ (ν − ν0 )/∆νD , ν0 being the resonant frequency. ∆νD = ν0 (vT /c) is the Doppler broadening, and vT = p 2kB T /m. Here, m is the mass of the atom and kB is the Boltzmann costant. Therefore, x measures the frequency deviation ∆ν = |ν − ν0 | in units of ∆νD = 1.06 × 1011 (T /104 )1/2 Hz. σ0 = πe2 f /me c = 1.10 × 10−2 cm2 is the cross section of the resonant scattering at frequency ν0 = 2.46 × 1015 s−1 . The function φ(x, a) in equation (7.1) is the normalized profile given by the Voigt function as ([19]) ∞ 2 a e−y Z φ(x, a) = dy . (7.2) π 3/2 −∞ (x − y)2 + a2 The parameter a in equation (7.2) is the ratio of the natural to the Doppler broad- ening. Natural broadening is due to the uncertainty of the energy due to the finite lifetime of an excited state, for e.g. due to spontaneous radioactive decay. Doppler effect is due to the fact that when atoms in a gas emit radiation, each photon is known to have a distribution of velocities. Depending on the velocity the emitted photon are either red or blue shifted. For the Lyα line, a = 2.35 × 10−4 (T /104 )−1/2 . The optical depth of the halo with column number density of neutral hydrogen NHI = nHI R is τ (x) = NHI σ(x) = τ0 φ(x), (7.3) 80 where  −1/2   −1 7 T NHI τ0 = NHI σ0 (∆νD ) = 1.04 × 10 . (7.4) 104 1020 cm2 7.1.2 Radiative transfer equation in spherical halo The radiative transfer of Lyα photons in spherical halo is described by the equation of the specific intensity I(η, r, x, µ) as ∂I ∂I (1 − µ2 ) ∂I ∂I +µ + −γ = ∂η ∂r r ∂µ ∂x Z −φ(x; a)I + R(x, x′ ; a)I(η, r, x′ , µ′ )dx′ dµ′ + S, (7.5) In this case we are using the source of photons and the redistribution function as in equations (2.4) and (2.5). Using the Eddington approximation and substitutions outlined in chapter 2 we get a set of equations, ∂j ∂f ∂j Z + = −φ(x; a)j + R(x, x′ ; a)jdx′ + γ + r2 S, (7.6) ∂η ∂r ∂x ∂f 1 ∂j 2 j + − = −φ(x; a)f. (7.7) ∂η 3 ∂r 3 r The mean intensity j(η, r, x) describes the x photons trapped in the halo by the resonant scattering, while the flux f (η, r, x) describes the photons in transit. Without resonant scattering, equations (7.6) and (7.7) yield ∂j ∂f ∂j + = −φ(x; a)j + γ + r2 S, (7.8) ∂η ∂r ∂x ∂f 1 ∂j 2 j + − = −φ(x; a)f. (7.9) ∂η 3 ∂r 3 r 81 For photons with the frequency shifted away from ν0 , the halo would be optically thin, and therefore, the Eddington approximation will no longer be proper when the variable x is large. However we are mainly interested in the profile around x = 0, hence the Eddington approximation is applicable. The source term S in the equations (7.6) and (7.7) can be described by a boundary condition of j and f at r = r0 . We can take r0 = 0, as r0 is not important if the optical depth of the halo is large. Thus, we have j(η, 0, x) = 0, f (η, 0, x) = S0 φs (x), (7.10) where S0 , and φs (x) are, respectively, the intensity and normalized frequency profile of the sources. Since equation (7.6) - (7.9) are linear, the intensity S0 can be taken as any constant. That is, the solution f (x) of S0 = S is equal to Sf1 (x), where f1 (x) is the solution :f S = 1. On the other hand, the equations (7.6) and (7.7) are not linear with respect to the function φs (x), i.e. the solution f (x) of φs (x) is not equal to φs (x)f1 (x), where f1 is the solution of φs (x) = 1. In the range outside the halo, r > R, no photons propagate in the direction µ < 0. R −1 Therefore, the boundary condition at r = R given by 0 µJ(η, R, x, µ)dµ = 0 is then ([36]) j(η, R, x) = 2f (η, R, x). (7.11) There is no photon in the field before t = 0. Therefore, the initial condition is j(0, r, x) = f (0, r, x) = 0. (7.12) 82 7.2 Solutions of time-independent sources We will present some results of the solution of equations (7.6) and (7.7). The source of the photon is, √ 2 φs (x) = (1/ π)e−x (7.13) 7.2.1 Time scale of local thermalization 100 100 η = 500 10-1 10-1 η = 300 η = 500 η = 300 η = 200 10-2 10-2 f(η,r,x) j(η,r,x) η = 200 10-3 10-3 10-4 10-4 -4 -2 0 2 4 -4 -2 0 2 4 x x 100 100 η = 500 η = 300 10-1 10-1 η = 200 η = 500 10-2 10-2 f(η,r,x) j(η,r,x) η = 300 η = 200 10-3 10-3 10-4 10-4 -4 -2 0 2 4 -4 -2 0 2 4 x x Figure 7.1: Top two panels are the solutions of j(η, r, x) (left) and f (η, r, x) (right) of equations (7.6) and (7.7) at r = R = 102 and η = 200, 300, 500 with the boundary condition equation (7.11). Bottom two panels are the solutions without the boundary condition equation (7.11). The source is taken to be (7.13 )with S0 = 1 . The parameters are taken to be a = 10−3 and γ = 0−5 . 83 7.2.2 Restoring ν0 photons Figure 7.1 shows that the flux has a valley at x = 0, while the mean intensity j at time η = 500 is a rather high at x = 0. This indicates that ν0 photons are restored by the resonant scattering. According to the re-distribution function R(x, x′ ), the probability of transferring a x′ photon to a |x| < |x′ | photon is larger than that from x′ to |x| > |x′ |. Therefore, the net effect of resonant scattering is to bring photons with frequency x′ 6= 0 to the central range x ∼ 0 of frequency space. Photons at x ∼ 0 are then effectively restored. The top right panel of figure 7.1 shows that at η = 500 the flux f at x = 0 can be as large as 1/10 of the flux at its two peaks of x ≃ ±(2 − 3) which would be observable. 1010 1010 0 0 10 10 f(η,r,x) X 10 10 j(η,r,x) X 10 10 10-10 10-10 -20 -20 10 10 -4 -2 0 2 4 -4 -2 0 2 4 x x Figure 7.2: Solutions of j(η, r, x) (left) and f (η, r, x) (right) of (7.8) and (7.9) at r = 102 and η = 200, 300, 500. The source and all parameters are the same as in figure 7.1. As a comparison to figure 7.1, we compute the mean intensity j(η, r, x) and flux f (η, r, x) with equations (7.8) and (7.9). The result is shown in figure 7.2, with the parameters are the same as those in figure 7.1. Both j or f in figure 7.2 are very different from those in figure 7.1. The profiles in figure 7.2 are typical absorption spectra, with a very deep valley at x = 0. The absorption of ν0 photons at r = 102 (or τ0 = 102 ) is e−100 = 10−43 , which is the amount of j and f shown in figure 7.2. 84 The curves in figure 7.2 do not show any evolution in the time range η = 200−500. This means that j and f at the radius r = 102 are already stable at the time η < 200. Therefore, the time-evolution of j and j after η = 200 shown in figure 7.1 is fully caused by the resonant scattering. The thermalization and the restoration of ν0 photons are from this delayed evolution. 7.2.3 Two peaks in the flux profile The profiles of both j and f shown in figure 7.1 have two peak structures. The flux f is dominated by photons with frequency x± ≃ ±(2 − 3). The two peak structure has been recognized in all the time-independent solutions of the Fokker- Planck approximation ([14], [24], [10]), and Monte Carlo simulation ([22], [39], [1], [3], [37]. The point we want to emphasize is that this structure is independent of the profile of the source S. It is because the initial properties of the central sources have been forgotten during the local thermalization. 100 10-1 η = 500 10-2 f(η,r,x) η = 300 η = 200 -3 10 10-4 -4 -2 0 2 4 x Figure 7.3: Flux f (η, r, x) at r = 102 at radius 102 and time η = 200, 300 and 500. The p r =−2x 2 frequency profile of the source is φs (x) = (1/ π/2)e . Other parameters are the same as in figure 7.1 85 Figure 7.1 presents the flux f given by equations (7.6) and (7.7) with the same pa- √ 2 rameters as the solution shown figure 7.1, and the source profile φs (x) = (1/ π)e−x 2 p replaced by φs (x) = (1/ π/2)e−2x . That is, the line width of the new one is equal 1 to 2 of that in figure 7.1. Figure 7.3 plots the solution of f , which is identical with the bottom-right panel of figure 7.1. It shows that the initial distribution of the photon frequency is already forgotten, and therefore, the flux of figures 7.1 and 7.2 at η > 200 are from the same thermalized sources. The flux of figures 7.1 and 7.2 √ 2 are also maintained if the line width is broader than φs (x) = (1/ π)e−x . We will show this to be true even when the source has a continuous spectrum. Since the shape of the local thermalized distribution is time-independent, the frequency of the two peaks, |x± |, at a given r is also time-independent. When the photons of the flux f mainly come from the local thermalized photons, the frequency of the two peaks, |x± |, should not be described by the relation |x± | = (aτ )1/3 , because this relation is from a solution of the time-independent Fokker-Planck equation, which does not describe the thermalization ([24]). 5 4.5 4 3.5 3 2.5 x+ 2 1.5 1 0 1 2 10 10 10 a*r Figure 7.4: The positions of the peaks |x+ | as a function of (ar)1/3 for the solutions of equations (7.6) and (7.7) with a = 10−2 . Other parameters are the same as in figure 7.1. The dashed line is for x± = ±(ar)1/3 . 86 Figure 7.4 presents the peak frequency |x± | as a function of ar for solutions given by equations (7.6) and (7.7) with a = 10−2 . We consider only r ≥ 102 , as in the case of r ≤ 102 photons do not undergo a large enough number of scattering, and therefore, are not completely thermalized. With the dimensionless variables, ar is equal to aτ . A line |x± | = (ar)1/3 is also shown in figure 7.4. It shows that our numerical result of |x± | is significantly different from the (aτ0 )1/3 -law if ar < 30. That is, in the range ar < 30, photons of the flux f are dominated by the local thermalized photons. The frequency x± actually is dependent on the width of the flat plateau or the local thermalized distribution of j. Therefore, x± is always larger than two. The curve of figure 7.4 is approximately available for a = 10−3 and 10−4 . Thus, |x± | is larger than (ar)1/3 if r < 3 × 105 and 3 × 106 for a = 10−3 and 10−4 , respectively. Figure 7.4 shows a slow increase of |x± | with r in the range ar ≤ 30, and then, it approaches (ar)1/3 for larger ar. When r is large, more photons of the flux are attributed to the resonant scattering by the Lorentzian wing of the Voigt function. |x± | is then approaching to (ar)1/3 . It should be emphasized once again that all these results are independent of the intrinsic width of Lyα emission from the central source. 7.2.4 Absorption spectrum If the radiation from the sources has the continuous spectrum, the effect of neutral hydrogen halos is to produce an absorption line at ν = ν0 . The profile of the absorption line can also be found by solving equations (7.6) and (7.7), but replacing 87 the boundary equation (7.10) by j(η, 0, x) = 0, f (η, 0, x) = S0 . (7.14) That is, we assume that the original spectrum is flat in the frequency space. 101 101 η = 500 0 10 100 η = 300 η = 200 -1 10 η = 500 f(η,r,x) j(η,r,x) 10-1 η = 300 η = 200 -2 10 10-2 -3 10 10-3 -4 -2 0 2 4 -4 -2 0 2 4 x x Figure 7.5: Solutions of j(η, r, x) (left) and f (η, r, x) (right) of equations (7.6) and (7.7) at r = 102 and η = 200, 300 and 500. The source is given by 7.14. The parameters are a = 10−3 and γ = 10−5 . 7.2.5 The estimation of the HI column density As an application of the absorption spectrum in figure 7.5, we study the profile of the red damping wing of the optically thick intergalactic medium or IGM at high redshifts. Since the variable x is independent of the redshift, the profile of a red damping wing is directly given by the flux f (η, r, x) at the wing of low frequency x ≤ 0. Red damping wing is sometimes modeled with a damped Lyα system (DLA), which assumes that the damping wing profile is fully determined by the absorption of an optical thick halos around the host objects. The DLA profile of a red damping 88 0 10 10-1 -2 10 f(x) f(η,r,x) with η = 20000 f(x) with τ0 = 10 4 f(x) with τ0 = 6X10 3 -3 10 f(x) with τ0 = 4X10 3 10-4 5 10 -x Figure 7.6: Red damping wing (a) the DLA model τ0 = 4 × 103 (dotdot dashed line), 6 × 103 (dot dashed line) and 104 (dashed line); (b) the solution f (x) of (7.15) for τ0 = 104 (solid line). For both (a) and (b), the parameter a = 10−2 . wing is then given by f (x) = e−τ (x) (7.15) and x < 0, where τ (x) is from the Voigt function equation (7.2) as ∞ 2 a e−y Z τ (x) = τ0 φ(x, a) = τ0 dy . (7.16) π 3/2 −∞ (y − x)2 + a2 The column number density of neutral hydrogen atoms, NHI , is generally estimated by fitting the profile equation (7.16) with observation, and NHI can be found from τ0 by equation (7.4). 89 7.3 Modeling of flash sources 7.3.1 Frequency-dependence of Lyα transfer If the light source this significant time-dependence, like the optical afterglow of GRBs, we can treat the source as a flash described by a boundary conditions as:   S0 φs (x), 0 < η < η0  j(η, 0, x) = 0, f (η, 0, x) = (7.17)  0,  η > η0 . It describes a flash within a time range 0 < η < η0 , or the time duration is ∆η = η0 . We consider the case of r ≥ η0 . That is, the size of the halo is larger than the spatial range covered by the flash, as photon spatial transfer in optical thick medium cannot be faster than the speed of light. The initial condition retains as equation (7.10). Figure 7.7 presents two time-dependent solutions of the mean intensity j and flux f : one is for a flash equation (7.17) with η0 = 100 at r = 102 ; the other is for a flash with η0 = 500 at r = 103 . Other parameters are S0 = 1, a = 10−2 and √ 2 φs (x) = (1/ π) e−x . We still see the typical flat plateau of j in all times, even when the original time duration of the flash is as short as ∆η = 100. 7.3.2 The evolution of time duration of a flash To further study the time evolution, we plot figure 7.8, which gives the light curve, or the time-dependence of the flux f at x = 0. It shows the light curve of rising and decaying phases of the flux f . With these light curves, one can find the maximum of the flux f at time ηmax ; the rising phase is then η < ηmax , and decaying phase is η > 90 0.2 10-1 η = 500 0.15 η = 600 η = 200 η = 300 η = 300 η = 400 0.1 η = 500 η = 600 10-2 f(η,r,x) j(η,r,x) 0.05 -3 η = 200 10 -2 0 2 -6 -4 -2 0 2 4 6 x x 0.055 0.05 0.045 η = 6000 10 -2 0.04 η = 8000 η = 5000 0.035 η = 4000 f(η,r,x) j(η,r,x) 0.03 10-3 0.025 η = 2000 η = 3000 0.02 η = 4000 η = 5000 10 -4 η = 6000 -2 0 2 -4 -2 0 2 4 x x Figure 7.7: The profiles of j(η, r, x) and f (η, r, x) with time-dependent source [eq.(7.14)]. Top panel: j(η, r, x) (left) and f (η, r, x) (right) of η0 = 100 at r = 102 . Bottom panel: η0 = 500 at r = 103 . The parameters a = 10−2 and γ = 10−5 . ηmax . For each curve, one can define a time duration ∆η = η2 − η1 , where η1 < ηmax and η2 > ηmax , and both are given by the condition f (η1,2 , r, 0) = 0.9f (ηmax , r, 0). The time duration is then ∆η = η2 − η1 . The two top panels of figure 7.8 are for a flash source with original time duration ∆η = 50. Figure 7.8 shows that the time duration of the flash, ∆η, is r-dependent. At r = 0, i.e. at the source, ∆η is 50. At r = 102 (top left of figure 7.8), ∆η is about 200, while at r = 103 , we have ∆η ≃ 2000. Therefore, the time duration ∆η is increasing with r. This result shows that the time duration seems to be dependent mainly on r (or optical depth τ ) of the halo, regardless of the initial time duration 91 0.0016 8.5E-06 8E-06 0.0014 7.5E-06 0.0012 7E-06 6.5E-06 f(η,r,x) f(η,r,x) 0.001 6E-06 0.0008 5.5E-06 5E-06 0.0006 4.5E-06 200 300 400 500 600 700 800 2000 3000 4000 5000 6000 7000 8000 x x 8.5E-05 1.6E-05 8E-05 7.5E-05 1.4E-05 7E-05 6.5E-05 f(η,r,x) f(η,r,x) 1.2E-05 6E-05 5.5E-05 1E-05 5E-05 4.5E-05 8E-06 2000 3000 4000 5000 6000 7000 8000 2000 3000 4000 5000 6000 7000 8000 x x Figure 7.8: The time-dependence of f (η, r, x) at x = 0. Top panels: f (η, r, 0) at r = 102 (left) and r = 103 (right) for flash source with η0 = 50. Bottom panel: f (η, r, 0) at r = 103 of flash source with η0 = 100 (left) and η0 = 500 (right). ∆η = η0 . The bottom two panels of figure 7.4 are for flash sources with original time duration ∆η = η0 = 100 (left) and 500 (right), and r = 103 . The light curves of the bottom two panels f100 (η) (left) and f500 (η) (right) of figure 7.8 are about the same, i.e. f500 (η) ≃ 5f100 (η). Although the two flash sources have very different time durations at r = 0, the two light curves and time duration actually are about the same. The coefficient 5 is due to the total number of photons given by the flash η0 = 500 is 5 times larger than that of η0 = 100. This result shows again that the time duration is only dependent on the size r, regardless of the original time duration. 92 7.4 Discussion In this chapter we have used our algorithm to calculate the light curves of the flux from the halo. It shows that the flux is still a flash. The time duration of the flash for the flux, however, is independent of the original time duration of the light source but depends on the optical depth of the halo. Therefore, the spatial transfer of resonant photons is a diffusion process, even though it is not a purely Brownian diffusion. This feature indicates that the spatial transfer of a flash essentially is a diffusion process. Chapter Eight Ongoing work 94 The radiative transfer of Lyα photons in the spherical halo is described by the equation of the specific intensity I(η, r, x, µ) as ∂I ∂I (1 − µ2 ) ∂I ∂I +µ + −γ = ∂η ∂r r ∂µ ∂x Z −φ(x; a)I + R(x, x′ ; a)I(η, r, x′ , µ′ )dx′ dµ′ + S, (8.1) Our ongoing work is to solve this equation in three dimensions ie, in space, frequency and angular direction. The main idea is to solve the PDE without the Eddington approximation. The extra dimension on µ renders the computation significantly more intensive and challenging. At first we look at the steady state solution is by solving the equation, ∂I =0 ∂η As, ∂I (1 − µ2 ) ∂I ∂I µ + −γ = ∂r r ∂µ ∂x Z −φ(x; a)I + R(x, x′ ; a)I(η, r, x′ , µ′ )dx′ dµ′ + S, (8.2) In order to make the computation efficient, we reformulate the problem such that another remaining variable takes the role of ”time” so that the PDE can evolve on this variable. ∂I (1 − µ2 ) ∂I γ ∂I + − = ∂r µr ∂µ µ ∂x Z −φ(x; a)I + R(x, x′ ; a)I(η, r, x′ , µ′ )dx′ dµ′ + S, (8.3) 95 The integral term couples the values over the range of the variables µ and x, this implying we can only evolve in the variable r. In the computational sense this means that we treat the variable r as the time variable. There are certain assumptions required in a computation of this nature, most importantly we use the boundary condition in r as the initial condition. Our future goal is to compare the results of this equation with that of equations (2.21) and (2.22) obtained by performing the Ed- dington approximation. For the computation of equation (8.1) parallel computation is necessary. Chapter Nine Conclusion 97 In this work we have outilined a higher order algorithm used for solving a Boltzmann type equation with a stiff source term. This algorithm is suitable for problems which propagate in time with the solution profile exhibiting a sharp front. Using this algorithm we solve a type of Boltzmann equation, modeling transfer of photons and its scattering by hydrogen. The scattering process undergoes several processes which have different physical interpretations. The algorithm is used to successfully answer some details of these mechanisms such as • Estimating the time scale of a coupling procedure between the spin temper- ature and the kinetic temperature of the atom - Wouthuysen-Field coupling compared to the lifetime of a star • The dependence of the time scale on physical quantities such as recoil of an atom or the expansion of the universe • Calculation of radiation and intensity in the 21 cm region of halo of a star • Radiative field undergoing changes due to a flash-source The success of this algorithm is the stability that is achieved in computation of these various quantities in high resolution. We also addressed the issue of compu- tational intensiveness of these aspects by implementing certain techniques such as adaptive mesh routine and an integration algorithm that reduces the CPU time from an O(N 2 ) to an O(N ) computation. Our ongoing work is to formulate a high order fast algorithm to compute the solution of the three dimensional equation (2.18). This equation models the transfer of photons inside the halo of the star, where the Eddington approximation has not 98 been applied. The extra dimension renders the problem much more computationally intensive. The main goal in solving this problem is the parallel implementation of the numerical schemes and we hope to answer the relevant physical questions with highest precision. Bibliography [1] S-H Ahn, H.W. Lee, and J. M. Lee. Lyα line formation in starbursting galaxies. II. extremely thick, dustless, and static HI media. The Astrophysical Journal, 567:922–930, 2002. [2] M. M. Basko. The thermalization length of resonance radiation with partial frequency redistribution. The Astrophysical Journal, 17:125–139, 1981. [3] S. Cantalupo, C. Porciani, S. J. Lilly, and F. Miniati. Fluorescent lyα emission from the high-redshift intergalactic medium. The Astrophysical Journal, 628:61– 75, 2005. [4] J. Carrillo, I. Gamba, A. Majorana, and C.-W. Shu. A WENO-solver for the 1d non-stationary boltzmann-poisson system for semiconductor devices. Journal of Computational Electronics, 1:365–370, 2002. [5] R. Cen. Detection and fundamental applications of individual first galaxies,. The Astrophysical Journal, 648:47–53, 2006. [6] S. Chandrashekhar. Stochastic problems in physics and astronomy. Rev. Mod. Phys., 15:1–89, 1973. [7] X. Chen and J. Miralda-Escude. The spin-kinetic temperature coupling and the heating rate due to lyα scattering before reionization: Predictions for 21 centimeter emission and absorption. The Astrophysical Journal, 602:1–11, 2004. [8] L. Chuzhoy and P. R. Shapiro. Heating and cooling of the early intergalactic medium by resonance photons. The Astrophysical Journal, 655:843–846, 2006. [9] B. Cockburn, C. Johnson, C.-W. Shu, and E. Tadmor. Essentially non- oscillatory and weighted essentially non-oscillatory schemes for hyperbolic con- servation laws, volume 1697. Springer, 1998. [10] M. Dijkstra, Z. Haiman, and M. Spaans. Lyα radiation from collapsing proto- galaxies. i. characteristics of the emergent spectrum. The Astrophysical Journal, 649:14–36, 2006. [11] G. B. Field. Excitation of the hydrogen 21-cm line. Proceedings of the IRE, 46:240–250, 1958. 99 100 [12] G. B. Field. The time relaxation of a resonance-line profile. The Astrophysical Journal, 129:551–564, 1959. [13] S. R. Furlanetto and J. R. Pritehard. The scattering of lyman-series photons in the intergalactic medium. Monthly Notices of the Royal Astronomical Society, 372:1093–1103, 2006. [14] J. P. Harrington. The scattering of resonance-line radiation in the limit of large optical depth. Monthly Notices of the Royal Astronomical Society, 162:43–52, 1973. [15] L. G. Henyey. The doppler effect in resonance lines. Proceedings of the National Academy of Sciences, 26:50–54, 1941. [16] C. M. Hirata. Wouthuysen-field coupling strength and application to high- redshift 21-cm radiation. Monthly Notices of the Royal Astronomical Society, 367:259–274, 2006. [17] D. G. Hummer and G. B. Rybicki. The sobolev approximation for line format ion with partial frequency redistribution. The Astrophysical Journal, 387:248– 257, 1992. [18] G. B. Hummer. Non-coherent scattering: I. the redistribution function with doppler broadening. Monthly Notices of the Royal Astronomical Society, 125:21– 37, 1962. [19] G. B. Hummer. The voigt function: An eight-significant-figure table and gen- erating procedure. Memoirs of the Royal Astronomical Society, 70:1–32, 1965. [20] G. B. Hummer. Non-coherent scattering-vi. solutions of the transfer problem with a frequency-dependent source function. Monthly Notices of the Royal As- tronomical Society, 145:95–120, 1969. [21] G-S Jiang and C-W Shu. Efficient implementation of weighted ENO schemes. Journal of Computational Physics, 126:202–228, 1996. [22] J. S. Lee. Monte carlo simulation of emission frequencies from partial frequency redistribution functions. The Astrophysical Journal, 192:465–474, 1974. [23] J Liu, J-M Qiu, L-L Feng, C-W Shu, and L-Z Fang. 21 cm signals from early ionizing sources. The Astrophysical Journal, 663:1–9, 2007. [24] D. Neufeld. The transfer of resonance-line radiation in static astrophysical me- dia. The Astrophysical Journal, 350:216–241, 1990. [25] D. E. Osterbrock. The escape of resonance-line radiation from an optically thick nebula. The Astrophysical Journal, 135:195–211, 1962. [26] J.-M. Qiu, C.-W. Shu, L.-L. Feng, and L.-Z. Fang. A WENO algorithm for the radiative transfer and ionized sphere at reionization. New Astronomy, 12:1–10, 2006. 101 [27] I Roy, J-M Qiu, C-W Shu, and L-Z Fang. A WENO algorithm for radiative transfer with resonant scattering: the time scale of the wouthuysen-field cou- pling. New Astronomy, 14:513–520, 2009. [28] I. Roy, C.-W. Shu, and L.-Z. Fang. Resonant scattering and lyα radiation emer- gent from neutral hydrogen halos. The Astrophysical Journal, 2010. Submitted. [29] I Roy, W Xu, J-M Qiu, C-W Shu, and L-Z Fang. Time evoluti on of wouthuysen- field coupling. The Astrophysical Journal, 694:1121–1130, 2009. [30] I Roy, W. Xu, J-M Qiu, C-W Shu, and L-Z Fang. Wouthuysen-field coupling in 21 cm region around high redshift sources. The Astrophysical Journal, 703:1992– 2003, 2009. [31] G. B. Rybicki and I. Dell’Antonio. The time development of a resonance line in the expanding universe. The Astrophysical Journal, 427:603–617, 1994. [32] J. Shen, C-W. Shu, and M. Zhang. A high order WENO scheme for a hierarchical size-structured population model. Journal of Scientific Computing, 33:279–291, 2007. [33] C.-W. Shu. Numerical experiments on the accuracy of ENO and modified ENO schemes. Journal of Scientific Computing, 5:127–149, 1990. [34] C.-W. Shu and S. Osher. Efficient implementation of essentially non-oscillatory shock-capturing schemes. Journal of Computational Physics, 77:439–471, 1988. [35] L. Spitzer and J. L. Greenstein. Continuous emission from planetary nebulae. The Astrophysical Journal, 114:407–420, 1951. [36] W. Unno. Theoretical line contour of the lyman α radiation of the ionized helium and the excitation of bowen lines in planetary nebulae. Publication of the Astronomical Society, 7:81, 1955. [37] A. Verhamme, D. Schaerer, and A. Maselli. 3D lyα radiation transfer. - I. under- standing lyα line profile morphologies. Astronomy and Astrophysics, 460:397– 413, 2006. [38] S. A. Wouthuysen. On the excitation mechanism of the 21-cm (radio-frequency) interstellar hydrogen emission line. The Astronomical Journal, 57:31–32, 1952. [39] Z. Zheng and J. Miralda-Escude. Monte carlo simulation of lyα scattering and application to damped lyα systems. The Astrophysical Journal, 578:33–42, 2002.