University of Southern California Dissertations and Theses

Molecular dynamics and virial expansions for equilibrium thermodynamics in the solar interior

(USC Thesis Other)

Molecular dynamics and virial expansions for equilibrium thermodynamics in the solar interior

PDF

Download a page range

Download transcript

Copy asset link

Request this asset

Transcript (if available)

Content MOLECULAR DYNAMICS AND VIRIAL EXPANSIONS FOR EQUILIBRIUM THERMODYNAMICS IN THE SOLAR INTERIOR by Kathleen Anita Mussack A Dissertation Presented to the FACULTY OF THE GRADUATE SCHOOL UNIVERSITY OF SOUTHERN CALIFORNIA In Partial Fulfillment of the Requirements for the Degree DOCTOR OF PHILOSOPHY (PHYSICS) August 2007 Copyright 2007 Kathleen Anita Mussack Dedication For Grandpa Bob ii Acknowledgments I begin by expressing my heartfelt gratitude to my wonderful advisor, Dr. Werner D¨ appen. Your continual wisdom, guidance and inspiration have helped me to become a better scientist, teacher, and fellow human being. Thank you for being such a great mentor. I only hope that I can one day live up to your example. I thank Dr. Asher Perez for teaching me some really cool math and bringing me in on the equation of state project. I have really enjoyed working on it and I look forward to continuing this work with you. I thank USC’s Collaboratory for Advanced Computing and Simulations, especially Dr. Rajiv Kalia, Dr. Aiichiro Nakano, and Dr. Priya Vashishta. Each of you made time to teach Dan and me about simulations and to talk with us about our research. Your guidance has been invaluable. Thanks also to Cheng Zhang for your patience with all my questions and your fabulous MD tech support. Thank you to Dr. Stephan Haas, Dr. Rajiv Kalia, Dr. Joe Kunc, and Dr. Nick Warner for serving on my guidance committee. I also thank my colleagues Dan Mao and Amy Cassidy for countless research dis- cussions that provided new perspectives, inspiration, and encouragement. Last, but certainly not least, I thank my loving family and friends for your endless support and understanding. I could not have come this far without all of you. iii Contents Dedication ii Acknowledgments iii List of Tables vii List of Figures viii Abstract xi 1 Introduction 1 1.1 Properties of the Sun . . . . . . . . . . . . . . . . . . . . . . . . . . 3 1.2 Understanding the Sun . . . . . . . . . . . . . . . . . . . . . . . . . 4 1.3 Equations of State . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5 1.4 Nuclear Reactions and Neutrinos . . . . . . . . . . . . . . . . . . . . 6 I Molecular Dynamics 9 2 Electrostatic Screening 9 2.1 Salpeter Approximation for Electrostatic Screening . . . . . . . . . . 11 2.2 Deviations from Salpeter Screening . . . . . . . . . . . . . . . . . . . 12 2.3 Dynamic Effects . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 13 2.4 Shaviv and Shaviv . . . . . . . . . . . . . . . . . . . . . . . . . . . . 13 3 Molecular Dynamics 16 3.1 Velocity-Verlet Algorithm . . . . . . . . . . . . . . . . . . . . . . . . 16 3.1.1 Verlet Algorithm . . . . . . . . . . . . . . . . . . . . . . . . 16 3.1.2 Velocity-Verlet Algorithm . . . . . . . . . . . . . . . . . . . 18 3.1.3 Energy Conservation . . . . . . . . . . . . . . . . . . . . . . 19 3.2 Periodic Boundary Conditions and Minimum Image Convention . . . 21 3.3 Smooth and Continuous Potential . . . . . . . . . . . . . . . . . . . . 23 3.4 Dealing with Long-Range Forces . . . . . . . . . . . . . . . . . . . . 24 3.4.1 Tail Corrections . . . . . . . . . . . . . . . . . . . . . . . . . 25 iv 3.4.2 Ewald Sums . . . . . . . . . . . . . . . . . . . . . . . . . . . 26 3.4.3 Particle-Mesh Method . . . . . . . . . . . . . . . . . . . . . 28 3.4.4 Effective Cutoff Radius . . . . . . . . . . . . . . . . . . . . . 29 3.5 Quantum Effects . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 30 3.6 Thermalization and Time Averages . . . . . . . . . . . . . . . . . . . 32 4 Setting Up Our System 33 4.1 Units . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 33 4.2 Initializing Coordinates and Velocities . . . . . . . . . . . . . . . . . 34 4.3 Assigning Charge and Mass . . . . . . . . . . . . . . . . . . . . . . . 34 4.4 Temperature and Density . . . . . . . . . . . . . . . . . . . . . . . . 35 5 Our Preliminary Tests 37 5.1 Conservation of Energy and Velocity Distributions . . . . . . . . . . . 37 5.2 Pair Distribution Functions . . . . . . . . . . . . . . . . . . . . . . . 38 6 Our Evaluation of Shaviv and Shaviv’s Techniques 43 6.1 Effective Radius of Interaction . . . . . . . . . . . . . . . . . . . . . 44 6.2 System Size . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 46 6.2.1 Number of Dimensions . . . . . . . . . . . . . . . . . . . . . 46 6.2.2 Number of Particles . . . . . . . . . . . . . . . . . . . . . . . 46 6.3 Effective Radius of Interaction Revisited . . . . . . . . . . . . . . . . 48 6.4 Quantum Effects . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 48 6.5 Thermostat . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 51 6.5.1 Nos´ e-Hoover Thermostat . . . . . . . . . . . . . . . . . . . . 51 6.5.2 Thermostat Tests . . . . . . . . . . . . . . . . . . . . . . . . 54 6.6 Simple Thermostat . . . . . . . . . . . . . . . . . . . . . . . . . . . 59 6.7 Temperature . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 62 II Virial Expansions 64 7 Solar Equations of State 64 7.1 The Mihalas-Hummer-D¨ appen EOS (MHD) . . . . . . . . . . . . . . 65 7.2 OPAL . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 68 7.3 Planck-Larkin Partition Function . . . . . . . . . . . . . . . . . . . . 70 7.4 Feynman-Kac . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 70 8 Path-integral Formalism in the Feynman-Kac Representation 72 8.1 Grand-Canonical Partition Function . . . . . . . . . . . . . . . . . . 73 8.2 Single Particle in the FK Representation . . . . . . . . . . . . . . . . 74 8.3 Maxwell-Boltzmann Statistics for SystemS . . . . . . . . . . . . . . 77 8.3.1 Maxwell-Boltzmann Grand Partition Function . . . . . . . . . 77 v 8.3.2 Virial Expansion for the Coulomb Potential . . . . . . . . . . 79 8.4 Quantum Statistics . . . . . . . . . . . . . . . . . . . . . . . . . . . 82 8.4.1 Quantum Grand-Canonical Partition Function . . . . . . . . . 82 8.4.2 Exchange Contributions . . . . . . . . . . . . . . . . . . . . 84 8.5 Density Expansions . . . . . . . . . . . . . . . . . . . . . . . . . . . 85 8.5.1 Equation of State up to Orderρ 5/2 . . . . . . . . . . . . . . . 86 9 Our Analysis of Applicability in the Solar Interior 89 9.1 Coupling . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 90 9.2 Degeneracy . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 91 9.3 Ionization and Recombination . . . . . . . . . . . . . . . . . . . . . 91 9.4 Domain of Applicability . . . . . . . . . . . . . . . . . . . . . . . . 94 10 Results of Application to Solar Interior 95 10.1 Developing the Model . . . . . . . . . . . . . . . . . . . . . . . . . . 96 10.2 Improving the Precision of the Model . . . . . . . . . . . . . . . . . 97 10.3 Applying the Formalism . . . . . . . . . . . . . . . . . . . . . . . . 99 10.4 Results for Pressure . . . . . . . . . . . . . . . . . . . . . . . . . . . 100 10.5 Free Energy, Internal Energy, and Entropy . . . . . . . . . . . . . . . 100 10.6 Results for Specific Heats . . . . . . . . . . . . . . . . . . . . . . . . 103 10.7 Results for the Adiabatic Exponents . . . . . . . . . . . . . . . . . . 105 11 Extensions of EOS 111 11.1 Relativistic Correction . . . . . . . . . . . . . . . . . . . . . . . . . 111 11.2 Saha for Ionization . . . . . . . . . . . . . . . . . . . . . . . . . . . 112 11.3 More Chemical Elements . . . . . . . . . . . . . . . . . . . . . . . . 113 11.4 Extension to Orderρ 3 . . . . . . . . . . . . . . . . . . . . . . . . . . 114 11.5 Stellar Modeling . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 114 Conclusion 115 Bibliography 117 Appendix 122 vi List of Tables 1.1 Currently accepted values for solar properties. . . . . . . . . . . . . . 4 2.1 Relevant Reactions: Gamov energy and the width of the Gamov peak with respect to the thermal energy. . . . . . . . . . . . . . . . . . . . 14 3.1 Thermal deBroglie wavelenghts for protons and electrons at a temper- ature of 47.5 E H . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 31 10.1 Mass, charge, and spin for the elements in our EOS model. . . . . . . 96 10.2 New values for the virial expansion constants. . . . . . . . . . . . . . 99 vii List of Figures 1.1 Pressure waves in the Sun reflect off the outer layers and refract as they travel through the deeper layers. (Figure provided by A.G. Koso- vichev). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5 2.1 Plasma polarized by ion. . . . . . . . . . . . . . . . . . . . . . . . . 9 2.2 Approaching ions are each surrounded by their own screening cloud. . 10 2.3 Comparison between the Coulomb potential and the screened poten- tial. Particles with a given relative energy can approach closer when subject to the screened potential than the Coulomb potential. . . . . . 10 2.4 Maxwell-Boltzmann velocity distribution for protons in the solar core (velocity is in atomic units). . . . . . . . . . . . . . . . . . . . . . . 14 3.1 Periodic boundary conditions: the simulation box (bold) is surrounded by an infinite periodic lattice of identical boxes. When a particle leaves one side of the simulation box, its image from a neighboring box enters on the opposite side. . . . . . . . . . . . . . . . . . . . . . . . . . . 21 3.2 The Lennard-Jones potential. . . . . . . . . . . . . . . . . . . . . . . 24 3.3 The Coulomb potential. . . . . . . . . . . . . . . . . . . . . . . . . 25 3.4 Representation of the Ewald sum method: each point-particle’sδ-function is represented as (a) aδ-function minus a gaussian (b) plus a gaussian. 27 3.5 Without quantum corrections, the energy jumps dramatically when particles get too close. . . . . . . . . . . . . . . . . . . . . . . . . . 30 4.1 Maxwell-Boltzmann velocity distribution for a temperature of 47.5 E H for (a) protons and (b) electrons. . . . . . . . . . . . . . . . . . . . . 36 5.1 Conservation of energy for the pilot MD system. . . . . . . . . . . . 38 viii 5.2 Speed distribution for electrons compared to the Maxwell-Boltzmann distribution. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 39 5.3 Speed distribution for protons compared to the Maxwell-Boltzmann distribution. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 39 5.4 Pair correlations for the same system at different stages of thermaliza- tion. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 41 5.5 Pair correlations of the MD system compared to Debye-based pair cor- relations (lines). . . . . . . . . . . . . . . . . . . . . . . . . . . . . 42 6.1 Only particles within a distance less than the cutoff radius are included in the potential sum. . . . . . . . . . . . . . . . . . . . . . . . . . . 44 6.2 Average potential energy per particle computed using varying cutoff radii for the Coulomb sums. Here, the Debye length is 1.94 Bohr radii. 45 6.3 Pair correlation functions for a two-dimensional system and a three- dimensional system. . . . . . . . . . . . . . . . . . . . . . . . . . . 47 6.4 Pair correlation functions for 2-D systems with 100 and 2500 particles. 47 6.5 Pair correlation functions for 2-D systems with 2500 particles and two different cutoff radii (5 and 7.07 Bohr radii). . . . . . . . . . . . . . 48 6.6 The energy is unstable in a classical system. . . . . . . . . . . . . . . 49 6.7 When quantum corrections are included, the energy is conserved. . . 49 6.8 Pair correlation functions for systems with and without quantum ef- fects included in the potentials. . . . . . . . . . . . . . . . . . . . . 50 6.9 Pair correlation functions from Debye theory compared to those from molecular-dynamics systems with and without quantum effects included in the potentials. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 51 6.10 Pair distribution functions for a micro-canonical system and a canoni- cal system with a Trotter-expansion based thermostat. . . . . . . . . 58 6.11 Pair correlation functions for a simple thermostat and for a Trotter ther- mostat. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 61 6.12 Comparison of pair correlation functions for systems of temperature 19.3 E H and 47.5 E H . . . . . . . . . . . . . . . . . . . . . . . . . . . 63 ix 8.1 A classical closed filament located at ~ r with a shape parameterized by ~ ξ(s). This represents the one-body diagonal matrix contribution. (Figure provided by A. Perez) . . . . . . . . . . . . . . . . . . . . . 76 8.2 Quantum point-particle systemS and the equivalent classical system of filamentsS ∗ . (Figure provided by A. Perez) . . . . . . . . . . . . 79 9.1 Coupling parameter throughout the Sun. . . . . . . . . . . . . . . . . 90 9.2 Degeneracy throughout the Sun for (a) electrons, (b) protons, and (c) helium nuclei. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 91 9.3 Hydrogen ionization throughout the Sun. . . . . . . . . . . . . . . . 92 9.4 Helium ionization throughout the Sun. . . . . . . . . . . . . . . . . . 93 9.5 Oxygen ionization throughout the Sun. . . . . . . . . . . . . . . . . 93 9.6 log(T) for each depth r/R in the Sun. . . . . . . . . . . . . . . . . . 93 10.1 Results for the pressureβP/n at each depth in the Sun. . . . . . . . . 100 10.2 Results for c v at each depth in the Sun. . . . . . . . . . . . . . . . . 105 10.3 Results forχ T at each depth in the Sun. . . . . . . . . . . . . . . . . 106 10.4 Results forχ ρ at each depth in the Sun. . . . . . . . . . . . . . . . . 106 10.5 Results for c p at each depth in the Sun. . . . . . . . . . . . . . . . . 107 10.6 Results forΓ 1 at each depth in the Sun. . . . . . . . . . . . . . . . . 108 10.7 Results for∇ ad at each depth in the Sun. . . . . . . . . . . . . . . . . 109 10.8 Results forΓ 2 at each depth in the Sun. . . . . . . . . . . . . . . . . 109 10.9 Results forΓ 3 at each depth in the Sun. . . . . . . . . . . . . . . . . 110 x Abstract In the deep solar interior, matter is mostly fully ionized. Precise helioseismic and neutrino data provide an incredibly accurate description of the Sun. The established MHD and OPAL equations of state do not match the solar data to the available pre- cision. A rigorous quantum-statistical formalism for Coulomb systems is used to de- velop new tools for equation of state modeling. Although the formalism is in principle a low-density development, the relevant parameters allow the formalism to be applied in the solar center. Solar models can also be improved through the study of nuclear reaction rates. The Coulomb potentials in the ionized plasma are screened, effectively truncating the interactions. This screening makes it easier for ions to tunnel through the potential barrier. In a seminal paper, Salpeter (1954) discussed this enhancement. He based his study on the approximation of a static screening potential. There is a legitimate concern that dynamic effects could alter the result. Since 1996, Shaviv and Shaviv have been examining this question by numerical molecular-dynamics simu- lations. The calculation is essentially classical, although electrons are treated with effective potentials to mimic some quantum corrections. Shaviv and Shaviv have re- ported dynamical effects for the fast protons which are the ions that are most likely to engage in nuclear reactions. Given the importance of these effects for solar and stellar modeling, molecular-dynamics simulations have been developed to verify the results of Shaviv and Shaviv in an independent calculation. These studies of hot, dense plasmas will help develop more accurate solar models. xi Chapter 1 Introduction The Sun is important to us, not only because of its myriad effects on our daily lives, but also because it provides us with a nearby laboratory in which we can study con- ditions that are not yet attainable on Earth. At its core we find an ionized plasma with a temperature of 15.7 million K and a density of 151 g/cc, conditions extreme enough for nuclear fusion. Since we are unable to go inside the Sun to directly mea- sure these properties, you may wonder how we can possibly know so much about it. As our nearest star, the Sun floods us with information about itself in the form of ra- diation, vibrations, and neutrinos. It is up to us to decode those messages in order to understand its inner workings. Helioseismology and solar neutrino studies provide in- credibly accurate solar data, demanding equally accurate solar theories. Current solar models are at least 90% accurate in describing the thermodynamics and fusion of the solar interior. The research presented in this thesis focuses on two different methods of examining the remaining subtleties of the equilibrium thermodynamics in the Sun in order to improve the accuracy of solar models. In the first approach, I use molecular dynamics techniques to address the question of whether the static approximation for Coulomb screening is sufficient to accurately 1 describe the Gamov tunneling probability, a crucial ingredient in nuclear-fusion cross sections. If the static approximation is inadequate, dynamic effects must be included. Pioneering work in this field was done by Shaviv and Shaviv [75, 76, 77, 78], who used the techniques of molecular dynamics to model the plasma of the solar core. Their simulations indicate that dynamic screening does affect solar nuclear reaction rates, however much controversy remains over their claims. In Part I, I develop molec- ular dynamics simulations of hot, dense, two-component Coulomb plasmas to test the work of Shaviv and Shaviv, analyze the assumptions inherent in their models, and ex- amine the influence of the simulation techniques they chose to implement. With this approach, I examine screening and screening enhancement of nuclear reactions in the solar core. This enables the evaluation of possible dynamic effects that would alter the widely accepted static approximation. In the second approach, I use virial expansions based on Feynman-path-integral methods to develop a solar equation of state with the systematic density expansion of the equilibrium thermodynamic functions presented by Alastuey and Perez [4, 2, 3, 5]. So far, solar models happily use the intuitive MHD and the fairly systematic OPAL models of the equations of state. However, discrepancies with observations remain, and more rigorous options are being pursued. While some of the expansion for- malisms, such as OPAL, choose to include only the graphs that give dominant con- tributions to the thermodynamic functions, the Alastuey-Perez formalism differs in that it includes the terms from all graphs up to orderρ 5/2 (including contributions from screening, long- and short-range quantum effects, and Fermi-Dirac and Bose-Einstein statistics) and is therefore exact up to that order. For example, diffraction terms are included here though they are usually left out of other treatments. In Part II, I analyze the applicability of this formalism in the Sun, determine the region in which its use is 2 appropriate, and then apply the treatment in the appropriate region in order to attain new high-precision equation of state diagnostics for solar modeling. I also describe ways to improve the model by adding such things as relativistic effects and embedding the Saha equation to compute the degree of ionization so that cooler regions can be included in the treatment. Both of these methods naturally lead to future extensions. There are several facets of the molecular dynamics simulations to examine such as seeking better ways to model two-component plasmas and to analyze dynamic effects within them. The equa- tion of state work can be adapted to include more chemical elements, and the calcula- tion can be extended to include the terms of order ρ 3 . In addition, this formalism can be used to examine implications of the new chemical abundances. Together, these two paths of research will lead to better, more accurate models of the solar interior. 1.1 Properties of the Sun Some of the macroscopic properties of the Sun can be obtained with the ubiquitous techniques of geometry, radar ranging, and gravitational analysis. First, the distance between the Sun and the Earth can be determined very accurately with radar echos between planets [81]. In the same way, the solar mass can be found using Kepler’s laws of planetary motion [81, 16]. Combining the known distance and angular size with a little geometry, we can find the solar radius [81, 14]. We can also determine the Sun’s luminosity from space-based radiance measurements, taking into account temporal variance due to solar cycles and the Sun’s geometrical asymmetry [86, 16]. With all this information, we can use Stefan-Boltzmann’s Law L = 4πR 2 σT 4 eff to find the effective temperature of the photosphere, where σ is the Stefan-Boltzmann constant. 3 Property Symbol Value mass M (1.9891± 0.0004)× 10 33 g radius R (6.9626± 0.0007)× 10 10 cm luminosity L (3.846± 0.01)× 10 33 erg/s surface temperature T eff 5777± 2.5K Table 1.1: Currently accepted values for solar properties. Finally, the chemical composition of the solar surface can be obtained through spectroscopy. The current solar surface abundance ratio in mass of heavy elements (ie: more massive than helium) to hydrogen is Z/X = 0.0245± 0.001, but it is slightly reduced to 0.0165± 0.0017 according to recent re-evaluation by Lodders and Asplund et al [46, 8]. Although helium was first discovered in and named after the Sun, the spectrum of atomic helium is too faint and complicated to determine the helium abun- dance from spectroscopy. Instead it must be determined through indirect methods such as helioseismic inversions. 1.2 Understanding the Sun Just as seismologists use virbations in the Earth to learn about what lies beneath the surface, helioseismologists use virbrations in the Sun to study its internal composition. Pressure waves, or p-waves, are standing waves that propagate through the solar media, leaving compression patterns on the surface to be detected with dopplergrams. These waves are believed to originate in the resonance cavity beneath the photo- sphere [19, 66]. As they approach the surface of the Sun, the density of the material through which they travel drops dramatically. Eventually the waves lose momentum and are reflected back into the Sun. As the waves penetrate deeper into the Sun, the temperature increases dramatically. Thus the sound speed rises sharply, gradually re- 4 Figure 1.1: Pressure waves in the Sun reflect off the outer layers and refract as they travel through the deeper layers. (Figure provided by A.G. Kosovichev). fracting the waves. Eventually they reach a turning point and begin propagating up- wards again (see Fig.1.1). Each wave is affected by the material in the regions through which it travels. Waves of different frequencies dive to different depths, probing differ- ent regions of the Sun. Therefore the combined behavior of all of the waves provides a detailed description of the solar interior. 1.3 Equations of State Equations of state are thermodynamic equations used to describe the relationships be- tween state variables in matter. The type of matter determines what physical contri- butions need to be included in the equation. In an ideal gas, the physics is relatively simple, yielding the ideal gas law as the appropriate equation of state pV = Nk B T . (1.1) 5 For the plasma of the solar interior, we need to include contributions from many other physical phenomena such as the screening that arises from the polarization of the plasma, bound and scattering states, and diffraction and exchange effects from the quantum nature of the plasma. At such high temperatures, relativistic effects may also provide a correction for the electron behavior. If all of these effects are correctly accounted for, a solar equation of state is a powerful tool for understanding the physics of the solar interior. It can be used to compute relevant properties throughout the Sun. There are currently two different versions of widely accepted solar equations of state, OPAL and MHD (see chapter 7). Though these equations of state provide in- credibly reliable descriptions of the solar interior, there is still room for improvement. Part II of this thesis will discuss the development and application of a new equation of state based on virial expansions of the equilibrium thermodynamic functions. This formalism has the potential to provide a more accurate description of the fully ionized plasma in the core of the Sun, improving our understanding of our closest star. 1.4 Nuclear Reactions and Neutrinos The main source of the Sun’s power is the fusion of hydrogen into helium in its core. The reaction can be summarized as 4 1 H→ 4 He+ 2 e + + 2 ν e . (1.2) The two positrons balance the change in charge from two of the protons converting into neutrons. The two neutrinos (leptons) balance out the positron (anti-lepton) creation in order to maintain lepton number conservation. By comparing the masses of the two sides of the reaction, we find that 26.73 Mev of energy is released each time four 6 hydrogen nuclei are converted into a helium nucleus. About 99% of this energy goes toγ-rays and the kinetic energy of the particles in the Sun (the remaining 1% is carried off by neutrinos which will be discussed at the end of this section). The simultaneous approach of four protons close enough to engage in nuclear fu- sion is so unlikely that the reaction described in equation (1.2) rarely occurs. Instead, helium is formed from hydrogen through several different chain reactions. The ba- sic reactions that take place in the Sun involve the direct fusion of protons, and are therefore called PP-chains. The basic PP-chain is PP-I: 1 H+ 1 H→ e + +ν e + 2 D 2 D+ 1 H→ γ+ 3 He 3 He+ 3 He→ 2 1 H+ 4 He (1.3) Once the 3 He is created in the second step of the PP-I chain, it can instead branch off to produce 4 He via the PP-II chain: 3 He+ 4 He→ γ+ 7 Be 7 Be+ e − → ν e + 7 Li 7 Li+ 1 H→ 4 He+ 4 He (1.4) A further branching of the PP-chains can occur once 7 Be is produced, leading to the PP-III chain: 7 Be+ 1 H→ γ+ 8 B 8 B→ e + +ν e + 8 Be 8 Be→ 4 He+ 4 He (1.5) Since we cannot go to the center of the Sun to watch these reactions take place, we 7 rely on neutrino evidence and solar luminosity for insight into the processes that fuel the Sun. Since the Sun is in equilibrium and its power is gereated by fusion in the core, the energy leaving the Sun (its luminosity) should equal the energy generated by the nuclear reactions. This analysis leads to an approximation of nuclear reaction rates. A check on these nuclear reaction rates can be obtained from solar neutrinos. The maximum amount of energy that can be generated by the above reactions is 26.73 MeV . Both muons and tauons are far too massive to be created in these reactions. Thus, due toμ- andτ- lepton number conservation, neitherμ norτ neutrinos are produced. Since positrons are created by these reactions, electron neutrinos are the only flavor produced in the Sun. However, neutrinos oscillate between all three flavors, so at Earth we detect only 1/3 of the solar neutrinos asν e , 1/3 asν τ , and 1/3 asν μ . Counts of solar neutrinos measured at the Earth confirm the nuclear reaction rates determined via solar luminosity. However, neutrinos are still not fully understood. The transitions between different flavors depend on the mixing angle and neutrino mass differences, which have yet to be determined. Thus there is room for improvement in calculating solar nuclear reaction rates. This prompted molecular dynamics studies of particle interactions in the solar core [75, 76, 77, 78]. Our analysis of the techniques used in these studies is described in part I of this thesis. 8 Part I Molecular Dynamics Chapter 2 Electrostatic Screening Under the extreme temperatures and densities of the solar core, the plasma is fully ionized. The free electrons and ions interact via the Coulomb potential U(r i j )= Z i Z j r i j . (2.1) In this Coulomb system, nearby plasma is polarized by each ion (see exaggerated illustration in Fig. 2.1). When two ions approach with the possibility of engaging in a nuclear reaction, each ion is surrounded by a screening cloud. Thus each ion is Figure 2.1: Plasma polarized by ion. 9 Figure 2.2: Approaching ions are each surrounded by their own screening cloud. Figure 2.3: Comparison between the Coulomb potential and the screened potential. Particles with a given relative energy can approach closer when subject to the screened potential than the Coulomb potential. attracted by the electrons and repelled by the protons in its partner’s cloud (Fig. 2.2). The combined effect of the particles in the screening clouds on the potential energy between the two ions is referred to as the screening effect. This screening effect reduces the standard Coulomb potential between approaching ions in a plasma to a screened potential which includes the contribution to the potential from the surrounding plasma. This reduced potential allows the ions to get closer (Fig. 2.3), enabling them to tunnel through the potential barrier more easily, thus enhancing fusion rates. 10 2.1 Salpeter Approximation for Electrostatic Screening In his 1954 paper, Salpeter derived the screening potential for electrons and ions in a plasma under the assumption of weak screening [71] Z i Z j r i j << k B T . (2.2) He began by assuming that the total potential energy is of the form U total (r)= Z i Z j e 2 φ total (r) (2.3) with contributions to the total potential from both the interacting ions and the plasma. φ total (r)= 1 r +φ screening (r) ! (2.4) The charge density and electrostatic potential are related via Poisson’s equation ∇ 2 (Z i eφ total (r))=−4π¯ ρ(r)− 4πZ i eδ 3 (r) . (2.5) The density for the electrons and ions in the plasma is the field-free density times the Boltzmann factor ¯ ρ(r)= ρN 0 ze A " exp −Z i ze 2 φ total (r) k B T ! − exp Z i e 2 φ total (r) k B T !# (2.6) where A and z are the atomic weight and charge of the main component of the gas and N 0 is Avagadro’s number. Combining equations 2.4, 2.5, and 2.6 and expanding the 11 exponentials yields a linear approximation of the Poisson-Boltzmann equation ∇ 2 φ(r)= 4πρN 0 z 2 + z A ! e 2 k B T " 1 r +φ(r) # . (2.7) The solution is φ(r)= 1 r exp −r R − 1 (2.8) where R 2 = k B T 4πρN 0 e 2 ! A z 2 + z . (2.9) This solution is equivalent to the formalism of the Debye-H¨ uckel theory [18], so the relevant potential is Yukawa φ(r)= z r e −αr (2.10) whereα is the inverse of the Debye length. 2.2 Deviations from Salpeter Screening Salpeter’s screening equation has been widely accepted and applied in solar equations of state and nuclear reaction rate calculations. However, the problem of “missing neutrinos” and, more recently, the quest to determine neutrino mass differences have inspired a second look at the question of electrostatic screening. A better approxima- tion for screening would lead to adjustments in solar nuclear reaction rate calculations which would alter the number of neutrinos that are expected to be produced in reac- tions in the solar core. Additionally, experimental and numerical results for screening (Assenbaum et al [9], Engstler et al [24], Shoppa et al [79], Rolfs [70], Mochkovitch and Hernanz [55], Ichimaru et al [37], and Arafune and Fukugita [7]) have challenged the applicability of Salpeter’s formalism. 12 To explain the experimental deviations, several authors derived alternate formu- lae for screening (Carraro et al [15], Opher and Opher [60], Shaviv and Shaviv [75], Savchenko [73], Lavagno and Quarati [45], and Tsytovich [83]). Their derivations were refuted in subsequent papers. Bahcall, Brown, Gruzinov, and Sawyer [10] sum- marized arguments in Salpeter’s defense. They also pointed out that all of the non- Salpeter screening corrections are different, claiming that “there is only one right an- swer, but there are many wrong answers.” 2.3 Dynamic Effects Although Salpeter’s formalism remains the preferred treatment for electrostatic screen- ing, the question of dynamic screening remains open. Since the protons in the plasma are much slower than the electrons, they are not able to rearrange themselves as quickly around fast moving ions. Nuclear reactions require greater energies than the average thermal energy in the solar core (see table 2.3 for the energies of relevant reactions in the Sun [75]), thus the ions that are able to engage in nuclear reactions are in the high end tail of the Maxwell-Boltzmann velocity distribution (Fig.2.4). The relevant ions are much faster than the thermal ions and would therefore not be accompanied by their full screening cloud. This dynamic effect would modify the screened potential for the interacting ions from the simplified Yukawa form and alter the nuclear reaction rates from current values. 2.4 Shaviv and Shaviv In order to determine the effects of dynamic screening, Shaviv and Shaviv [75], [76], [77], [78] have developed numerical simulations using the techniques of molecular 13 Reaction E Gamov /kT Δ/kT p+ p 4.8 4.93 p+ Be 7 11.7 8.50 He 3 + He 3 14.0 9.40 He 3 + He 4 14.6 9.62 p+ N 14 17.3 10.47 Table 2.1: Relevant Reactions: Gamov energy and the width of the Gamov peak with respect to the thermal energy. Figure 2.4: Maxwell-Boltzmann velocity distribution for protons in the solar core (ve- locity is in atomic units). 14 dynamics to follow the time evolution of a system under solar conditions. By directly calculating the motion of ions and electrons due to Coulomb interactions, they are able to avoid the mean field assumptions inherent in Salpeter’s derivation. The results of their simulations differ from those predicted by Salpeter’s theory. The pair correlation functions of their system do not match those computed using Salpeter’s formalism [75]. They have also demonstrated that the energy transferred to a pair of ions by the plasma depends on the relative kinetic energy of the pair [78]. Shaviv and Shaviv attribute these differences to dynamic effects. 15 Chapter 3 Molecular Dynamics Molecular dynamics is a class of simulation techniques that can be used to numerically solve the dynamics of a classical many-body system. Although the plasma of the solar core is fully ionized and has no molecules, these are general techniques that can be applied to the two-component plasma of electrons and ions. This chapter describes the tools from molecular dynamics [6], [27], [57], [84] that are used in our simulations. 3.1 Velocity-Verlet Algorithm 3.1.1 Verlet Algorithm Molecular dynamics simulations are deterministic, using Newton’s equation of motion m i d 2 ~ r i dt 2 = ~ f i (3.1) 16 to compute the dynamics of each particle in the system at each time step. The relevant potential for our system is the Coulomb potential U(r)= q i q j r i j (3.2) which depends only on the positions of all the particles~ r N . Therefore, given a set of initial positions and velocities for all particles in the system{~ r i (0),~ v i (0)|i = 1, N}, the positions and velocities at a later time t can be obtained by integrating the following differential equation: m i ~ .. r i =− ∂U(~ r N ) ∂~ r i . (3.3) In order to solve this problem numerically, time must be discretized into small steps Δt, so that for each time step, the problem becomes: given~ r N (t),~ v N (t), compute~ r N (t+ Δt),~ v N (t+Δt). The Verlet discretization scheme for computing~ r N (t+Δt),~ v N (t+Δt) is derived via Taylor expansions of~ r(t+Δt) and~ r(t−Δt) ~ r(t+Δt)=~ r(t)+~ v(t)Δt+~ a(t) (Δt) 2 2 + ~ ... r(t) (Δt) 3 6 + O(Δt 4 ) (3.4) ~ r(t−Δt)=~ r(t)−~ v(t)Δt+~ a(t) (Δt) 2 2 − ~ ... r(t) (Δt) 3 6 + O(Δt 4 ). (3.5) Equations 3.4 and 3.5 can be added to give ~ r(t+Δt)+~ r(t−Δt)= 2~ r(t)+~ a(t)Δt 2 + O(Δt 4 ), (3.6) so that the position of each particle is ~ r i (t+Δt)= 2~ r i (t)−~ r i (t−Δt)+~ a i (t)Δt 2 + O(Δt 4 ). (3.7) 17 The velocities can then be obtained by subtracting equations 3.4 and 3.5 to give ~ r(t+Δt)−~ r(t−Δt)= 2~ v(t)Δt+ O(Δt 3 ), (3.8) so the velocity of each particle is ~ v i (t)= ~ r i (t+Δt)−~ r i (t−Δt) 2Δt + O(Δt 2 ). (3.9) Equations 3.7 and 3.9 for~ r i (t+Δt) and~ v i (t) can be used to step through the motion of the system. This version requires~ r i (t) and~ r i (t+Δt) to be calculated and~ r i (t−Δt) to be stored in order for~ v i (t) to be computed. 3.1.2 Velocity-Verlet Algorithm A more convenient form, the velocity-Verlet discretization scheme, can be derived from Taylor expansions of both the position and the velocity. The Taylor expansion of the position is shown in equation 3.4, while that of the velocity is ~ v(t+Δt)=~ v(t)+~ a(t)Δt+ O(Δt 2 ). (3.10) Better accuracy can be obtained by replacing the acceleration with the average [~ a(t)+ ~ a(t+Δt)]/2, yielding the velocity-Verlet discretization scheme: ~ r i (t+Δt)=~ r i (t)+~ v i (t)Δt+~ a i (t) Δt 2 2 (3.11a) ~ v i (t+Δt)=~ v i (t)+ [~ a i (t)+~ a i (t+Δt)] Δt 2 , (3.11b) 18 which is equivalent to the Verlet discretization scheme but enables the velocity and position to be computed at the same time step. By defining ~ v i t+ Δt 2 ! ≡~ v i (t)+~ a i (t) Δt 2 , (3.12) equations 3.11a and 3.11b can be rewritten as ~ r i (t+Δt)=~ r i (t)+~ v i t+ Δt 2 ! Δt (3.13a) ~ v i (t+Δt)=~ v i t+ Δt 2 ! +~ a i (t+Δt) Δt 2 . (3.13b) This scheme can be easily implemented with the following algorithm: Given initial conditions~ r i (t),~ v i (t), 1. compute~ a i (t) for i= 1, N from{~ r i (t)} 2.~ v i (t+ Δt 2 )←~ v i (t)+ Δt 2 ~ a i (t) for i= 1, N 3.~ r i (t+Δt)←~ r i (t)+~ v i (t+ Δt 2 ) for i= 1, N 4. compute~ a i (t+Δt) for i= 1, N from{~ r i (t+Δt)} 5.~ v i (t+Δt)←~ v i (t+ Δt 2 )+ Δt 2 ~ a i (t+Δt) for i= 1, N 6. go to 2 This velocity-Verlet algorithm is used to calculate the trajectories{~ r i (t),~ v i (t)} of the system. 3.1.3 Energy Conservation If Newton’s equation (3.1) is integrated exactly, the total energy is conserved. Defining the total energy E as the sum of the kinetic energy K and the potential energy V, the 19 energy can be writen as E = K+ V = N X i=1 m i 2 . r i 2 + X i< j U(r i j ). (3.14) To prove that energy is conserved, its derivative is shown to be zero: . E = X m i 2 2| . r i |·| .. r i |+ X . r i · ∂U ∂r i ! (3.15a) = X m i . r i · .. r i + X . r i · ∂U ∂r i (3.15b) = X . r i · m i .. r i + ∂U ∂r i ! (3.15c) = X . r i · (m i .. r i − F i ) (3.15d) where m i .. r i − F i = 0 (3.16) by Newton’s equation (3.1). Hence the total energy of the system is conserved by exact integration. In practice, however, time is discretized for numerical integration. With this ap- proximation, the total energy is no longer rigorously conserved. Smaller time steps better approximate exact integration and energy conservation. The time step should therefore be chosen to be small enough to adequately conserve the total energy yet not so small that the system evolves too slowly. 20 Figure 3.1: Periodic boundary conditions: the simulation box (bold) is surrounded by an infinite periodic lattice of identical boxes. When a particle leaves one side of the simulation box, its image from a neighboring box enters on the opposite side. 3.2 Periodic Boundary Conditions and Minimum Im- age Convention The system used in our molecular dynamics simulations is a box filled with electrons and protons. Since the number of particles is not very large (N = 1000 for our prelim- inary tests), many of the particles are at or near the surface of the box. The purpose of the model is to simulate properties of a macroscopic plasma. The environment of the surface atoms in the simulation differs from that of the bulk plasma, so surface effects must be eliminated to accurately model the conditions of the solar plasma. Surface effects are easily removed by implementing periodic boundary conditions. With this technique, the original box is surrounded by an infinite periodic lattice of identical boxes (see Fig.3.1). The particles in the box interact with the entire infinite system, allowing the image particles in the surrounding boxes to mimic the effect of being surrounded by an infinite homogeneous plasma. 21 Since all of the infinite boxes are identical, only the positions and velocities of particles in the original box are needed. The motion of the image particles in the surrounding boxes imitates the motion in the original box. When a particle leaves the original box, it’s image leaves a neighboring box and enters the original box from the other side, as seen in Fig.3.1. This condition is enforced in each dimension by adjusting the coordinate x i of each exiting particle by the box length L x so that the departing particle becomes its image particle that just entered the box on the other side: x i ← x i − sign L x 2 , x i ! − sign L x 2 , L x ! , (3.17) where the function sign(A, B)=                A if B > 0 −A if B≤ 0 . (3.18) Although the particles theoretically interact with all of the infinite images of every particle, in practice the original box should be made large enough that contributions to the potential from beyond L x 2 are negligible. This enables us to consider the potential on each particle from only the nearest image of the other particles. Since the location of the border of the original box is arbitrary, this minimum image convention amounts to shifting the boundary of the box so that the ith particle is at the center (note that the volume and number of particles within the box remain fixed). Then all of the ith particle’s nearest images are contained within the shifted central box. The box is re- centered around the ith particle by adjusting the distance between it and each of the other particles in every dimension. x i j ← x i j − sign L x 2 , x ij + L x 2 ! − sign L x 2 , x ij − L x 2 ! . (3.19) 22 This is done for each particle at every time step in the simulation so that it only interacts with the nearest image of each particle. Other boundary conditions such as reflective or free boundaries can also be used in molecular dynamics simulations. However, our system is small (N = 1000 for the test system), so a reflective boundary would introduce unwanted surface effects. A free boundary would also be inappropriate because it would fail to preserve the volume and density of the system. 3.3 Smooth and Continuous Potential The minimum image convention described in the previous section introduces a dis- continuity in the potential at the edge of the box. Since only the closest image of each particle is considered, all interactions are cut off beyond a distance of L x 2 . In order to derive the force from the potential, the potential must be smooth and continuously differentiable. The discontinuity can be repaired by shifting the potential so that both U(r) and dU(r) dr are continuous at the cutoff r cut U 0 (r i j )= U(r i j )− U(r cut )− (r i j − r cut ) dU dr r=r cut (3.20) for r i j < r cut and zero elsewhere. The smooth and continuous version of the Coulomb potential is thus U 0 (r i j )=                Z i Z j 1 r i j + r i j r 2 cut − 2 r cut if r i j < r cut 0 if r i j ≥ r cut . (3.21) 23 Figure 3.2: The Lennard-Jones potential. This form of the potential eliminates the discontinuity created by truncating the poten- tial. Analysis of an appropriate cutoff radius for the Coulomb system is discussed in chapter 6. 3.4 Dealing with Long-Range Forces The long-range nature of the Coulomb potential introduces additional challenges to a molecular dynamics simulation. For many of the materials commonly modeled using molecular dynamics, short-range potentials such as the Lennard-Jones potential of a simple fluid u(r)= 4 " σ r 12 − σ r 6 # (3.22) are appropriate. This potential decreases drastically with the separation of the two particles (see Fig.3.2). It can therefore be truncated quite reasonably at a distance of 24 Figure 3.3: The Coulomb potential. r cut = (2.5− 3)σ. The Coulomb potential, on the other hand, decays slowly with particle separation (see Fig.3.3). Particles at much greater distances must be included to accurately compute the force on each particle. Ideally, the contribution of each of the infinite image particles should be included in the potential sums and the potential should never be truncated. Unfortunately this would take infinite computing time and is therefore impractical. Fortunately, more efficient techniques have been developed to compute long-range interactions [27] [6]. In the following sections, several of these techniques are described and evaluated for appropriateness in our simulations. 3.4.1 Tail Corrections When short-ranged potentials are truncated as described in section 3.3, the systematic error in the potential energy of the system can be corrected by adding a tail contribution 25 to the total potential energy U tot = X i< j u c (r i j )+ N n 2 Z ∞ r c u(r)4πr 2 dr, (3.23) where u c is the truncated potential energy and n is the average number density. The tail correction is only finite if the potential energy U(r) decays faster than r −3 . Since this is not the case for the Coulomb potential, truncations with tail corrections should not be used in Coulomb systems. 3.4.2 Ewald Sums Ewald summation [25, 27] is the most widely used technique for computing long- range forces. This method relies on the periodic nature of the infinite lattice to sum the interaction between a particle and all of the infinite periodic images. The summation of the potential energy over the infinite lattice can be written as U = 1 2 X → n 0 N X i=1 N X j=1 Z i Z j | → r i j + → n| , (3.24) where the sum over → n is the sum over n x L, n y L, n z L (for a three dimensional system of side length L), where n x , n y , and n z are integers and the prime indicates that i = j is omitted for → n = 0 (ie: the particle interacts with all of its periodic images but not with itself). The sum in equation 3.24 is conditionally convergent, depending on the order in which the terms are added. The Ewald method improves convergence of the sum by rewriting the charge den- sity of the system. Here, the point charges are represented by delta functions. Since the system is electrically neutral, each particle is surrounded by a charge distribution 26 Figure 3.4: Representation of the Ewald sum method: each point-particle’sδ-function is represented as (a) aδ-function minus a gaussian (b) plus a gaussian. of equal and opposite net charge. The Ewald method represents this effect by adding a Gaussian charge distribution centered on each particle with opposite sign to the parti- cle. ρ i ( → r )= Z i α π 3/2 exp(−αr 2 ) (3.25) where α is an arbitrary parameter describing the width of the distribution. The added cloud effectively neutralizes the long-range nature of the interaction by screening neighboring particles. The system, however, consists of point charges not screened charges, so the screening clouds that were added to each particle are later removed by adding a second set of Gaussian charges centered on each particle, this time with the same charge as the particle, as shown in Fig.3.4. The combination of the original point charges and the first set of Gaussian charges that screen them is now short ranged and can be summed in real space. The second set of Gaussian charges added to cancel the 27 screening clouds can then be summed as a Fourier series in reciprocal space. Finally, to remove the self-interaction of each charge cloud with itself, the self-energy correction (which is just a constant) is subtracted. The complete Ewald sum becomes U coul = 1 2L 3 X → k,0 4π k 2 |ρ( → k)| 2 exp( −k 2 4α )− α π 1/2 N X i=1 q 2 i + 1 2 N X i, j q i q j erfc( √ αr ij ), (3.26) where → k is the reciprocal vector 2π → n L 2 and erfc is the complimentary error function. Lattice methods like the Ewald sum rely heavily on the periodic nature imposed by the model. This is problematic for our purposes, because the relevant particles for nu- clear reactions are the few at the high-end tail of the velocity distribution. Shaviv and Shaviv artificially increase the number of fast particles in their simulations in order to attain better statistics for fast particle interactions, disrupting the velocity distribution of their system, an effect that is exaggerated by the infinitely repeating lattice. Assum- ing that each of the fast particles has an infinite number of periodically spaced images is wrong. Even without the addition of fast particles, using periodic images replicates fluctuations in the charge distribution instead of averaging them out as a truly large system would. If the number of particles is increased to reduce these problems, the computational effort for Ewald sums which scales as order N 3/2 becomes expensive for large systems, with the sum in reciprocal space scaling as order N 2 . Since lat- tice methods overemphasize the periodic nature of the system, Ewald sums are not appropriate for use in our simulations. 3.4.3 Particle-Mesh Method The particle-mesh method [34, 27] separates the net force on a particle into short- and long- range contributions. The short-range interaction is summed directly. The long- 28 range interactions are then computed with the particle-mesh method. This method divides the simulation box into a finely-spaced mesh, taking advantage of the fact that it is more efficient to solve the Poisson equation when the charges are distributed on a grid. The charge density of the system is approximated by assigning charges or fractions of charges to the grid points. The Poisson equation is then solved for the potential due to the approximated charge distribution on the grid (via fast Fourier transform), giving the potential at each mesh point. The potential at each mesh point is differentiated numerically so that the forces from the long-range interactions can be assigned to each particle. For large N, this method scales as order N, making it useful for large systems. The particle-mesh method needs too many mesh points per particle to give good accuracy when particles are close together. This makes it computationally expensive for collisional systems, and therefore not ideal for use in our system. 3.4.4 Effective Cutoff Radius Instead of handling the long-range force with the techniques described in the previous sections, Shaviv and Shaviv implement a simple cutoff radius [78]. They only include electrons and protons within a given distance from each particle in the sum for that particle’s potential. No tail corrections are included, though a thermostat is used to prevent the discontinuities in the potential from disrupting the temperature of the sys- tem. This technique is based on the assumption that if a large enough cutoff radius is used more distant interactions will be screened by the surrounding plasma. Shaviv and Shaviv limit their calculations to particles within a few Debye radii, including 200- 420 interactions for each particle. They claim that this yields acceptable results. The validity of this assumption is examined in chapter 6. 29 (a) classical system (b) quantum system Figure 3.5: Without quantum corrections, the energy jumps dramatically when parti- cles get too close. 3.5 Quantum Effects In order to test the assumptions and approximations of Shaviv and Shaviv, our initial test systems were stripped-down classical versions of their simulations. One approxi- mation we attempted to avoid was the use of effective potentials to account for quan- tum effects. Without including quantum effects, our system was not able to handle the close approach of oppositely charged particles. Whenever a proton and electron got too close, the system’s energy jumped. Fig.3.5 shows 3D systems of 1000 particles with T = 47.5 E H , number density n 0 = 6× 10 24 cm −3 , and identical thermostats. This is indicative of the fact that when two particles in a two-component plasma approach each other closer than the thermal deBroglie wavelength Λ α = h √ 2πm α k B T , (3.27) the effects of quantum diffraction and symmetry become significant. The thermal de- Broglie wavelengths for our system are shown in table 3.1. The average interparticle 30 particle mass (m e ) Λ α (a 0 ) proton 1836 .008488 electron 1 .3637 Table 3.1: Thermal deBroglie wavelenghts for protons and electrons at a temperature of 47.5 E H . distance in our model is 1 a o . SinceΛ p << 1, a classical treatment is generally suf- ficient for the protons. However, it will be common to have two electrons that are separated by a distance less thanΛ e , requiring quantum mechanical treatment. In order to fix this problem, we added the quantum approximations used by Shaviv and Shaviv [75], i.e. the effective pair potentials derived for a hydrogen plasma by Barker [11], Deutsch [20], and Deutsch, Gombert, and Minoo [21][22]. Quantum diffraction effects are described by v (d) α,β = e α e β r " 1− exp −r o α,β !# , (3.28) where o= ~ p 2πμ α,β k B T (3.29) and μ α,β = 1 1 m α + 1 m β (3.30) is the reduced mass of the pair. The exclusion principle between electrons is included by adding v (s) ee = k B T ln(2) exp −r 2 πo 2 ee ln(2) ! . (3.31) These potentials only differ from the classical Coulomb potential at short distances, where quantum effects are relevant. Once the effective pair potentials were included, 31 the simulations ran smoothly without disruptions to the energy conservation, as seen in Fig3.5(b). 3.6 Thermalization and Time Averages After a system is initialized (see chapter 4), the integration scheme must run for a few times the amount of time it takes for a thermal proton to move the average interparticle distance in order to remove any artifacts of the initial conditions. As the particles inter- act, they exchange energy, equilibrating the system. Once the properties of the system no longer change with time (see section 5.2 for thermalization checks), the system is considered thermalized and measurements can be taken. Since the measurements are subject to statistical noise, properties are averaged over many time steps to improve accuracy. 32 Chapter 4 Setting Up Our System 4.1 Units Choosing appropriate units can simplify simulations and increase the speed by reduc- ing the number of calculations necessary. Our system treats protons and electrons separated by distances on the order of the Bohr radius and includes quantum effects. Atomic units work well in this case, since the mass and charge of the electron are both equal to one, as are the Bohr radius a 0 and the reduced Planck constant~ (see appendix for complete description of atomic units). In these units, the Coulomb potential energy is U(r)= q i q j r i j (4.1) where q i and q j are plus or minus one. The temperature of the solar core, expressed as k B T in Hartrees, is 47.5 E H . The density is 0.4− 1.0 m e /a 3 0 . These numbers are more manageable for our simulations than the 15× 10 6 K and 150 g/cc generally used in solar tables. 33 4.2 Initializing Coordinates and Velocities The initial positions and velocities should be chosen to be compatible with the struc- ture that the system is designed to simulate. Generally, the particles are assigned initial positions on a lattice and allowed to relax to a more realistic distribution during ther- malization (initial runs which are thrown away before property calculation can begin). Since the relative velocities of electrons and protons scale with the square root of their masses p m p /m e ≈ 43, our two-component plasma is plagued by the different time scales of electron and proton motions. The electrons equilibrate quickly while the protons drag behind. To reduce the amount of time needed for the full system to reach thermal equilibrium, we start with initial conditions that more closely resemble the plasma structure instead of beginning with a lattice. The initial coordinates are as- signed randomly in each dimension as the box length times a random number between 0 and 1. The initial proton and electron velocities are assigned via Maxwell-Boltzmann velocity distributions f (v)d 3 r d 3 v= N V m 2πk B T ! 3/2 v 2 exp −mv 2 2k B T ! d 3 r d 3 v (4.2) with the center of mass velocity set equal to zero. The initial structure and any artifacts of the initialization should disappear during equilibration. 4.3 Assigning Charge and Mass During the set up, each particle is determined to be either a proton or an electron. To ensure that the protons and electrons are randomly distributed, each particle is assigned a random number between 0 and 1. If the number is less than 0.5, the particle is given 34 the charge and mass of an electron. Otherwise it is given the charge and mass of a proton. To maintain the necessary overall charge neutrality, the numbers of protons and electrons are counted during assignment. If the numbers do not match, the charges and masses are reassigned. The end result is a plasma with an even number of randomly distributed protons and electrons. 4.4 Temperature and Density In order to simulate the plasma in the core of the Sun, the system must model the temperature and density of that region. This is done by scaling the velocities to repre- sent the desired temperature and scaling the particle separation to recreate the desired density. The temperature is related to the kinetic energy of a system by the equipartition theorem. Each degree of freedom N f in an ideal gas contributes 1 2 k B T to the average kinetic energy, so that the average kinetic energy of the system is given by < K >= * N X i=1 p 2 i 2m i + = N f 2 k B T, (4.3) where k B is the Boltzmann constant, and T is the temperature of the system. This makes the average velocity of each species ¯ v= r N f k B T m . (4.4) (As indicated in section 4.2, the average electron velocity will be about 43 times faster than the average proton velocity.) For a system in statistical equilibrium, the velocities 35 (a) protons (b) electrons Figure 4.1: Maxwell-Boltzmann velocity distribution for a temperature of 47.5 E H for (a) protons and (b) electrons. of the particles are distributed according to the Maxwell-Boltzmann distribution f (v)d 3 r d 3 v= N V m 2πk B T ! 3/2 v 2 exp −mv 2 2k B T ! d 3 r d 3 v. (4.5) See Fig.4.1 for the Maxwell-Boltzmann velocity distributions for protons and elec- trons at a temperature of 15 million K. To reproduce the desired density, the average interparticle spacing should be ¯ r= 1 3 √ n (4.6) where n is the number density. The box length is scaled by the average interparti- cle distance to ensure that the system has the correct size. Assigning the appropriate speeds and separations creates a system that simulates the temperature and density of the plasma in the core of the Sun. 36 Chapter 5 Our Preliminary Tests We conducted preliminary tests to evaluate the accuracy of the numerical simulations [47, 56]. The initial test system consisted of 1000 protons and electrons in a cube of side length 10 Bohr radii at a temperature of 15 million K, subject to periodic boundary conditions. Beginning with purely classical calculations, the particles interact via the Coulomb potential. Molecular-dynamics routines compute the time evolution of the positions and velocities of the particles in the system. Properties of the system can be calculated or the motion of individual particles can be followed. With this test system, we check the simulations before applying them to analyze Shaviv and Shaviv’s initial hypotheses. 5.1 Conservation of Energy and Velocity Distributions In order to ensure that the simulations are working properly, we first test basic prop- erties such as the conservation of energy and the velocity distributions of the protons and electrons. As seen in Fig.5.1, the energy of our system remains constant over time. Fig.5.1 shows the total energy computed during a sample time period in which a 37 Figure 5.1: Conservation of energy for the pilot MD system. thermal proton would move ten times the average inter-particle distance. The velocities of particles in an ideal gas in thermal equilibrium are given by the Maxwell-Boltzmann distribution f (v)d 3 r d 3 v= N V m 2πk B T ! 3/2 v 2 exp −mv 2 2k B T ! d 3 r d 3 v. (5.1) As seen in Fig.5.2 and Fig.5.3, the speed distribution of electrons and protons in our relaxed system matches the expected Maxwell-Boltzmann distribution. These tests show our model’s ability to reproduce the energy conservation and velocity distribution of a physical system with conditions similar to the solar core. 5.2 Pair Distribution Functions Pair distribution functions [6] provide a way to check the thermalization of the system. This function gives the probability of finding a particle within a distance r of another 38 Figure 5.2: Speed distribution for electrons compared to the Maxwell-Boltzmann dis- tribution. Figure 5.3: Speed distribution for protons compared to the Maxwell-Boltzmann distri- bution. 39 particle compared to the same probability for a system of the same density with a completely random distribution. The pair distribution function of a thermalized system is computed by integrating the normalized configurational distribution function over the positions of all but two of the particles g(r 1 , r 2 )= N(N− 1) ρ 2 Z Z dr 3 dr 4 ... dr N exp −U(r 1 , r 2 ,... r N ) k B T ! . (5.2) When computing the pair distribution function numerically, it is more convenient to use an equivalent definition which takes the ensemble average over pairs of particles g(r)= V N 2 * X i X j,i δ(r− r i j ) + (5.3) (here V is the volume). In practice, this is done by counting all pairs separated by a distance between r and r + dr and compiling a normalized histogram of all pair separations for each range of distances. For the two-component plasma, pair distribution functions can be computed sepa- rately for proton-proton, proton-electron, and electron-electron pairs. The relative ve- locities of electrons and protons scale with the square root of their masses p m p /m e ≈ 43. Thus the lighter electrons move much faster than the protons, engaging in more frequent close encounters through which energy is exchanged. The electrons therefore thermalize with each other more quickly than the protons. The proton-electron pair distribution functions indicate thermalization soon after that of the electrons, but the protons take much longer to reach a thermal distribution. This effect can be seen in the comparison of Fig.5.4. The numerically determined pair distribution functions can be compared to theo- 40 (a) early stages of thermalization (b) at a later time Figure 5.4: Pair correlations for the same system at different stages of thermalization. retical pair correlation functions derived from the Debye theory [44] w α,β = 1− Z α Z β k B T e −χr r (5.4) where χ= s 4π k B T X a n a Z 2 a . (5.5) Fig.5.5 shows the pair correlations of our system compared to those derived from De- bye theory. There are clear differences between the theoretical and numerical pair correlation functions. These differences may be caused by dynamic effects, as sug- gested by Shaviv and Shaviv [75] (although their pair correlation functions differ from both our MD results and the theoretical results). 41 Figure 5.5: Pair correlations of the MD system compared to Debye-based pair corre- lations (lines). 42 Chapter 6 Our Evaluation of Shaviv and Shaviv’s Techniques Shaviv and Shaviv apply the techniques of molecular dynamics to the Sun’s core in order to numerically determine the effect of screening. By directly calculating the motion of ions and electrons due to Coulomb interactions in the sun’s core, they are able to compute the effect of screening without the mean field assumption inherent in Salpeter’s approximation. However, Shaviv and Shaviv’s work includes a variety of other assumptions that need to be evaluated before their results can be accepted. This chapter starts with our pilot study of 1000 (half electrons, half protons) under conditions similar to the solar core (temperature T = 1.5× 10 7 K, number density n 0 = 6× 10 24 cm −3 ) to examine the assumptions and approximations from Shaviv and Shaviv’s work. We used modified systems to test the effect of various changes to the model. 43 (a) r cut = 1 (b) r cut = 2 (c) r cut = 3 Figure 6.1: Only particles within a distance less than the cutoff radius are included in the potential sum. 6.1 Effective Radius of Interaction A preliminary test is carried out to determine the effective radius of interaction in our system, beyond which the Coulomb potential of ions and electrons is screened out by the plasma. Shaviv and Shaviv limit their calculations to within a few Debye radii, in- cluding 200 - 420 interactions for each particle, which they claim is a sufficient number to yield acceptable results [78]. To test this hypothesis, we compute the average poten- tial energy per particle when various cutoff radii are imposed on the Coulomb sums. Only electrons and protons within a distance of the cutoff radius from each particle are included in the sum for that particle’s potential (see Fig. 6.1 for illustration). When a large enough cutoff radius is used, the potential of more distant particles will be screened by the plasma. Increasing the cutoff radius beyond this point would not result in a difference in the potential energy. Thus including particles beyond the effective radius would be irrelevant and computationally wasteful, while truncating the potential calculations at a shorter radius would introduce inaccuracies in the force calculations. In order to avoid discontinuities at the cutoff radius, we implement the smooth and 44 Figure 6.2: Average potential energy per particle computed using varying cutoff radii for the Coulomb sums. Here, the Debye length is 1.94 Bohr radii. continuous version of the Coulomb potential described in section 3.3 U 0 (r i j )= Z i Z j " 1 r i j + r i j r 2 cut − 2 r cut # (6.1) The positions and properties of each ion and electron are taken from a thermalized molecular dynamics simulation so that the particle distribution realistically represents the desired system. In a static calculation, the average potential energy per particle is computed using the Coulomb potential with various cutoff radii. Results of this cutoff radius test for are shown in Fig.6.2. Although the average potential energy per particle appears to level off within a few Debye radii, the potential energy does not become perfectly constant within the size of our box. Tests in larger systems are necessary to determine the appropriate cutoff radius for the Coulomb potential (see section 6.3). 45 6.2 System Size The cutoff radius test in the pilot study is inconclusive due to the size limitation of the system used. This highlights the need to develop tests for larger systems without increasing the compute time necessary for the system to reach equilibrium. Since the potential cannot be truncated at a distance of 5 a 0 or less, other time-saving methods must be introduced. Based on the desire to increase the cutoff radius without dramati- cally increasing the compute time, a two-dimensional system was developed. 6.2.1 Number of Dimensions To justify the two-dimensional approach, the effect of the number of dimensions on the system was checked. Since the relevant feature of the system for screening effect studies is the distribution of particles, pair correlation functions for a 2-D system with 2500 particles and a 3-D system with 1000 particles are compared. Both systems have the same temperature, density, and cutoff radius as the pilot study (T = 1.5× 10 7 K, n 0 = 6× 10 24 cm −3 , r cut = 5 a 0 ) and include quantum effects and a thermostat (as explained in the following sections). As seen in Fig.6.3, the similar pair correlations for the 2-D and 3-D systems indicate that 2-D systems can be used to save compute time without sacrificing accurate particle distributions. This enables not only further cutoff radius testing (see section 6.3) but also the development of a larger number of systems to test different parameters. 6.2.2 Number of Particles In order to check the effect of system size on the pair correlation functions, the 2-D system of 2500 particles described in section 6.2.1 is compared to a 2-D system of 100 46 Figure 6.3: Pair correlation functions for a two-dimensional system and a three- dimensional system. Figure 6.4: Pair correlation functions for 2-D systems with 100 and 2500 particles. particles. Other than the number of particles, all other system variables are identical. The comparison of the pair correlation functions in Fig.6.4 indicate that the number of particles in the system does not affect the distribution of the particles. It should be noted that both systems used in this study contain enough particles that the width of the simulation box is greater than or equal to twice the cutoff radius. In a system that is too small to meet this requirement, a new cutoff radius would be necessary to maintain spherical symmetry. In that case, the effect of this shorter cutoff radius would need to be tested. 47 Figure 6.5: Pair correlation functions for 2-D systems with 2500 particles and two different cutoff radii (5 and 7.07 Bohr radii). 6.3 Effective Radius of Interaction Revisited Now that a larger system has been developed in 2-D, the cutoff radius analysis can be continued. To check the effect of the cutoff radius on the distribution of particles, pair correla- tion functions are compared for two systems with two dimensions, 2500 particles, the same temperature and density as the pilot study (T = 1.5× 10 7 K, n 0 = 6× 10 24 cm −3 ), quantum effects, and a thermostat. One system has a cutoff radius of 5 a 0 ; the other has a cutoff radius of 7.07 a 0 . As seen in Fig.6.5, increasing the cutoff radius beyond that of the pilot study does not alter the particle distribution. This indicates that the cutoff radius of 5a 0 is sufficient for our system. 6.4 Quantum Effects As discussed in section 3.5, the simulations were initially restricted to classical poten- tials. However, classical mechanics could not appropriately handle the treatment of approaching particles with opposite charge. Whenever a proton and electron got too 48 Figure 6.6: The energy is unstable in a classical system. Figure 6.7: When quantum corrections are included, the energy is conserved. close, the energy of the system jumped irreparably (see Fig.6.6). When the potential was adjusted to account for quantum effects, energy conservation was restored (see Fig.6.7). This indicates the necessity of a quantum mechanical treatment for the hot, dense, two-component plasma. Comparison of the pair correlation functions for 3-D systems of 1000 particles with and without quantum effects are shown in Fig.6.8. The pair correlation for the classical system was calculated during a time in which the energy of the system was conserved and equal to that of the quantum system. The dramatic differences in the 49 (a) g pp and g ee (b) g pe Figure 6.8: Pair correlation functions for systems with and without quantum effects included in the potentials. pair correlation functions support the need for the inclusion of quantum effects in the simulations. This specifically illustrates that the quantum effects enable the protons to get closer to each other, increasing the nuclear reaction rates over the classical ap- proximation. When compared with the theoretical value derived from Debye theory (Fig.6.9), we see that neither the classical nor quantum simulations match the predicted result. This indicates that the system is affected by physics that is not addressed in the static screening theory. Shaviv and Shaviv include a quantum correction in their calculations. They do this via a widely used effective pair potential that approximates the effects of diffraction and symmetry (see section 3.5). Since direct quantum calculations are computation- ally expensive, such approximations of quantum effects are far more convenient. We elected to model quantum effects in the same way, however this effective pair potential should be tested and compared with alternative formulations before it is accepted for the simulations. 50 (a) g pp and g ee (b) g pe Figure 6.9: Pair correlation functions from Debye theory compared to those from molecular-dynamics systems with and without quantum effects included in the po- tentials. 6.5 Thermostat As described in section 3.4.4, Shaviv and Shaviv chose to deal with the long range Coulomb force by truncating the interactions. This truncation causes the system to heat up, preventing it from modeling a plasma of the desired temperature. In order to compensate for this problem, they implement a thermostatic force (Hoover [33] [35]) to maintain a constant temperature as measured through the instantaneous kinetic energy of the system. 6.5.1 Nos´ e-Hoover Thermostat Standard molecular dynamics techniques model a micro-canonical ensemble with con- stant number of particles N, volume V, and energy E. The Hamiltonian of the system can be altered to produce a new set of equations of motion that model the desired canonical system with constant N, V, and temperature T coupled to a heat bath with which it can exchange energy. Nos´ e [59] accomplished this by introducing a time- 51 scaling variable s, its conjugate momentum p s , and a parameter Q that acts as the effective mass of the thermostat. Nos´ e’s new Hamiltonian for a system with X degrees of freedom is H Nose =Φ(q)+ X p 2 2ms 2 + (X+ 1)k B Tln(s)+ p 2 s 2Q . (6.2) The equations of motion from this Hamiltonian are ˙ q= p ms 2 ˙ p= F(q) ˙ s= p s Q ˙ p s = X p 2 ms 3 − (X+ 1)k B T s . (6.3) Nos´ e showed that smooth, deterministic, time-reversible trajectories generate a phase- space distribution that is canonical in the variables q and p/s. Hoover [33] adjusted Nos´ e’s equations to produce the same effect using a thermo- dynamic friction coefficient instead of time scaling. Hoover’s formulation makes it possible to estimate finite-size effects on dynamical averages, which are of interest in Shaviv and Shaviv’s simulations. Hoover begins by writing equations 6.3 in a simpler form, reducing the time scale by s so that dt old = sdt new and expressing the rates as 52 derivatives with respect to t new ˙ q= p ms ˙ p= sF ˙ s= sp s Q ˙ p s = X p 2 ms 2 − (X+ 1)k B T . (6.4) The time-scaling variable s is then eliminated by rewriting the equations in terms of q and its first and second derivatives ¨ q= ˙ p ms − p ms ˙ s s = F m − ˙ qp s Q ≡ F(q) m −ζ ˙ q , (6.5) whereζ = p s /Q is the thermodynamic friction coefficient which evolves in time as ˙ ζ = P m˙ q 2 − (X+ 1)k B T Q . (6.6) The phase-space distribution from these equations can be made canonical without resorting to time-scaling by redefining p ≡ m˙ q and replacing Nos´ e’s X + 1 by X. Hoover’s resulting equations of motion for a canonical distribution are ˙ q= p m ˙ p= F(q)−ζ p ˙ ζ = P p 2 m − Xk B T Q . (6.7) 53 6.5.2 Thermostat Tests The Nos´ e-Hoover algorithm generates the correct canonical distribution if there is only one conserved quantity. Since the total energy of the extended system (including the heat bath) is conserved, conservation of momentum can only be included if the center of mass remains fixed. For a system in which there is more than one conservation law, Nos´ e-Hoover chains (Martyna, Klein, Tuckerman, [48]) must be used to produce the correct canonical distribution. A discussion of two different integration schemes for implementing Nos´ e-Hoover chains is presented in the following subsections. Nos´ e-Hoover Chains (NHC) The first integration scheme was developed by Martyna, Tuckerman, and Klein (MTK) following the formalism of Nos´ e and Hoover. MTK generated canonically distributed positions and momenta using a chain of Nos´ e-Hoover thermostats to drive a dynamical system (Martyna et al [48] [49]). The equations of motion they proposed for this implementation are ˙ q i = p i m i ˙ p i =− ∂V(q) ∂q i − p i p η 1 Q 1 ˙ η i = p η i Q i ˙ p η 1 =        N X i=1 p 2 i m i − Nk B T        − p η 1 p η 2 Q 2 ˙ p η j =        p 2 η j−1 Q j−1 − k B T        − p η j p η j+1 Q j+1 ˙ p η M =       p 2 η M−1 Q M−1 − k B T       . (6.8) 54 M is the number of chains of thermostats η j , with masses Q j and momenta p η j . N is the number of particles, N f is the number of degrees of freedom, and F i =−∇ i φ(r) is the force derived from the potential. The conserved quantity for these equations of motion is H 0 (p, q,η, p η )= V(q)+ N X i=1 p 2 i 2m i + M X i=1 p 2 η i 2Q i + Nk B Tη 1 + M X i=2 k B Tη i . (6.9) Trotter Expansions Energy conservation was not maintained in the system in which we implemented the MTK Nos´ e-Hoover chain method just described. Instead, the following alternate in- tegration scheme derived from classical propagators using operator factorization tech- niques (Martyna et al [50]) was applied. This method yields an explicit reversible integrator for the canonical ensemble. The time evolution of a system of coupled, first-order differential equations can be described by applying the evolution operator ~ Γ(t)= exp(iLt) ~ Γ(0) (6.10) iL= ~ ˙ Γ· ~ ∇ Γ , (6.11) where ~ Γ is a multidimensional vector of coordinates and velocities and iL is the Liou- ville operator. Generally, the evolution operator’s action on the coordinates cannot be computed analytically. Instead the system is evolved by P successive applications of a short-time approximation of the operator that is accurate at timeΔt= t/P ~ Γ(t)= P Y i=1        Y s exp(iL s Δt)        ~ Γ(0) . (6.12) 55 For an nth order factorization, this approximation has an overall error of t n /P n−1 . The Liouville operator for the equations of motion in the canonical ensemble is iL= N X i=1 ~ v i · ~ ∇ r i + N X i=1       ~ F i (~ r) m i       · ~ ∇ v i = N X i=1 v η 1 ~ v i · ~ ∇ v i + M X i=1 v η i ∂ ∂η i + M−1 X i=1 (G i − v η i v η i+1 ) ∂ ∂v η i + G M ∂ ∂v η M . (6.13) Here, a chain of M Nos´ e-Hoover thermostats is coupled to an N particle system with G 1 = 1 Q 1        N X i=1 m i v 2 i − N f k B T        (6.14) G i = 1 Q i (Q i−1 v 2 η i −1 − k B T) i > 1 . (6.15) From a generalization of the Trotter formula, the evolution operator can be rewritten as exp(iLΔt)= exp iL NHC Δt 2 ! exp iL 1 Δt 2 ! exp(iL 2 Δt) × exp iL 1 Δt 2 ! exp iL NHC Δt 2 ! + O(Δt 3 ) (6.16) with iL 1 = N X i=1       ~ F i (~ r) m i       · ~ ∇ v i (6.17) iL 2 = N X i=1 ~ v i · ~ ∇ r i , (6.18) and iL NHC = iL− iL 1 − iL 2 (6.19) 56 contains the rest of the iL terms. The Nos´ e-Hoover chain part of the evolution operator in a multiple time step approach (n c > 1) can be simplified as exp iL NHC Δt 2 ! = N C Y i=1 exp iL NHC Δt 2n c ! (6.20) with exp iL NHC Δt 2 ! = exp Δt 4n c G ∂ ∂vη M ! exp − Δt 8n c v η M v η M−1 ∂ ∂v η M−1 ! × exp Δt 4n c G M−1 ∂ ∂v η M−1 ! exp − Δt 8n c v η M v η M−1 ∂ ∂v η M−1 ! × exp        − Δt 2n c N X i=1 v η 1 ~ v i ·∇~ v i        exp        Δt 2n c M X i=1 v η i ∂ ∂η i        × exp − Δt 8n c v η M v η M−1 ∂ ∂v η M−1 ! exp Δt 4n c G M−1 ∂ ∂v η M−1 ! × exp − Δt 8n c v η M v η M−1 ∂ ∂v η M−1 ! exp Δt 4n c G M ∂ ∂v η M ! . (6.21) Typically (if the frequency of the chain is not high), n c can be taken to be one. This expression for the NHC part of the operator can be substantially reduced with a higher order algorithm such as exp iL NHC Δt 2 ! = N C Y i=1         n ys Y j=1 exp iL NHC w j Δt 2n c !         . (6.22) where the w j are chosen such that the error is order (Δt/2n c ) 5 . The values used for n ys and w j are n ys = 3 w 1 = w 3 = 1 2− 2 1/3 w 2 = 1− 2w 1 . (6.23) 57 Figure 6.10: Pair distribution functions for a micro-canonical system and a canonical system with a Trotter-expansion based thermostat. This apparently complicated approach is actually straightforward to implement and computationally inexpensive. First, the operator exp(iL NHC Δt/2) is applied to update the variables v, v η , andη. Then, the standard velocity-Verlet algorithm is applied to the updated velocities. The operator exp(iL NHC Δt/2) is applied again to the output of the velocity-Verlet algorithm. To ensure that the thermostat does not alter the distribution of the particles, pair correlation functions for a standard micro-canonical system and a canonical system with a Trotter thermostat are compared. Both are 2-D systems of 100 particles with quantum effects included and the same temperature, density, and cutoff radius as the pilot study (T = 1.5× 10 7 K, n 0 = 6× 10 24 cm −3 , r cut = 5 a 0 ). As seen in Fig.6.10, the thermostat does not alter the distribution of particles and is therefore acceptable for studying screening. 58 6.6 Simple Thermostat “It is well known that due to the truncation in the calculation of the force between the particles, the system heats up. This is particularly true for the light electrons. To overcome this problem, we add a thermostatic force (Hoover 1985, 1991) which forces the temperature of the plasma to be constant during the calculations. The relaxation from the initial condi- tions is therefore isothermal rather than adiabatic.” [75] This is the only information Shaviv and Shaviv provide about their use of a ther- mostat. They do not discuss how they chose to implement the thermostatic force. Fur- thermore, the only references they list on the matter are Hoover’s 1985 article [33] and subsequent thermodynamics text [35]. These references only present the mechanics of the thermostat, not its implementation. This leaves Shaviv and Shaviv’s implementa- tion of the thermostat open to question. Implementing the thermostat can be quite complex, as seen in section 6.5. The main challenge is that the thermostatic force is velocity dependent. The standard velocity-Verlet equations for a micro-canonical ensemble are derived from Taylor ex- pansions, resulting in expressions for the positions and velocities of the form ~ r i (t+Δt)=~ r i (t)+~ v i (t)Δt+~ a i (t) Δt 2 2 (6.24a) ~ v i (t+Δt)=~ v i (t)+~ a i (t)+~ a i (t+Δt) Δt 2 . (6.24b) If the same technique is applied to the Nos´ e-Hoover equations of motion, ~ r i (t+Δt)=~ r i (t)+~ v i (t)Δt+ [~ a i (t)−ξ(t)~ v i (t)] Δt 2 2 (6.25a) ~ v i (t+Δt)=~ v i (t)+ [~ a i (t)+~ a i (t+Δt)−ξ(t)~ v i (t)−ξ(t+Δt)~ v i (t+Δt)] Δt 2 . (6.25b) 59 The expression for the velocity is self-referential, with ~ v i (t+Δt) appearing on both sides of the equation. This prevents these equations from being integrated exactly [27]. Instead, the Nos´ e-Hoover thermostat is usually solved iteratively or implemented using a predictor-corrector scheme, making the solution no longer time reversible. This problem can be solved by using the MTK implementation as previously described. Since Shaviv and Shaviv do not address this issue, it is possible that they imple- mented the thermostat incorrectly. They most likely began by adapting the micro- canonical velocity-Verlet integrator by replacing the force derived from the potential with the total force including the thermostatic force. ~ r i (t+Δt)=~ r i (t)+~ v i (t)Δt+        ~ f i (t) m i −ξ(t)~ v i (t)        Δt 2 2 (6.26a) ~ v i (t+Δt)=~ v i (t)+ [~ a i (t)+~ a i (t+Δt)−ξ(t)~ v i (t)−ξ(t+Δt)~ v i (t+Δt)] Δt 2 . (6.26b) This naturally produces the same result as the Taylor-expansion derived equations 6.25a and 6.25b, however the implementation of this integration scheme is not clear. Since Shaviv and Shaviv do not reference MTK, Trotter-expansion based, or other time-reversible integrators, it is likely that they applied the thermostatic force in the simplest way, replacing the velocity-Verlet algorithm Given initial conditions~ r i (t),~ v i (t), 1.compute~ a i (t) for i= 1, N from~ r i (t) 2.~ v i (t+ dt 2 )←~ v i (t)+ dt 2 ~ a i (t) for i= 1, N 3.~ r i (t+ dt)←~ r i (t)+~ v i (t+ dt 2 ) for i= 1, N 4.compute~ a i (t+ dt) for i= 1, N from~ r i (t+ dt) 5.~ v i (t+ dt)←~ v i (t+ dt 2 )+ dt 2 ~ a i (t+ dt) for i= 1, N 6. go to 2 60 Figure 6.11: Pair correlation functions for a simple thermostat and for a Trotter ther- mostat. with a similar version in which the force is amended to include the thermostatic force Given initial conditions~ r i (t),~ v i (t), 1.compute~ a i (t) for i= 1, N from~ r i (t) 2.~ v i (t+ dt 2 )←~ v i (t)+ dt 2 [~ a i (t)−ξ(t)~ v i (t)] for i= 1, N 3.~ r i (t+ dt)←~ r i (t)+~ v i (t+ dt 2 ) for i= 1, N 4.compute~ a i (t+ dt) for i= 1, N from~ r i (t+ dt) 5.~ v i (t+ dt)←~ v i (t+ dt 2 )+ dt 2 [~ a i (t+ dt)−ξ(t)~ v i (t+ dt 2 ) for i= 1, N 6. updateξ(t+ dt) from KE(t+dt) 7. go to 2. To test whether an error of this nature is the cause of Shaviv and Shaviv’s molec- ular dynamics results, a system was developed using this simplified version of the thermostat. The pair correlation functions of this system were compared with a sys- tem generated from identical conditions but with the Trotter thermostat described in section 6.5.2. As seen in Fig.6.11, the particle distribution from the system with the simplified thermostat does not deviate from the time-reversible distribution. The sim- 61 ple thermostat also does not reproduce the pair correlation functions of Shaviv and Shaviv [75]. This only means that the differences in Shaviv and Shaviv’s pair correlation func- tions are not due to this type of thermostat implementation. It is still possible that they have used a simplified thermostat and that it has affected their other results. 6.7 Temperature While testing the effect of different simulation choices on the distribution of the par- ticles, the influence of physical parameters is also of interest. Here, the effect of the system’s temperature is examined. As seen in section 5.2, the classical expression for the pair correlation functions depends on the temperature and density of the system. To check the influence of tem- perature in our simulations, the pair correlation functions for a system with temperature T = 19.3 E H were compared with those of a system with temperature T = 47.5 E H (solar core temperature). Both systems have the same density, number of particles, number of dimensions, and type of thermostat (n 0 = 6×10 24 cm −3 , N = 100, ndim= 2, Trotter thermostat) and include quantum effects. As expected and shown in Fig.6.12, the lower temperature results in greater distances between particles of like charge and closer approaches for particles with opposite charge. 62 Figure 6.12: Comparison of pair correlation functions for systems of temperature 19.3 E H and 47.5 E H . 63 Part II Virial Expansions Chapter 7 Solar Equations of State Equations of state (hereafter EOS) are thermodynamic equations used to describe the relationships between state variables in matter. The simplest EOS is the ideal gas law p= Nk B T V (7.1) for gases in which the interaction between particles is negligibly weak. For non-ideal plasmas, as found in the Sun, more complex expressions are needed in order to include the effect of the interactions. Solar equations of state traditionally express either the pressure or the free energy throughout the Sun as a function of the state variables of the plasma at each depth, F(T, V, N) or p(T, V,μ). Once the EOS is obtained, it is straightforward to derive expressions for other thermodynamic variables (internal energy, entropy, adiabatic coefficients, etc.), thus obtaining a complete thermodynamic description of the solar interior. Over the last twenty years, two different versions of solar equations of state (MHD and OPAL, described in the following sections) were developed as part of the two most recent opacity recalculations. 64 7.1 The Mihalas-Hummer-D¨ appen EOS (MHD) Named for its creators, Mihalas, Hummer, and D¨ appen, the MHD EOS [36], [53], [17], [58], [82]) was developed as part of the international Opacity Project [74], [13]. This approach relies on the intuitive chemical picture, assuming that the building blocks in the plasma are nuclei, electrons, atoms, ions, and molecules and then writing an expression for the free energy in terms of the temperature, volume, and number of each species F total (T, V, N e , N p , N H , N He + , N He ++ ,... ). Since the chemical picture deals with numbers of species (N e , N p , N H ,...), it is best treated in the canonical ensemble (constant N, V, T). The canonical partition function is defined classically as Z c = X i e −βE i (7.2) or quantum mechanically as Z c = Tr e −βH , (7.3) where H is the Hamiltonian of the system, E i are energy eigenvalues of H, and β= 1 k B T . (7.4) The free energy can then be derived from the canonical partition function as F =−k B T ln Z c . (7.5) Assuming that F total exists as a function of T, V, N e , N p , N H ,... and that Z c can be writ- ten as a product of partition functions Z total c = Z e · Z p · Z H · ... · Z interaction , (7.6) 65 then F total can be written as the sum of simple parts F total = F ideal + F interaction (7.7) where F ideal (T, V, N e , N p , N H ,... )= F e (T, V, N e )+ F p (T, V, N p )+ F H (T, V, N H )+... (7.8) and F interaction = F DH (T, V, N e , N p , N H ,... )+ F con f ined atom (T, V, N e , N p , N H ,... ). (7.9) F DH is the well-known Debye-H¨ uckel free energy from screening theory F DH = −k B TV 12πr 3 D . (7.10) The main task of MHD is creating a specific model for the confined atom part of the free energy. For simplicity, the basic method of determining F con f ined atom is described for the case of a hydrogen plasma. In this case, F = F(T, V, N e , N p , N H ) . (7.11) There are two useful constraints on the system, conservation of mass N p + N H = constant (7.12) 66 and charge neutrality N p = N e . (7.13) The effective degree of freedom is the reaction parameter x, which describes how much hydrogen has been ionized. We can write N p = x, N e = x, N H = N 0 H − x . (7.14) The free energy can then be written as F total = F total (T, V, N 0 H , x) (7.15) where N 0 H is a constant. In order to determine the reaction parameter of the system (which depends on the temperature and volume), MHD uses the method of free energy minimization to find x such that F is minimum ∂F ∂x ! T,V = 0 . (7.16) The pressure is then just the derivative of the free energy p=− ∂F ∂V , (7.17) and other thermodynamic variables can be computed in a similar fashion. 67 7.2 OPAL The OPAL EOS (Rogers [67], Rogers et al [69], Rogers and Nayfonov [68]) was de- veloped as part of Livermore’s OPAL opacity recalculation project (Iglesias et al [38], Iglesias and Rogers [39], [40], [41], [42]. Instead of the chemical picture, this ap- proach relies on the physical picture. The physical picture does not require descriptions of atoms. It works with only elementary constituents (electrons, protons, and nuclei), recreating atoms through the formalism. The goal is to derive an expression for the pressure as a function of temperature, volume, and chemical potential p(T, V,μ). The physical picture allows the ensemble to be chosen by convenience. OPAL uses the grand canonical ensemble (constant T, V,μ). The grand canonical partition function is Z GC = ∞ X N=0 ∞ X i=1 e −β(E i −μN) . (7.18) For the case of an ideal gas, the pressure can be easily computed in the grand canonical ensemble p(T, V,μ)= k B T V ln Z GC (T, V,μ). (7.19) For a set of ideal particles, the sum over N in equation 7.18 can be rewritten as a sum of the occupations n i of the individual energy states i E N i = X i n i i . (7.20) When equation 7.18 is written as Z GC = ∞ X N=0 ∞ X i=1 e β(μ− i ) N , (7.21) 68 we are left with a geometric series, which reduces to Z GC = X i 1 1− e β(μ− i ) . (7.22) This leads to an expression for the pressure of an ideal gas p(T, V,μ)= k B T V        ∓ ∞ X i=1 ln h 1∓ e β(μ− i ) i        (7.23) (Note: upper sign for Bosons; lower sign for Fermions). For the plasma of the solar interior, OPAL systematically adds perturbations to the ideal case, computing the non- ideal terms from the grand canonical partition function. This leads to an expansion in which perturbation terms are grouped in clusters. OPAL replaces the chemical potential with the activity z j = (2s j + 1)λ −3 j e βμ j (7.24) which includes the effects of the chemical potential, temperature, and degeneracy g i = (2s j + 1), and where λ j = h √ 2πmk B T . (7.25) The resulting EOS expanded in terms of the activity p= k B T n z+ b 2 z 2 + b 3 z 3 +... o (7.26) is also called ACTEX (activity expansion). The first term represents the ideal gas solution, with the second-order term adding all effects of two-body interactions, the third-order term for three-body interactions, and so forth. 69 7.3 Planck-Larkin Partition Function One important feature that arises naturally from the physical picture but not the chem- ical picture is the Planck-Larkin partition function (PLPF) Z PL = X i g i e −β i − 1+β i . (7.27) This version of the internal partition function results from a regrouping of terms in the physical picture expansions to remove the detrimental divergence of the internal partition function Z int = X n g n e −βE n . (7.28) Since the chemical picture treats this divergence by appropriately weighting the terms of the internal partition function Z int = X n g n w n e −βE n , (7.29) the PLPF is not built into MHD. The success of OPAL over MHD in matching observational data is partially at- tributed to the presence of the PLPF. When MHD is modified to include the PLPF, its accuracy is dramatically improved (Aihua Liang, 2003). 7.4 Feynman-Kac While both MHD and OPAL agree with solar observations to three digits, OPAL proves to be more accurate in the fourth digit. However, precise solar data shows that there is still room for improvement in both models, meriting new EOS attempts. 70 The EOS described in the next chapter is based on the Feynman-Kac representa- tion. Like OPAL, it is a physical-picture EOS evaluated in the framework of the grand canonical ensemble and naturally contains the crucial PLPF. However this method is based on virial density expansions evaluated by the path-integral formalism of Feyn- man and Kac [26, 80]. The Feynman-Kac representation equates a quantum point- charge system to a classical system of extended filaments. This allows all the non-ideal contributions to the EOS to be calculated systematically, exactly, and analytically. The downfall of the virial density expansion method is the inherent limitation on the do- main of validity. The Feynman-Kac approach (FK) can only be used under temperature and density conditions close to full ionization. In spite of this limitation, FK can be applied to large parts of the solar interior, as shown in chapter 9. 71 Chapter 8 Path-integral Formalism in the Feynman-Kac Representation In the Sun, matter can be described as a quantum plasma of electrons and nuclei. This chapter presents an equation of state formalism for such quantum plasmas based on the formulation of the Feynman-Kac (FK) representation [4, 2, 3]. The FK formalism equates a point-charge quantum system to a classical system of extended filaments, enabling the quantum Coulomb plasma to be described with classical techniques. The path-integral formalism is applied in two distinct steps. First, only Maxwell- Boltzmann statistics are considered (all exchange effects are omitted in this treatment and will be dealt with in step two). The FK representation is used to introduce an equivalent classical system of closed filaments interacting via two-body forces. The methods of Abe-Meeron [1, 52] are used to obtain a formal diagrammatic represen- tation of the Maxwell-Boltzmann quantities for this classical system. Second, the exchange contribution is evaluated perturbatively with Slater sums. In the FK repre- sentation, the exchange effects appear as open filaments within the system of closed filaments described by Maxwell-Boltzmann statistics. All of the Maxwell-Boltzmann 72 and exchange terms are collected from both steps, resulting in virial expansions of the thermodynamic functions in powers of the total density. This treatment includes the ef- fects of screening, diffraction, bound states and exchange contributions as corrections to the ideal gas EOS. Here, an expression is given for the pressure up to order ρ 5/2 . (In the theoretical development provided in this chapter, ρ refers to the number density in order to agree with the litterature. In the quantitative applications of chapters 9 and 10,ρ will refer to total density and n will be used for number density.) At this order, the EOS describes all two-body bound-state, diffusive, and exchange processes systematically, analytically, and exactly. Future work will extend the calculation to higher order, including many- body processes. 8.1 Grand-Canonical Partition Function To describe the model for the FK representation, we introduce a systemS consisting of an arbitrary number of elements. Each element α has a mass m α , a charge e α , and a spinσ α . Two charges a distance r apart interact via the standard Coulomb potential V(r)= e α e β r . (8.1) The Hamiltonian for such a system of N particles in a box with volumeΛ is H N,Λ =− X i ~ 2 2m i Δ i + 1 2 X i, j e i e j |~ r i −~ r j | (8.2) where Δ i is the Laplacian and i = [ k α ] denotes a double index. k runs from 1 to the number of charges N α of element α. α runs from 1 to the number of elements n s 73 (N = X α N α ). Dirichlet boundary conditions are imposed to ensure that the eigenfunc- tions vanish at the surface of the box. In the Grand-Canonical ensemble, the grand partition function of systemS is Ξ Λ = Tr Λ exp[−β(H N − X α μ α N α )] . (8.3) The trace is taken over all states that satisfy the boundary conditions and are sym- metrized according to the spin of each element. In the space configurational repre- sentation, the diagonal terms describe the Maxwell-Boltzmann statistics and the off- diagonal terms describe contributions from the Fermi or Bose statistics. In the thermodynamic limit (Λ → ∞ while β and μ α are kept fixed), the bulk pressure and number density of each element are βP= lim T L 1 Λ lnΞ Λ (8.4) ρ α = z α ∂ ∂z α h lim T L 1 Λ lnΞ Λ i | β (8.5) where z α = exp(βμ α ) is the activity of elementα. The grand partition function expressed in terms of the Green function is not useful for practical perturbative calculations because the kinetic and potential operators do not commute. This problem can be avoided by using the FK representation of the path integral formalism [26] described in the following sections. 8.2 Single Particle in the FK Representation To illustrate the formalism of the FK representation, we consider a simple case of a single particle of mass m, charge e, and spin σ in an external potential V(r). In or- 74 der to calculate the grand partition function, we must calculate the matrix element <~ r| exp(−βH)|~ r 0 > where|~ r > and|~ r 0 > define the ket positions, and H is the Hamilto- nian of the particle. In the FK representation, the path integral is defined as the statistical equivalent to the Feynman path integral with the time variable t replaced by an imaginary time defined as t= i~β. Then the matrix elements become <~ r| exp h −β(− ~ 2 2m Δ+ V) i |~ r 0 >= X all paths exp − S (~ r,~ r 0 ) ~ ! (8.6) where S (~ r,~ r 0 )= Z β~ 0 dt 0 n m 2 h d~ r(t 0 )/dt i 2 + V(~ r(t 0 )) o (8.7) is the classical action for the potential−V for a path~ r(t 0 ) going from~ r to~ r 0 during the imaginary time t= i~β. By making the change t = sβ~ and parameterizing the path according to ~ r (t) = ~ r+λ ~ ξ(s) withλ= (β~ 2 /m) 1/2 , equation (8.6) leads to the Feynman-Kac representation [80]: <~ r| exp(−β[− ~ 2 2m Δ+ V(r)])|~ r >= 1 (2πλ 2 ) 3/2 Z D( ~ ξ) exp[−βV ? (~ r, ~ ξ)] . (8.8) In equation (8.8),the new internal degree of freedom ~ ξ(s) is introduced in order to obtain the potential V ? (~ r, ~ ξ) from the external potential V(~ r) V ? (~ r, ~ ξ)= Z 1 0 dsV(~ r+λ ~ ξ(s)) . (8.9) With a specific choice of the new dimensionless variable ~ ξ(s), we can define a Gaussian measureD( ~ ξ) of the Brownian bridge process which defines the functional 75 Figure 8.1: A classical closed filament located at~ r with a shape parameterized by ~ ξ(s). This represents the one-body diagonal matrix contribution. (Figure provided by A. Perez) integration over all dimensionless paths ~ ξ(s) subject to the constraint ~ ξ(0)= ~ ξ(1)= ~ 0. This measure is normalized to 1 and its covariance is given by Z D( ~ ξ)ξ μ (s)ξ ν (t)= δ μν ×                s(1− t), s≤ t t(1− s), t≤ s (8.10) As seen in equation (8.8), the parameterization of the space variable ~ r (t) with the dimensionless path ~ ξ replaces the kinetic part of the Hamiltonian with a constant factor. The resulting expression is composed of only a functional of the classical potential V ? . In this context, exp[−βV ∗ ( ~ ξ )] can be interpreted as the Boltzmann factor associated with a classical closed filament located at~ r with a shape parameterized by ~ ξ(s) (see Fig.8.1). Equivalently, the quantum matrix element in the left-hand side of equation (8.8) can be expressed in terms of a functional of the classical potential. This allows us to rewrite a perturbative series of the grand-partition functionΞ which can be used for more general purposes. 76 8.3 Maxwell-Boltzmann Statistics for SystemS The FK formalism presented for the case of a single particle in section 8.2 can now be generalized to describe our mulitcomponent systemS. As described at the beginning of this chapter, this is done in two parts. In the first part, only Maxwell-Boltzmann statistics are considered. The second part, the exchange contribution, will be treated in section 8.4. 8.3.1 Maxwell-Boltzmann Grand Partition Function The grand partition function Ξ MB ∧ of systemS is given by the space-configurational expression Ξ MB Λ = ∞ X N α =0 Y α z N α α N α ! (2σ α + 1) N α × Z Λ N Y i d~ r i D ~ R N | exp(−βH N )| ~ R N E (8.11) where| ~ R N > represents the ket of the N individual states|~ r i >, | ~ R N >=⊗ i |~ r i > . (8.12) As mentioned in the previous section, equation (8.8) can be easily generalized to the N−body case. Then the FK representation of the diagonal matrix element of reads < ~ R N | exp(−βH N )| ~ R N >= Y α 1 (2πλ 2 α ) 3N α /2 Z Y i D( ~ ξ i )× exp " − β 2 X i, j e i e j v(E,E 0 ) # (8.13) where v(E,E 0 ) is the overall potential of two filaments in statesE andE 0 , andE denotes 77 (α,~ r, ~ ξ ). In this notation, the potential is v(E,E 0 )= Z 1 0 ds v c (|~ r+λ α ~ ξ(s)−~ r 0 −λ α 0 ~ ξ 0 (s)|) . (8.14) This potential is not the same as the electrostatic interaction energy between two uni- formly charged filaments, because the average of v c is taken over positions at the same “time” s. However, at large distances the modification of the shapes of the filaments during a time interval is not distinguishable and they can then be replaced by points. In this case, the potential reduces to the regular Coulomb potential, i.e. v(E,E 0 )∼ 1/|~ r−~ r 0 | when|~ r−~ r 0 |→∞. The grand-partition functionΞ MB ∧ of equation (8.11) can be rewritten using the FK matrix element expression of equation (8.13). This leads to a sum over the states of S ∗ weighted by Boltzmann factors associated with the filaments’ interactions. In this sum, the phase-space measure for each filament can be defined as dE=“dα”d~ rD( ~ ξ) with z(E)= (2σ α + 1)z α /(2πλ 2 α ) 3/2 as the fugacity. The above sum is then identified as the grand-partition functionΞ ∗ ∧ of the equivalent systemS ∗ made of classical filaments Ξ MB Λ = ∞ X N=0 1 N! Z N Y k=1 dE k z(E k ) Y k<l [1+ f (E k ,E l )] (8.15) where f (E k ,E l ) = exp[−βe α k e α l v(E k ,E l )] denotes the Mayer-bond associated with v(E k ,E l ). Notice the total equivalence between the grand partition function calcu- lated in the quantum point-particle system and the grand partition function defined within the classical filament system. In the system of classical filaments, the quantum mechanical aspect ofS MB is hidden in the complex nature of the filaments. These ex- tended objects describe quantum fluctuations of the point particles. According to this formalism, the interactions between the filaments are strictly of the two-body type. 78 Figure 8.2: Quantum point-particle systemS and the equivalent classical system of filamentsS ∗ . (Figure provided by A. Perez) This situation allows us to proceed with a perturbative expansion of the grand parti- tion function in powers of the fugacity for classical objects interacting exclusively via two-body type potentials. 8.3.2 Virial Expansion for the Coulomb Potential Now that we have an ordinary classical system of point particles with two-body in- teractions that is isomorphic to our system of closed filaments (Fig.8.2), we can study the equilibrium properties of the quantum systemS by applying classical statistical mechanics techniques to the equivalent systemS ∗ . InS ∗ , the position~ r is replaced by the generalized coordinateE, and the quantum mechanical operators are replaced by classical numbers. Since all quantities inS ∗ commute, we can apply the calculation rules of Abe [1] and Meeron [52]. Due to the long-range nature of the filament-filament potential, all the Mayer-like graphs of our system diverge. Alastuey and Perez [4, 2, 3, 5] have shown that the corresponding series can be rearranged in series of finite re-summed graphs. The main steps of the virial-like resummation method are outlined here. 79 First, the Ursell function h(E a ,E b ) ofS ∗ , which is defined by ρ(E a )ρ(E b )h(E a ,E b )= z(E a )z(E b ) lim T L δ 2 lnΞ ? Λ δz(E a )δz(E b ) , (8.16) is calculated in terms of the one-filament densities. It can be represented by a series of Mayer graphsΓ in terms of the filament density ρ(E)= z(E) lim T L δ lnΞ ? Λ δz(E) (8.17) The graphsΓ are defined via familiar topological prescriptions, where the usual points are now replaced by filaments. Each graph is built with the two root filamentsE a andE b and n field filamentsE 1 ,...,E n which are integrated over. Each field filament E i has a statistical weight ρ(E i ). Two filaments are linked by at most one f -bond. According to this topological prescription, each graph Γ is connected and does not contain articulation filaments. It appears that the contribution of each graph Γ is divergent because of the non- integrable 1/r-decay behaviour of the Mayer bond f . Following the procedure in- troduced by Abe and Meeron, the resummation procedure starts with the following decomposition f (E,E 0 ) = f C (|~ r−~ r 0 |)+ 1 2 f 2 C (|~ r−~ r 0 |) + Z 1 0 ds[λ α ~ ξ(s). ~ ∇+λ α 0 ~ ξ 0 (s). ~ ∇ 0 ] f C (|~ r−~ r 0 |)+ f T (E,E 0 ) (8.18) where f C is the shape-independent Coulomb bond f C (|~ r−~ r 0 |) =−βe α e α 0v C (|~ r−~ r 0 |). From this useful decomposition, the truncated bond f T appears. By construction f T decays as 1/|~ r−~ r 0 | 3 when|~ r−~ r 0 |→∞. Inserting Equation (8.18) in each graphΓ, 80 we obtain a new representation of the Ursell function h(E a ,E b ) in terms of graphs ˜ Γ built with bonds ˜ f which may be either f C , f 2 C / 2,λ ~ ξ. ~ 5 f C or f T . Since f T is almost integrable, the divergences in the ˜ Γ-diagrams are induced only by the other non-integrable bonds ˜ f . They are of the same type as those encountered in the purely classical case and can be eliminated by using the mathematical recipe in- troduced by Mayer [51] and Salpeter [72], i.e., the resummation of all the convolution chains built with the Coulomb bond f C . This procedure transforms the whole set of ˜ Γ-diagrams into a set of new graphsΠ built with re-summed bonds named generically F. It turns out that for anyΠ graph, only four kinds of re-summed bonds F have to be considered. The single chains with ending bonds which are either f C or λ i ~ ξ i . ~ 5 i f C (λ j ~ ξ j . ~ 5 j f C ), lead to three re-summed bonds F D , (λ ~ ξ. ~ 5 F D ) and F dip and the other structures, involving several chains and/or the bonds f T and f 2 C /2, give raise to a fourth bond F R . The resummation procedure automatically excludes specific association of any combination of F re-summed bonds. At this level of the diagrammatic expansion, the re-summed bonds F can be calculated explicitly in terms of the MB particle densities. Thus, the summation of all the convolution chains can be performed in terms of the familiar Debye potential φ D (r) = exp (−κr)/r with κ = (4πβ X α e 2 α ρ MB α ) 1/2 . In particular, the fist re-summed bond F D reads: F D (E i ,E j )=−βe α i e α j φ D (|~ r i −~ r j |) (8.19) In contrast to F D andλ ~ ξ. ~ 5F D which decay exponentially fast (likeφ D essentially), the bonds F dip and F R are found to decay algebraically as 1/r 3 when r→∞. The re-summed diagrammatic expansion of the two-point correlations ofS is im- mediately obtained by inserting the previousΠ-representation of the Ursell functional 81 h(E a ,E b ) in the identity ρ MB T (α a ~ r a ;α b ~ r b )= Z D( ~ ξ a )D( ~ ξ b )ρ(E a )ρ(E b )h(E a ,E b ) (8.20) In Equation (8.20), each graph Π is multiplied by ρ(E a ) ρ(E b ) and integrated over the shapes ~ ξ a and ~ ξ b of the root filamentsE a andE b . If we let ~ → 0, the bonds λ~ξ. ~ 5F D and F dip disappear, while the bond F R reduces to (exp(F D )− 1− F D ) and the statistical weights ρ(E) are replaced by the particle densities. We then recover the nodal expansion of the classical correlations first derived by Meeron [52]. 8.4 Quantum Statistics Now that we have a formal diagrammatic representation of the Maxwell-Boltzmann quantities of the classical systemS ∗ , we can proceed with step two - evaluating the exchange contribution. In order to do this, we must properly take into account the contributions to the thermodynamic functions from the Fermi and Bose statistics. 8.4.1 Quantum Grand-Canonical Partition Function Consider Equation (8.3) of the grand-cananical partition function over the properly symmetrized states| ~ R N σ N > s in configuration and spin spaces. These states are defined via the usual Slater sums | ~ R N σ z N > S = 1 Q α (N α !) 1/2 X P α Y α α (P)⊗ i |~ r P α (i) σ z P α (i) > S (8.21) P α is a permutation of (1...N α ),P α (i) = (P α (k),α), andE α (P α ) is either 1 if the par- ticles of species α are bosons (σ α integer) or the signature (±1) ofP α if the particles 82 are fermions (σ α half integer).⊗ means a tensional product over the one-body states |~ rσ z > describing a particle localized at~ r with the projection of its spin along a given Z-axis equal toσ z (σ z may take (2σ α + 1) values). Using equation (8.21), we obtain Ξ Λ = ∞ X N α =0 Q α z N α α [ Q α (N α !)] 2 × X P α ,P 0 α Y α α (P α ) α (P 0 α ) X {σ z i } Y i < σ z P 0 α (i)|σ z P α (i) > × Z Λ N Y i d~ r i <~ r P 0 α (i)|⊗ i | exp(−βH N,Λ )|⊗ i |~ r P α (i) > (8.22) Since the Hamiltonian H N does not depend on the spins, the spin-part of the matrix elements contributes the trivial degeneracy factor X σ z i Y i < σ z P 1 α (i) |σ z P α (i) > which only depends on the pairs of permutations (P α ,P 0 α ). The Slater-sum representation of Equation (8.22) of Ξ ∧ allows a natural identifi- cation of the exchange effects. Indeed, the square terms (P α =P 0 α for any α), where the diagonal matrix elements of exp(−βH N ) in configuration space appear, obviously correspond to MB statistics. An off-diagonal term (P α ,P 0 α for at least one species) involves the exchange of n particles (n≥ 2). The corresponding matrix elements of exp(−βH N ) are off-diagonal with respect to the positions of the exchanged particles. In the language of the filaments, the structure of the FK representation of these off- diagonal matrix elements can be interpreted as opened filamentsF α kl associated with the exchange of a particleα from position~ r k to position~ r l . The shape ofF α kl is param- eterized by ~ ω α kl (s)= (1− s)~ r k + s~ r l +λ α ~ ξ(s) which describes a path of the exchanged particle in the genuine Feynman path integral (~ ω kl (0)=~ r k and ~ ω kl (1)=~ r l ). 83 The above closed filamentsE are again associated with the non-exchanged parti- cles. For instance, if one considers the off-diagonal matrix element <~ r 2 ~ r 1 ~ r 3 ...~ r N | exp(−βH N )|~ r 1 ~ r 2 ~ r 3 ...~ r N > (8.23) which corresponds to the exchange of two particles, two opened filamentsF α 12 andF α 21 and (N− 2) closed filamentsE z ,...,E N appear. This inhomogeneous situation can be dealt with using standard perturbative techniques where the reference system is the homogeneous bathS ∗ described in the previous section. 8.4.2 Exchange Contributions The exchange effects can now be taken into account by using the Slater expansion of Ξ ∧ from equation (8.3) in the grand-canonical expression of βP from equation (8.4). This gives βP= βP MB + ∞ X n=2 E n , (8.24) where E n is the contribution of n exchanged particles. In the impurities point of view, E n can be expressed in terms of the densityρ ∗ (E a |F kl ) ofS ∗ in the presence of n opened filamentsF kl . To compute the exchange contribution up to order ρ 5/2 (two-body inter- actions), we only need to calculate the first term E 2 E 2 = − 1 2 X α (−1) 2σ α +1 (2σ α + 1) z 2 α (2πλ 2 α ) 3 × Z Λ d~ r 12 Z D( ~ ξ 1 )D( ~ ξ 2 )× exp h −r 2 12 λ 2 α −βv(F α 12 ,F α 21 ) i × exp h −β Z 1 0 dg Z dE a ρ(E a |g;F α 12 ,F α 21 ) v(E a ,F α 12 )+ v(E a ,F α 21 ) i (8.25) 84 The first exponential represents the exchange contribution in a vacuum. The second exponential describes the many-body effects on the two-particle exchange. All other E n ’s have a similar structure. The diagrammatical method exposed in section 8.3 can be extended to the inho- mogeneous systemS ∗ in the presence of theF ’s. This provides a representation of ρ(E a |F kl ) in terms of graphs made of the n opened filamentsF kl and closed filaments E a . The statistical weight of a closed filament is the density of the homogeneous sys- temρ(E). Two closed filaments are linked by at most one re-summed bond F. A closed filamentE and an opened filamentF are linked by at most one external Mayer bond (exp [−βv(E,F )]− 1). The density expansion of E n follows from the use of the previous diagrammatic representation of ρ(E a |F kl ). The integrals over the positions of the opened and closed filaments can be evaluated via a scaling analysis with respect to κ D . The resulting expression for E n takes the general form z n multiplied by a double integer series in (ρ MB ) 1/2 and in lnρ MB . 8.5 Density Expansions We can now collect all the Maxwell-Boltzmann and exchange contributions for the virial-like expansion of the thermodynamic functions. We go back to the expression for the total pressure βP= βP MB + ∞ X n=2 E n . (8.26) The Maxwell-Boltzmann pressure can be found from βP MB = X α ρ MB ∂ ∂ρ MB lim T L h βF MB Λ Λ i | β − lim T L h βF MB Λ Λ i (8.27) 85 where the free energy per unit volume can be expressed in terms of the two-point correlations as βF MB Λ = βF MB id Λ + β 2 X α,β Z 1 0 dg Z d~ rρ MB T,g (α ~ 0,β~ r)e α e β v C (r) . (8.28) ρ MB T,g is calculated for a fictitious systemS g where all the interactions are multiplied by the dimensionless coupling parameter g. In addition, the temperature and the Maxwell- Boltzmann particle densities ofS g are identical to those ofS. A diagrammatic expan- sion of the free energy is obtained by inserting theΠ-representation ofρ MB T,g . We can now combine the Maxwell-Boltzmann pressure from equation 8.27 with the exchange contribution E n to create an expression for the total pressure. 8.5.1 Equation of State up to Orderρ 5/2 We now arrive at the goal of this application of the FK formalism - an equation of state expressing pressure as a function of temperature and density. We limit the current treatment to two-body interactions, truncating the virial expansion of the pressure at 86 orderρ 5/2 [4] βP = X α ρ α − κ 3 D 24π + π 6 (ln 2− 1) X α,β β 3 e 3 α e 3 β ρ α ρ β − π √ 2 X α,β ρ α ρ β λ 3 αβ Q(x αβ )− π 3 β 3 X α,β ρ α ρ β e 3 α e 3 β ln (κ D λ αβ ) + π √ 2 X α (−1) 2σ α +1 (2σ α + 1) λ 3 αα ρ 2 α E(x αα )− 3π 2 √ 2 β X α,β e α e β κ D ρ α ρ β λ αβ 3 Q(x αβ ) − π 2 β 4 X α,β ρ α ρ β e 4 α e 4 β κ D ln (κ D λ αβ )+ 3π 2 √ 2 β X α (−1) 2σ α +1 (2σ α + 1) λ 3 αα ρ 2 α e 2 α κ D E(x αα ) + 1 16 X α β 2 ~ 2 e 2 α m α κ 3 D ρ α +π( 1 3 − 3 4 ln 2+ 1 2 ln 3)× X α,β β 4 e 4 α e 4 β κ D ρ α ρ β +C 1 X α,β,γ β 5 e 3 α e 4 β e 3 γ κ −1 D ρ α ρ β ρ γ + C 2 X α,β,γ,δ β 6 e 3 α e 3 β e 3 γ e 3 δ κ −3 D ρ α ρ β ρ γ ρ δ (8.29) with κ D = (4πβ X α e 2 α ρ α ) 1/2 , l αβ = βe α e β , m αβ = m α m β /(m α + m β ), λ αβ = (β~ 2 /m αβ ) 1/2 , x αβ = − √ 2 l αβ λ αβ , C 1 = 15.205 ± .001, C 2 = −14.733 ± .001 [62], and the Euler- Mascheroni constant C = 0.577216. Q(x αβ ) is the quantum second-virial coefficient first introduced by Ebeling et al [23, 43]. Q(x αβ )= 1 ( √ 2πλ 3 αβ ) lim R→∞ (Z r<R d~ r h (2πλ 2 αβ ) 3/2 <~ r|e −βh αβ |~ r >−1+ βe α e β r − β 2 e 2 α e 2 β 2r 2 i + 2π 3 β 3 e 3 α e 3 β [ln( 3 √ 2R λ αβ )+ C] ) (8.30) while E(x αα ) is defined as the two-body exchange integral, E(x αα )= (2 √ π) lim R→∞ Z r<R d~ r <−~ r|e −βh αα |~ r > (8.31) 87 In Equations (8.30) and (8.31), h αβ is the one-body Hamiltonian of the relative particle with mass m αβ subject to the Coulomb potential e α e β /r. The functions Q and E only depend on the temperature via the single dimensionless parameter (−l/λ). 88 Chapter 9 Our Analysis of Applicability in the Solar Interior The formalism for the FK EOS is now in place. Before applying it to the Sun, its domain of validity must be examined. To facilitate this analysis, the following four characteristic lengths are defined: d α = 3 4πn α ! 1/3 mean interparticle distance (9.1) λ αβ = β~ 2 m αβ ! 1/2 de Broglie wavelength (9.2) l αβ = βe α e β Landau length (9.3) a 0 = ~ 2 me 2 Bohr radius (9.4) The indicesα andβ indicate electrons, protons, and helium nuclei. For the quantitative solar applications that follow, number density is denoted by n instead of ρ, which was used in chapter 8 for the theoretical development. 89 Figure 9.1: Coupling parameter throughout the Sun. 9.1 Coupling The first condition for application of the FK EOS is that the plasma must be weakly coupled Γ α << 1 (9.5) where Γ α = l αα d α (9.6) is the coupling parameter. Though it is well-known that the plasma inside the Sun is weakly coupled, this result is verified using a standard solar model (Model S from Christensen-Daalsgaard, as described by Christensen-Daalsgaard et al [16]). As seen in Fig.9.1, the coupling parameter throughout the Sun is safely below the critical value of one. 90 (a) electrons (b) protons (c) helium Figure 9.2: Degeneracy throughout the Sun for (a) electrons, (b) protons, and (c) he- lium nuclei. 9.2 Degeneracy The second condition for application of the FK EOS is that the plasma is non-degenerate. This is the case if the degeneracy parameter n α Λ 3 α is sufficiently low n α Λ 3 α << 1 . (9.7) Λ α is the thermal de Broglie wavelength defined by Λ α = h √ 2πm α k B T . (9.8) As seen in Fig.9.2, the degeneracy of protons, electrons, and helium nuclei is low throughout the solar interior. 9.3 Ionization and Recombination The third and most restrictive condition for application of the FK EOS is full or nearly full ionization. Recombination prevents use of the low-density virial expansion in the 91 Figure 9.3: Hydrogen ionization throughout the Sun. outer layers of the Sun. In the interior, the region in which the plasma is sufficiently ionized for the virial expansion is determined by computing the standard Saha equation throughout the Sun. N z+1 N e − N z = g e g z+1 g z ! V Λ 3 e exp −χ z k B T ! (9.9) N α is the number of particles of species with ion charge α. g α is the statistical weight. χ α is the ionization potential for the reaction Z→ Z+1. Figures 9.3, 9.4, and 9.5 show the ionization throughout the Sun for hydrogen, helium, and the first two ionization stages of oxygen (representing a typical heavy element). In these figures the Saha equation is solved for the ratio of the particle densities of two different ionization levels x z = N z+1 N 0 z . (9.10) In the cool outer layers, the virial expansions break down due to recombination. How- ever, the plasma is at least 90% ionized for log(T) > 5. Fig.9.6 of the temperature at each depth shows that this temperature gange includes most of the Sun. Thus the FK EOS can be used throughout a large part of the solar interior. 92 Figure 9.4: Helium ionization throughout the Sun. Figure 9.5: Oxygen ionization throughout the Sun. Figure 9.6: log(T) for each depth r/R in the Sun. 93 9.4 Domain of Applicability As seen in the previous sections, the conditions of weak coupling and low degeneracy are met throughout the Sun. The ionization condition is also met in the hot inner layers, though the exact boundary between ionized plasma and a region of recombination will vary slightly depending on which heavy elements are included in the analysis. Despite the limitation of nearly full ionozation, the FK formalism can be applied to develop equations of state that are valid in a large part of the solar interior. For completelness, the FK EOS applied to the inner regions can be blended with a more conventional EOS such as OPAL or MHD applied in the outer layers. The localizing nature of helioseismology enables such a patchwork method, which would allow each region of the Sun to be described by the formalism best suited to the conditions of the region. This will lead to better solar models and better constraints on the solar EOS. 94 Chapter 10 Results of Application to Solar Interior Now that we have determined the appropriate region for this formalism, the FK EOS can be applied to the solar interior. We begin with the EOS expression for the pressure up to order 5/2 in density from equation (8.29) βP = X α n α − κ 3 D 24π + π 6 (ln 2− 1) X α,β β 3 e 3 α e 3 β n α n β − π √ 2 X α,β n α n β λ 3 αβ Q(x αβ )− π 3 β 3 X α,β n α n β e 3 α e 3 β ln (κ D λ αβ ) + π √ 2 X α (−1) 2σ α +1 (2σ α + 1) λ 3 αα n 2 α E(x αα )− 3π 2 √ 2 β X α,β e α e β κ D n α n β λ αβ 3 Q(x αβ ) − π 2 β 4 X α,β n α n β e 4 α e 4 β κ D ln (κ D λ αβ )+ 3π 2 √ 2 β X α (−1) 2σ α +1 (2σ α + 1) λ 3 αα n 2 α e 2 α κ D E(x αα ) + 1 16 X α β 2 ~ 2 e 2 α m α κ 3 D n α +π( 1 3 − 3 4 ln 2+ 1 2 ln 3) X α,β β 4 e 4 α e 4 β κ D n α n β +C 1 X α,β,γ β 5 e 3 α e 4 β e 3 γ κ −1 D n α n β n γ + C 2 X α,β,γ,δ β 6 e 3 α e 3 β e 3 γ e 3 δ κ −3 D n α n β n γ n δ . (10.1) 95 electron proton helium nucleus oxygen nucleus m (kg) 9.1095D−31 1.6726D−27 6.6448D−27 2.6561D−26 q (C) −1.6022D−19 1.6022D−19 3.2043D−19 1.2818D−18 spin .5 .5 0 0 Table 10.1: Mass, charge, and spin for the elements in our EOS model. For these quantitative applications, n α is used to represent number density in order to distinguish it from the mass density which will be denoted asρ α . 10.1 Developing the Model In order to evaluate this EOS throughout the Sun, we must first define the elements α for our model. We will take the constituents of our plasma to be electrons and nuclei of hydrogen (about 70% of the Sun’s mass), helium (about 28% of the Sun’s mass), and oxygen (representing the remaining heavy elements in the Sun). We define the charge e α , mass m α , and spinσ α of each element in table 10.1. Then we can calculate the number densities n α , de Broglie wavelengths λ αβ , and Debye radius κ D at each depth in the Sun. To do this, we use the standard solar model provided by Christensen-Daalsgaard [16] for the temperature T, mass density ρ, hy- drogen abundance X, and heavy element abundance Z for each value of the radius r/R (R =solar radius) in 2,481 increments between 0 and 1.0007123. The helium abun- dance Y at each depth is computed from the normalization condition X+ Y + Z = 1. The number density of each nucleus for a given depth is n α = ρ A α m α (10.2) where A α is the mass abundance X, Y, or Z. Since the plasma is electrically neutral, the 96 number density of electrons is n e = n H + 2n He + 8n O . (10.3) The de Broglie wavelengths at each depth are then computed via λ αβ = β~ 2 m αβ ! 1/2 (10.4) where m αβ is the reduced mass of a pair of particles. The Debye radiusκ D at each depth is computed via κ= (4πβ X α e 2 α n α ) 1/2 . (10.5) 10.2 Improving the Precision of the Model The precision of the FK EOS relies on the precision of the parameters used in the model. Perez [62] computed the constants for the virial expansion using Simpson’s method out to three decimal places. This set the precision limit for all of the calcu- lations done in the model. Solar data, however, is precise to four decimal places. In order to improve the precision of the model to meet the precision of the data, I used Maple to calculate the multidimensional integrals that contribute to these constants. For the constants C 1 and C 2 , we need C 1 = 3 2 I 1 (10.6) C 2 = 3 2 I 2 (10.7) 97 where I 1 = I d + I f 1 + 2I 1,1 + I 2,1 (10.8) I 2 = I e + I g1 + 4I 1,2 + 2I 2,2 + I 0,1,1 + 2I 1,1,1 + I 2,1,1 (10.9) require the calculation of the integrals I d = 8 9 π Z ∞ 0 dk arctan 2 k 2 ! k 2 (k 2 + 1) 2 ! (10.10) I e = −32 9 π 2 Z ∞ 0 dk arctan 2 k 2 ! k 2 (k 2 + 1) 3 ! (10.11) I f 1 = 8 9 π Z ∞ 0 dk 4 (k 2 + 1) arctan k 2 ! arctan(k) (10.12) I g1 = −32 9 π 2 Z ∞ 0 dk 2 (k 2 + 1) 2 arctan k 2 ! arctan(k) (10.13) I m,n = (−1) m+n+1 π n−1 16 9 Z ∞ 0 dk Z ∞ 0 dk 0 1 k 1 (k 2 + 1) m 1 (k 02 + 1) n (10.14) × arctan k 0 2 ! ln 1+ (k+ k 0 ) 2 1+ (k− k 0 ) 2 ! . (10.15) I m,n,p = (−1) m+n+p+1 π 3 64 9 ∞ X 0 (2b+ 1) Z ∞ 0 dk 1 Z ∞ 0 dk 2 Z ∞ 0 dk 3 (10.16) × 1 k 2 1 (k 2 1 + 1) m (k 2 2 + 1) n (k 2 3 + 1) p Q b (x 12 )Q b (x 13 )Q b (x 23 ) . (10.17) Q b (x i j ) are Legendre funtions of order b and x i j is a dimensionless parameter related to the integration variables by x i j = 1+ k 2 i + k 2 j 2k i k j . (10.18) 98 Perez value Mussack value I d 1.967 1.96736308 I e -2.195 -2.19514005 I f 1 11.218 11.21819449 I g1 -11.156 -11.15584777 I 1,1 −2.206 -2.20950196 I 1,2 1.999 1.99889526 I 2,1 1.366 1.36822544 I 2,2 −1.348 -1.34870382 C 1 15.205 15.20216865 Table 10.2: New values for the virial expansion constants. Values for I d , I e , I f 1 , I g1 , I m,n , and C 1 are shown in table 10.2. Evaluation of the three-dimensional integrals I m,n,p is still necessary to compute C 2 . The value of C 2 = −14.733 computed by Perez only includes the first five terms of the sum in equation (10.16). Higher order terms will need to be included in order to improve the precision of the C 2 calculation. 10.3 Applying the Formalism Although the formalism has been developed to treat a multi-component plasma, we begin by testing it on a simpler model of only electrons and hydrogen nuclei (X = 1) with solar temperatures and densities. In the sections that follow, we present results for thermodynamic quantities of this hydrogen plasma as determined by our application of the FK EOS. 99 Figure 10.1: Results for the pressureβP/n at each depth in the Sun. 10.4 Results for Pressure Now that we have the relevant parameters for the EOS, we can compute the pressure at each depth from equation (10.1). The results forβP/n are shown in Fig.10.1. 10.5 Free Energy, Internal Energy, and Entropy The first thermodynamic function we will address is the free energy per unit volume f = F/V which can be found by direct integration of the identity βP= X α n α ∂(β f ) ∂n α −β f . (10.19) 100 This results in the following expression for f up to orderρ 5/2 : β f = βF V = X α n α ln h n α (2πλ 2 α ) 3/2 i − 1 − κ 3 D 12π + π 6 ln 2 X α,β β 3 e 3 α e 3 β n α n β − π √ 2 X α,β n α n β λ 3 αβ Q(x αβ )− π 3 β 3 X α,β n α n β e 3 α e 3 β ln (κ D λ αβ ) + π √ 2 X α (−1) 2σ α +1 (2σ α + 1) λ 3 αα n 2 α E(x αα )− π √ 2 β X α,β e α e β κ D n α n β λ αβ 3 Q(x αβ ) − π 3 β 4 X α,β n α n β e 4 α e 4 β κ D ln (κ D λ αβ )+ π √ 2 β X α (−1) 2σ α +1 (2σ α + 1) λ 3 αα n 2 α e 2 α κ D E(x αα ) + 1 24 X α β 2 ~ 2 e 2 α m α κ 3 D n α +π( 1 3 − 1 2 ln 2+ 1 3 ln 3) X α,β β 4 e 4 α e 4 β κ D n α n β + 2 3 C 1 X α,β,γ β 5 e 3 α e 4 β e 3 γ κ −1 D n α n β n γ + 2 3 C 2 X α,β,γ,δ β 6 e 3 α e 3 β e 3 γ e 3 δ κ −3 D n α n β n γ n δ . (10.20) We can also find an expression for the internal energy per unit volume from the identity u= U V = ∂ ∂β βF V ρ . (10.21) 101 Differentiating equation (10.20), we find βu= βU V = 3 2 X α n α − κ 3 D 8π + π 2 ln 2− 2 3 ! X α,β β 3 e 3 α e 3 β n α n β − π √ 2 X α,β n α n β λ 3 αβ " 3 2 Q(x αβ )+β ∂ ∂β Q(x αβ ) # −πβ 3 X α,β n α n β e 3 α e 3 β ln (κ D λ αβ ) + π √ 2 X α (−1) 2σ α +1 (2σ α + 1) λ 3 αα n 2 α " 3 2 E(x αα )+β ∂ ∂β E(x αα ) # − π √ 2 β X α,β e α e β κ D n α n β λ αβ 3 " 3Q(x αβ )+β ∂ ∂β Q(x αβ ) # − 3π 2 β 4 X α,β n α n β e 4 α e 4 β κ D ln (κ D λ αβ ) + π √ 2 β X α (−1) 2σ α +1 (2σ α + 1) λ 3 αα n 2 α e 2 α κ D " 3E(x αα )+β ∂ ∂β E(x αα ) # + 7 48 X α β 2 ~ 2 e 2 α m α κ 3 D n α +π( 7 6 − 9 4 ln 2+ 3 2 ln 3) X α,β β 4 e 4 α e 4 β κ D n α n β +3C 1 X α,β,γ β 5 e 3 α e 4 β e 3 γ κ −1 D n α n β n γ + 3C 2 X α,β,γ,δ β 6 e 3 α e 3 β e 3 γ e 3 δ κ −3 D n α n β n γ n δ . (10.22) We can also find the entropy per unit volume from the identity s= S V = k B β 2 ∂ f ∂β ! ρ = k B " β ∂(β f ) ∂β −β f # (10.23) 102 which leads to the expression s k B = X α n α 5 2 ln h n α (2πλ 2 α ) 3/2 i ! − κ 3 D 24π + π 3 (ln 2− 1) X α,β β 3 e 3 α e 3 β n α n β − π √ 2 X α,β n α n β λ 3 αβ " 1 2 Q(x αβ )+β ∂ ∂β Q(x αβ ) # − 2π 3 β 3 X α,β n α n β e 3 α e 3 β ln (κ D λ αβ ) + π √ 2 X α (−1) 2σ α +1 (2σ α + 1) λ 3 αα n 2 α " 1 2 E(x αα )+β ∂ ∂β E(x αα ) # − π √ 2 β X α,β e α e β κ D n α n β λ αβ 3 " 2Q(x αβ )+β ∂ ∂β Q(x αβ ) # − 7π 3 β 4 X α,β n α n β e 4 α e 4 β κ D ln (κ D λ αβ ) + π √ 2 β X α (−1) 2σ α +1 (2σ α + 1) λ 3 αα n 2 α e 2 α κ D " 2E(x αα )+β ∂ ∂β E(x αα ) # + 5 48 X α β 2 ~ 2 e 2 α m α κ 3 D n α +π( 5 6 − 7 4 ln 2+ 7 6 ln 3) X α,β β 4 e 4 α e 4 β κ D n α n β + 7 3 C 1 X α,β,γ β 5 e 3 α e 4 β e 3 γ κ −1 D n α n β n γ + 7 3 C 2 X α,β,γ,δ β 6 e 3 α e 3 β e 3 γ e 3 δ κ −3 D n α n β n γ n δ . (10.24) 10.6 Results for Specific Heats Specific heats describe how much heat must be added in order to raise the temperature of the system by one unit under given conditions. Specific heats generally take the form c x = dQ dT ! x (10.25) where x is kept fixed while T changes. Q is related to the internal energy, pressure, and specific volume (V ρ = 1/ρ) by the first law of thermodynamics dQ= dU+ PdV ρ (10.26) 103 so that the specific heat for a constant specific volume is c v ρ = dQ dT ! ρ = ∂U ∂T ! ρ . (10.27) Differentiating equation (10.22), we find c v k B V = 3 2 X α n α − κ 3 D 16π +π 5 3 − ln 2 ! X α,β β 3 e 3 α e 3 β n α n β − π √ 2 X α,β n α n β λ 3 αβ " − 3 4 Q(x αβ )− 3β ∂ ∂β Q(x αβ )−β 2 ∂ 2 ∂β 2 Q(x αβ ) # +2πβ 3 X α,β n α n β e 3 α e 3 β ln (κ D λ αβ ) + π √ 2 X α (−1) 2σ α +1 (2σ α + 1) λ 3 αα n 2 α " − 3 4 E(x αα )− 3β ∂ ∂β E(x αα )−β 2 ∂ 2 ∂β 2 E(x αα ) # − π √ 2 β X α,β e α e β κ D n α n β λ αβ 3 " −6Q(x αβ )− 6β ∂ ∂β Q(x αβ )−β 2 ∂ 2 ∂β 2 Q(x αβ # − 21π 4 β 4 X α,β n α n β e 4 α e 4 β κ D ln (κ D λ αβ ) + π √ 2 β X α (−1) 2σ α +1 (2σ α + 1) λ 3 αα n 2 α e 2 α κ D " −6E(x αα )− 6β ∂ ∂β E(x αα )−β 2 ∂ 2 ∂β 2 E(x αα ) # − 35 96 X α β 2 ~ 2 e 2 α m α κ 3 D n α +π( 31 12 − 63 8 ln 2+ 21 4 ln 3) X α,β β 4 e 4 α e 4 β κ D n α n β − 21 2 C 1 X α,β,γ β 5 e 3 α e 4 β e 3 γ κ −1 D n α n β n γ − 21 2 C 2 X α,β,γ,δ β 6 e 3 α e 3 β e 3 γ e 3 δ κ −3 D n α n β n γ n δ . (10.28) Values for c v throughout the Sun are shown in Fig.10.2. We can then find the specific heat for a constant pressure from values we have already computed by using the relation c P = c v + P ρT χ 2 T χ ρ . (10.29) 104 Figure 10.2: Results for c v at each depth in the Sun. whereχ T andχ ρ are defined as χ T = ∂ ln P ∂ ln T ! ρorV = T P ∂P ∂T ! ρorV ρ (10.30) χ ρ = ∂ ln P ∂ lnρ ! T = ρ P ∂P ∂ρ ! T (10.31) In a monoatomic ideal gas, χ T = χ ρ = 1. Values for χ T and χ ρ throughout the Sun are shown in figures 10.3 and 10.4. Values for c P are shown in Fig.10.5. 10.7 Results for the Adiabatic Exponents The dimensionless adiabatic exponents measure the thermodynamic response of a sys- tem to adiabatic changes. They are used extensively in stellar modeling and are defined 105 Figure 10.3: Results forχ T at each depth in the Sun. Figure 10.4: Results forχ ρ at each depth in the Sun. 106 Figure 10.5: Results for c p at each depth in the Sun. by the following expressions: Γ 1 = ∂ ln P ∂ lnρ ! ad (10.32) Γ 2 Γ 2 − 1 = ∂ ln P ∂ ln T ! ad = 1 ∇ ad (10.33) Γ 3 − 1= ∂ ln T ∂ lnρ ! ad (10.34) where∇ ad is the adiabatic gradient defined by this equation that can also be expressed as ∇ ad = P Tc p χ T χ ρ . (10.35) For a monoatomic ideal gas,Γ 1 =Γ 2 =Γ 3 = 5/3. First we examineΓ 1 throughout the Sun. Instead of using equation (10.32),Γ 1 can 107 Figure 10.6: Results forΓ 1 at each depth in the Sun. be computed from quantities we have already found by using the expression Γ 1 = χ T P ρT χ T c vρ ! +χ ρ . (10.36) The results forΓ 1 throughout the Sun are shown in Fig.10.6. Second we evaluate∇ ad throughout the Sun using equation (10.35). The results are shown in Fig.10.7. Third we examineΓ 2 throughout the Sun. From equation (10.33) we have Γ 2 = 1 1−∇ ad , (10.37) The results ofΓ 2 are shown in Fig.10.8. Finally we evaluate the last adiabatic exponentΓ 3 throughout the Sun. Since the 108 Figure 10.7: Results for∇ ad at each depth in the Sun. Figure 10.8: Results forΓ 2 at each depth in the Sun. 109 Figure 10.9: Results forΓ 3 at each depth in the Sun. adiabatic exponents are not all independent,Γ 3 can be expressed in terms ofΓ 1 andΓ 2 . Γ 3 − 1 Γ 1 = Γ 2 − 1 Γ 2 =∇ ad (10.38) so that Γ 3 = 1+ (Γ 2 − 1)Γ 1 Γ 2 . (10.39) The results ofΓ 3 are shown in Fig.10.9. 110 Chapter 11 Extensions of EOS The solar EOS developed within the FK formalism is still fairly new. Although it has the potential to more accurately describe the central region of the Sun, the FK EOS does not yet match solar data as well as the established MHD and OPAL equations of state. There is much work to be done to improve the model. A few planned extensions are discussed in the following sections. 11.1 Relativistic Correction At the extreme temperatures of the solar core, relativistic effects become relevant. Although thermal protons have an average speed of only v p = .002c, the average speed of electrons in the core is v e = .087c. The current model does not treat the electrons relativistically. Adding relativistic effects to the EOS should yield a correction of a few percent. 111 11.2 Saha for Ionization The main limitation of the FK EOS is the restriction to fully- or nearly fully-ionized plasma. This limits its application to the center region of the Sun. The range can be extended by building the Saha equation N z+1 N e − N z = g e g z+1 g z ! V Λ 3 e exp −χ z k B T ! (11.1) into the formalism to compute the degree of ionization at each depth. Here N α is the number of particles of the element with ion chargeα. g α is the statistical weight. χ α is the ionization potential for the reaction Z→ Z+ 1. To see how the inclusion of chemical reactions would alter the pressure and internal energy of a system, consider the simple example of an ideal hydrogen gas [31]. In this case, the Saha equation can be written as N + N e − N 0 = V Λ 3 e exp −χ H k B T ! (11.2) where N + is the number density of ionized hydrogen and N 0 is the number density of neutral hydrogen such that the total number density n = N + + N 0 , and N e − = N + to preserve charge neutrality. We can define the total ion plus neutral atom number density per unit mass as N = n/ρ, which is independent of the density. Then defining y= N + /n= N e −/n, we can write the pressure as P= (N e −+ N + + N 0 )k B T = (1+ y)Nρk B T (11.3) 112 and the internal energy as U = (1+ y) n ρ 3k B T 2 + y n ρ χ H (11.4) or U = (1+ y)N 3k B T 2 + yNχ H . (11.5) With these expressions for pressure and internal energy, we can now compute the ther- modynamic derivatives as done in chapter 9. A more complete development will need to be done in order to adapt this extension to the multi-component solar plasma. 11.3 More Chemical Elements The current model includes only hydrogen (about 70% of solar mass), helium (about 28% of solar mass), and oxygen (as a representative of the remaining heavy elements in the Sun). Once we reach the desired precision for this model, we can extend the system to include more chemical elements. Candidates for inclusion are elements that are sufficiently abundant in the Sun, such as carbon, oxygen, nitrogen, neon, and iron, as well as elements that are relevant for interesting nuclear reactions, such as beryllium. The Saha equation would play an important role in this extension to determine the ionization levels of each new element. Adapting the model to include more elements can also be a useful tool for ex- amining the implications of the new solar chemical abundances [46, 8]. There is a discrepency between the recently suggested values and solar theories develeped with the previously accepted values [30]. This extension of the FK EOS could address these issues and help resolve the controversy. 113 11.4 Extension to Orderρ 3 The most complex extension of the FK EOS is to include higher order terms in the expansions. The current development includes all two-body interactions and is thus exact up to order ρ 5/2 . There are many more three-body effects, so extending the EOS to order ρ 3 would yield an explosion in the number of terms. At this point, it would be wise to sacrifice exactness in favor of practicality. The contribution of each of the third-order terms would need to be analyzed. Only dominant contributors would be kept. This extension would increase the precision of our calculations, hopefully improving agreement with solar data. 11.5 Stellar Modeling The FK EOS formalism is designed to treat fully-ionized quantum plasmas. Thus it is ideally suited for stellar matter. Though the Sun is our closest and therefore most familiar star, similar physics is happening in stars throughout the universe. Just as he- lioseismology provides reliable data on the Sun, the relatively new field of asteroseis- mology has brought a wealth of data on other stars. Once the EOS is fully developed for the Sun, a natural extension is to apply this treatment to other stars, leading to more accurate descriptions of stellar interiors. 114 Conclusion The research presented in this thesis has furthered our understanding of the solar interior. The molecular dynamics work has brought new insight to the issue of screen- ing in nuclear reactions and the equation of state work has provided a new tool for precision solar modeling. Conclusions from the molecular dynamics work come from analyzing the pair cor- relation functions of the simulation systems. We found that our pair correlation func- tions differ greatly from those derived from the Debye-H¨ uckel theory. Much of this difference can be attributed to quantum effects. However, even the quantum correc- tions do not fully account for the differences. These differences may be attributed to dynamic screening as suggested by Shaviv and Shaviv [75]. Our pair correlation functions match those of Shaviv and Shaviv more closely than those derived from the Debye-H¨ uckel theory, though differences remain between ours and the Shavivs’. In an attempt to account for these differences, we tested their use of a thermostat to maintain constant temperature in their system and their assumption that the long-range Coulomb potential from distant particles is screened sufficiently to be truncated beyond a dis- tance of a few Debye lengths. We found that these two assumptions do not affect the pair correlation functions and are therefore not the source of our disagreement with their results. Thus we have found no problem with the techniques used by Shaviv and Shaviv and have obtained results that are similar to theirs. Though further testing is necessary, our results support the need for corrections to the current screening theory. 115 The path-integral based virial expansions developed by Alastuey and Perez [4] were evaluated for use in solar modeling. By analyzing the limitations of the formal- ism, we determined that it is appropriate for use in much of the Sun. We then adapted the formalism for practical application and used it to develop a solar equation of state. A test of our model in a sample Sun made of pure hydrogen showed that the FK EOS provides the correct results for a mono-atomic plasma. Now that the formalism has proven itself in this test case, it is ready to be extended to a more realistic solar model. This formalism will provide a powerful tool for studying the solar interior. Both of these lines of work improve descriptions of hot, dense plasmas. These results address not only the equilibrium thermodynamics of the solar interior but also our understanding of all stellar matter. By improving the accuracy of solar models, we further our knowledge of nuclear fusion with the hope that it could one day lead to practical fusion applications on Earth. With the coming of the National Ignition Facility, we may soon be able to test this physics in a terrestrial lab. Until then, we rely on the laboratories of the Sun and the stars to explore the nature of fusion and the behavior of hot, dense plasma. 116 Bibliography [1] Abe, R. (1959). Prog. Theor. Phys., 22:213. [2] Alastuey, A., Cornu, F., and Perez, A. (1994). Phys. Rev. E, 49:1077. [3] Alastuey, A., Cornu, F., and Perez, A. (1995). Phys. Rev. E, 51:1725. [4] Alastuey, A. and Perez, A. (1992). Europhys. Lett., 20:19. [5] Alastuey, A. and Perez, A. (1996). Phys. Rev. E, 53:5714. [6] Allen, M.P. and Tildesley, D.J. (1987) Computer Simulation of Liquids. Oxford University Press. [7] Arafune, J. and Fukugita, M. (1992) Prog. Theor. Phys., 87:1467. [8] Asplund, M., Grevesse, N., and Sauval. A.J. (2005). ASPC, 336:25 (AGS05). [9] Assenbaum, H.J., Langanke, K., and Rolfs, C. (1987). Z. Phys. A 327:461. [10] Bahcall, J. N., Brown, L. S., Gruzinov, A, Sawyer, R. F. (2002). A&A, 383:291. [11] Barker, A. A. (1971). J. Chem. Phys., 55:1751. [12] Beatson, Rick and Greengard, Leslie. A Short Course on Fast Multipole Methods. [13] Berrington, K. A. 1997, The Opacity Project, vol. II, Institute of Physics Publish- ing. Bristol [14] Brown, T.M. and Christensen-Dalsgaard, J. (1998). ApJL, 500:195. [15] Carraro, C., Schffer, A.,and Koonin, S. E.( 198). ApJ, 331:565. [16] Christensen-Daalsgaard, J., D¨ appen, W. and the GONG Team (1996). Science 272:1286. [17] D¨ appen, W., Mihalas, D., and Hummer, D. G. (1988). ApJ, 332:261. [18] Debye, P. and H¨ uckel, E. (1923). Z. Phys., 24:305. 117 [19] Deubner, F. L. (1975). Astron. Astrophys., 44:371 – 375. [20] Deutsch, C. (1977). Phys. Lett., 60:317. [21] Deutsch, C., Gombert, M. M., and Minoo, H. (1978). Phys. Lett. 66A:381. [22] Deutsch, C., Gombert, M. M., and Minoo, H. (1979). Phys. Lett. 72A:481. [23] Ebeling, W., Kraeft, W. D., Kremp, D. (1976) Theory of Bound States and Ion- ization Equilibrium in Plasmas and Solids, Berlin, DDR: Akademie Verlag [24] Engstler, S., Krauss, A., Neldner, K., Rolphs, C., Schr¨ oder, U., and Langanke, K. (1988). Phys. Lett. B, 202:197. [25] Ewald, P.P. (1921). Ann. Phys., 64:253-287. [26] Feynman R. P. and Hibbs, A.R. (1965). Quantum Mechanics and Path Integrals Mc Graw Hill, New York. [27] Frenkel, Daan and Smit, Berend (2002) Understanding Molecular Simulation: from Algorithms to Applications. Academic Press. [28] Greengard, L. and Rokhlin, V . (1987). J. Comp. Phys. 73:325-348. [29] Grevesse, N. and Sauval, A. J. (1998). Space Science Review, 85:161. [30] Guzik, J.A. (2006). in Proc. SOHO 18/GONG 2006/HelAs I Workshop, (ESA SP 624, 2006). [31] Hansen, C.J. and Dawaler, S.D. (1994). Stellar Interiors: Physical Principles, Structure, and Evolution, Springer-Verlag. [32] Hansen, J.P. and McDonald, I.R. (1981). Phys Rev A, 23:2041-2059. [33] Hoover, Wm.G. (1985). Phys Rev A, 31:1695. [34] Hockney, R.W. and Eastwood, J.W. (1981). Computer Simulations Using Parti- cles, McGraw-Hill, New York. [35] Hoover, Wm.G. (1991) Computational Statistical Mechanics. Elsevier Science Publishers, Amsterdam. [36] Hummer, D. G. and Mihalas, D. (1988). ApJ, 331:794. [37] Ichimaru, S., Tanaka, S., and Iyetomi, H. (1984). Phys. Tev. A, 29:2033. [38] Iglesias, C.A., Rogers, F.J. and Wilson, B.G. (1987) ApJ, 322:L45. 118 [39] Iglesias, C.A. and Rogers, F.J. (1991). ApJ, 371:408. [40] Iglesias, C.A. and Rogers, F.J. (1993). ApJ, 412:752. [41] Iglesias, C.A. & Rogers, F.J. (1995) ApJ, 443:460. [42] Iglesias, C. A. and Rogers, F. J. (1996). ApJ, 464:943. [43] Kraeft W. D., Kremp D., Ebeling W., R¨ opke G. (1986). Quantum Statistics of Charged Particle Systems Plenum, New York [44] Landau, L. D. and Lifshiftz, E. M. (1980). Statistical Physics. Butterworth- Heinemann. [45] Lavagno, A. and Quarati, P. (2000). Nucl. Phys. B (Proc. Suppl.), 87:209. [46] Lodders, K. (2003). ApJ, 591:1220. [47] Mao, D., Mussack, K. and D¨ appen, W. (2004). “Development of Molecular- Dynamics Tools to Study Hot Dense Coulomb Systems,” in Proc. SOHO 14/GONG 2004 Workshop, (ESA SP 559, October 2004). [48] Martyna, G.J., Klein, M.L., Tuckerman, M. (1992). J. Chem. Phys., 97:2635. [49] Martyna, G.J., Tobias, D. J.,and Klein, M.L. (1994). J. Chem. Phys., 101:4177. [50] Martyna, G.J., Tuckerman, M., Tobias, D. J.,and Klein, M.L. (1996). Mol. Phys., 87:1117. [51] Mayer, J.E. (1950). J. chem. Phys., 18:1426. [52] Meeron, E. (1958). J. Chem. Phys. 28:630. [53] Mihalas, D., D¨ appen, W., and Hummer, D. G. (1988). ApJ, 331:815. [54] Mitler, H. E. (1977). ApJ, 212:513. [55] Mochkovitch, R., and Hernanz, M. (1985). in Nucleosynthesis and Its Implica- tions on Nuclear and Particle Physics, ed. J. Audouze and N. Mathieu (Dordrecht: Reidel), 407. [56] Mussack, K., Mao, D. and D¨ appen, W. (2006). “Molecular-Dynamics Simula- tions of Hot Dense Coulomb Systems,” in Proc. SOHO 18/GONG 2006/HelAs I Workshop, (ESA SP 624, 2006). [57] Nakano, A. (2004). Methods of Computational Physics class notes. 119 [58] Nayfonov, A., D¨ appen, W., Hummer, D. G., and Mihalas, D. M. (1999). ApJ, 526:451. [59] Nos´ e, S. (1984). Mol. Phys., 52:255. [60] Opher, M., and Opher, R. (2000) ApJ, 535:473. [61] Pang, Tao (1997). An Introduction to Computational Physics. Cambridge Uni- versity Press. [62] Perez, A. (1994). “Developpements Diagrammatiques pour un Plasma Quantique dans la Repres´ entation de Feynman-Kac,” [63] Perez, A. and Chabrier, G. (2004). “Quantum statistical corrections on the en- hancement factor in solar fusion reactions”, in Equation-of-State and Phase- Transition Issues in Models of Ordinary Astrophysical Matter, edited by V . Cele- bonovic, W. D¨ appen & D. Gough, AIP Conference Proceedings 731, Melville, New York, 2004, p.208 [64] Perez, A., and D¨ appen, W. (1995). “Rigorous constraints on the ionization of ele- ments in the solar center”, in Proc. of Fourth SOHO Workshop Helioseismology, (ESA-SP-376, June 1995), 15-18. [65] Perez, A., Dppen, W., Mao, D., Mussack, K. (2007). “The Stellar Equation of State with the Path Integral Formalism,” A&A, submitted. [66] Rhodes, E. J., Ulrich, R. K., and Simon, G. W. (1977). ApJ, 218:901 –919. [67] Rogers, F.J. (1986). partition function”, ApJ, 310:723-728. [68] Rogers, F.J. and Nayfonov, A. (2002). helioseismology”, ApJ, 576:1064. [69] Rogers, F.J., Swenson, F.J. and Iglesias, C.A. (1996). Astrophys. J., 456:902. [70] Rolfs, C. (1994). Nucl. Phys. Proc. Suppl. 35:334. [71] Salpeter, E. E. (1954). Australian J. Phys., 7:373. [72] Salpeter, E. E. (1958). Ann. Phys., 5:183. [73] Savchenko, V . I. (1999). Preprint (astro-ph/9904289). [74] Seaton, M. J. (1995). The Opacity Project, volume I. Institue of Physics Publish- ing, Bristol. [75] Shaviv, N. J. and Shaviv, G. (1996). ApJ, 468:433. [76] Shaviv, N. J. and Shaviv, G. (1997). Ap&SS Lib., 214:43. 120 [77] Shaviv, N. J. and Shaviv, G. (2000). ApJ, 529:1054. [78] Shaviv, N. J. and Shaviv, G. (2001). ApJ, 558:925. [79] Shoppa, T.D., Koonin, S.E., Langanke, K., and Seki, R. (1993). Phys. Rev. C, 48:837. [80] Simon, B. (1979). “Functional Integration and Quantum Physics,” Academic, New York. [81] Stix, M. (1989). The Sun. Springer-Verlag, Berlin. [82] Trampedach, R., D¨ appen, W., & Baturin, V .A. (2006) Astrophys. J., 646:560. [83] Tsytovich, V . N. (2000). A&A, 356:L57. [84] Vashishta, P. (2003, 2004). Atomistic Simulation of Materials class notes. [85] Weiss, A., Flaskamp, M., and Tsytovich, V . N. (2001) A&A, 371:1123. [86] Wilson, R. C. and Hudson, H. S. (1988). Nature, 332:810. 121 Appendix: Atomic Units The molecular dynamics simulations were done in atomic units: e= 1 m e = 1 ~= 1 1/(4π 0 )= 1 a 0 = 1 k= 1 c= 1/α≈ 137 so that the reduced units (primed) are: q= e * q’ m= m e * m’ r= a 0 * r’ E= E H * E’ t= ~/E H ∗ t 0 v= α∗ c∗ v 0 T= E H ∗ T 0 122

Abstract (if available)

Abstract In the deep solar interior, matter is mostly fully ionized. Precise helioseismic and neutrino data provide an incredibly accurate description of the Sun. The established MHD and OPAL equations of state do not match the solar data to the available precision. A rigorous quantum-statistical formalism for Coulomb systems is used to develop new tools for equation of state modeling. Although the formalism is in principle a low-density development, the relevant parameters allow the formalism to be applied in the solar center. Solar models can also be improved through the study of nuclear reaction rates. The Coulomb potentials in the ionized plasma are screened, effectively truncating the interactions. This screening makes it easier for ions to tunnel through the potential barrier. In a seminal paper, Salpeter (1954) discussed this enhancement. He based his study on the approximation of a static screening potential. There is a legitimate concern that dynamic effects could alter the result. Since 1996, Shaviv and Shaviv have been examining this question by numerical molecular-dynamics simulations. The calculation is essentially classical, although electrons are treated with effective potentials to mimic some quantum corrections. Shaviv and Shaviv have reported dynamical effects for the fast protons which are the ions that are most likely to engage in nuclear reactions. Given the importance of these effects for solar and stellar modeling, molecular-dynamics simulations have been developed to verify the results of Shaviv and Shaviv in an independent calculation. These studies of hot, dense plasmas will help develop more accurate solar models.

Linked assets

University of Southern California Dissertations and Theses

Conceptually similar

PDF

Improved thermodynamics of the dense solar plasma and molecular-dynamics simulations of the nuclear-reaction rates

PDF

Modeling the hot-dense plasma of the solar interior in and out of thermal equilibrium

PDF

On the Feynman path into the sun

PDF

Laboratory investigations of the near surface plasma field and charging at the lunar terminator

PDF

Quantum molecular dynamics and machine learning for non-equilibrium processes in quantum materials

PDF

Multi-scale quantum dynamics and machine learning for next generation energy applications

PDF

Step‐wise pulling protocols for non-equilibrium dynamics

PDF

Neural network for molecular dynamics simulation and design of 2D materials

PDF

Molecular simulations of water and monovalent ion dynamics in the electroporation of phospholipid bilayers

PDF

High throughput computational framework for synthesis and accelerated discovery of dielectric polymer materials using polarizable reactive molecular dynamics and graph neural networks

PDF

Out-of-equilibrium dynamics of inhomogeneous quantum systems

PDF

Photoexcitation and nonradiative relaxation in molecular systems: methodology, optical properties, and dynamics

PDF

Dynamics of water in nanotubes: liquid below freezing point and ice-like near boiling point

PDF

Shock-induced poration, cholesterol flip-flop and small interfering RNA transfection in a phospholipid membrane: multimillion atom, microsecond molecular dynamics simulations

PDF

The effects of divalent counterions on the formation and stabilization of RNA tertiary structure

PDF

Adiabatic and non-adiabatic molecular dynamics in nanoscale systems: theory and applications

PDF

Experimental and numerical investigations of charging interactions of a dusty surface in space plasma

PDF

The space environment near the Sun and its effects on Parker Solar Probe spacecraft and FIELDS instrumentation

PDF

Advancing ab initio QM/MM free energy calculations: refining, validating and quantifying the reference potential approach

PDF

The effect of microhydration on ionization energy and proton transfer in nucleobases: analysis and method development

Asset Metadata

Creator Mussack, Kathleen Anita (author)

Core Title Molecular dynamics and virial expansions for equilibrium thermodynamics in the solar interior

School College of Letters, Arts and Sciences

Degree Doctor of Philosophy

Degree Program Physics

Publication Date 06/19/2007

Defense Date 05/07/2007

Publisher University of Southern California (original), University of Southern California. Libraries (digital)

Tag equation of state,hot, dense plasma,molecular dynamics,OAI-PMH Harvest,screening,solar model,virial expansions

Language English

Advisor Däppen, Werner (committee chair), Haas, Stephan (committee member), Kalia, Rajiv (committee member), Kunc, Joseph A. (committee member), Warner, Nicholas P. (committee member)

Creator Email mussack@usc.edu

Permanent Link (DOI) https://doi.org/10.25549/usctheses-m534

Unique identifier UC164409

Identifier etd-Mussack-20070619 (filename),usctheses-m40 (legacy collection record id),usctheses-c127-506306 (legacy record id),usctheses-m534 (legacy record id)

Legacy Identifier etd-Mussack-20070619.pdf

Dmrecord 506306

Document Type Dissertation

Rights Mussack, Kathleen Anita

Type texts

Source University of Southern California (contributing entity), University of Southern California Dissertations and Theses (collection)

Repository Name Libraries, University of Southern California

Repository Location Los Angeles, California

Repository Email cisadmin@lib.usc.edu