University of Southern California Dissertations and Theses

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

(USC Thesis Other)

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

PDF

Download a page range

Download transcript

Copy asset link

Request this asset

Transcript (if available)

Content MODELING THE HOT-DENSE PLASMA OF THE SOLAR INTERIOR IN AND OUT OF THERMAL EQUILIBRIUM Copyright 2012 by Hsiao- Hsuan Lin A Dissertation Presented to the FACULTY OF THE USC GRAD UATE SCHOOL UNIVERSITY OF SOUTHERN CALIFORNIA In Partial Fulfillment of the Requirements for the Degree DOCTOR OF PHILOSOPHY (PHYSICS) May 2012 Hsiao-Hsuan Lin Dedication To my parents, sister and husband, 11 Acknowledgments I would like to express deep gratitude to my advisor Werner Dappen, who makes it possible for me to have a wonderful life as a graduate student. As an research advisor he conveyed his wisdom wholeheartedly, and as an role model, I learned from him that there is always a way to see a more beautiful world. I also thank his wife, Veronika, who always being so concerned with my life and my feelings, and with the most sensitive heart. I thank Professors Tu -nan Chang, Aiichino Nakano, Stephan Haas and Joseph Kunc, , for serving my PhD advisory committee and being fri endly and supportive. I also feel gratitude that my parents never doubt my decisions and always provide definite support and endless love. I hope to devote my later life with them as much as I can. Of course I thank Yun g-Ching for being with me all the way, and became my husband one year ago. We share the sweetest memori es during all these years and we will continue for the years to come. I also feel lucky to have two female peers of the same year with me, Yi-Chen and Yaqi, thank you both for sharing your life with me the se years. 111 Table of Contents Dedication 11 Acknowledgments 111 List of Figures v11 Abstract xm Chapter 1: Introduction 1 Chapter 2: Helioseimology 5 2. 1 Basic structure of our sun 5 2.2 The standard solar model 7 2.3 Helioseismology . . 10 2.3.1 History . . . . . 10 2.3.2 Observation . . . 11 2.3.3 Asymptotic behavior of stellar oscillat ions . 14 2.3 .4 Adiabatic Stratification of Convection Zone 17 2.3.5 Inversion . . . . . . . . . . . . . . . . . . 20 Chapter 3: Equation of state 24 3.1 Phy sical Picture . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 24 3.1.1 The Activity Expansion of the Grand Canonical Part ition Func- tion . . . . . . . . . . . . . . . . . . . . . . . . . . . . 24 3.1.2 Mayer's Regroupig and ABE's Nodal expansion method 32 3.1.3 Equation of State With Activity Expansion . 36 3.1.4 Quantum cluster expansion . 41 3.1.5 The OPAL equation of state 49 3.2 Chemical picture . . . . . . . . . . 54 IV 3.2.1 Free-energy-minimization method 3.2.2 The MHD equation of state . 3.2.3 EFF & CEFF . . . . . . . . . . . 54 57 59 Chapter 4: MHD and CEFF EOS and their modifications 61 4. 1 MHD ................................... 61 4.2 Modified MHD . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 64 4.2.1 Planck-Larkin partition function as internal partition function 64 4.2.2 Scattering terms in Free Energy 66 4.2.3 More scattering terms . . . . . . 69 4.3 EFF and CEFF . . . . . . . . . . . . . 71 4.3.1 Approximation to Fermi-Dirac Integrals 71 4.3.2 An Approximate Treatment of Pressure Ionization . 75 4.3.3 Some issues . . . . 77 4.3.4 CEFF . . . . . . . . . . 78 4.4 Deta ils of modified CEFF . . . . 79 4.4.1 Variable transformation . 80 4.4.2 Transformations in first order derivatives . 81 4.4.3 Transformations in second order derivatives 87 4.4.4 Switching off pressure ionization . . . . . . 99 Chapter 5: Results and Discussion 100 5.1 Hydrogen Plasma . . . . . . . . . . . . 100 5.1.1 PLPF and the Scattering terms . 100 5.1.2 MHD and modified MHD . . . 102 5.1.3 CEFF and modified CEFF . . . 104 5.1.4 Compari son between modified MHD and Modified CEFF . . 108 5.2 Hydrogen-Helium Plasma . . . . . . . . . . . . . . 11 0 5.2.1 MHD and modified MHD . . . . . . . . . 11 0 5.2.2 More scattering terms for modified MHD . 111 5.2.3 CEFF and modified CEFF . . . . . . . . . 11 6 5.2.4 More scattering terms for modified CEFF . 118 5.2.5 Compari son between modified CEFF and CEFF(more ) , and modified MHD and MHD(more) . 121 5.3 Outlook .................................. 125 v Chapter 6: The Shear-driven Soret Effect 126 6. 1 The solar tachocline . . 126 6.2 Soret effect . . . . . . . . . . 131 6.3 MD simulation . . . . . . . . 133 6.3.1 Molecular dynamics . 133 6.3.2 Integration method . . 134 6.3.3 Periodic Boundary Conditi ons . 136 6.3.4 Equation of motion . . . . . . . 137 6.3.5 Dimension less unit . . . . . . . 138 6.4 Reverse Non-equilibrium Molecular dyna mics . 139 6.4.1 Method:ca lculating shear viscosity as an example . 140 6.4.2 Methods for thermal conduction and thermal diff usion . 142 6.5 Our RNEMD for sheared mixture system . 6. 6 Program check . . . . . . . 6. 7 Soret effect in gas phase . . . . . . . . 6.8 Soret effect in liquid phase ..... . 6.9 Shear-driven Soret effect in gas phase 6.10 Shear-driven Soret effect in liquid phase 6.11 Outlook ................ . Chapter 7: Conclusion Bibliography . 142 . 143 . 147 . 152 . 155 . 164 . 168 170 172 Vl List of Figures 2. 1 The interior structure of the sun. (Figure provided by J. W. Leibacher) 6 2.2 Power spectrum of velocity observations from the SOI/MDI experi ment on the SOHO spacecraft. The ridges of power concen tration cor respond to separate radial orders, starting at the lowest frequency with the f mode, with n = 0. . . . . . . . . . . . . . . . . . . . . . . . 13 2.3 Reflection, refracti on, and turning points of sound waves inside the sun. (Figure provided by A. G. Kosovichev) . . . . . . . . . . . . . . . 14 2.4 A rough temperature gradient �� profile of the Sun, this doesn't include the complicated information of the near-surface. . . . . . . . . . 19 2.5 Int rinsic difference between 11 obtai ned from an inversion of helio seismological data [2], and 11 of 4 MHD models (M1-4: filled points) and 4 OPAL models (M5-8: empty poin ts), respectively. All results are in the sense Sun - Model. . . . . . . . . . . . . . . . . . . . . . . . 23 3.1 Upper : relation between f- bond and u- bond. Lower: the type of graphs contribute to Debye q- bond. . . . . . . . . . . . . . . . . . . . 33 5.1 Compari son of relative difference (to OPAL model) quant ities between the original MHD (MHD)and the modified MHD with PLPF and extra scattering term s(PL+ex). The result of the modified MHD has been brought closer to OPAL signifi cantly. . . . . . . . . . . . . . . . . . . 103 Vll 5.2 Compari son of relative difference (to OPAL model) quant ities between the CEFF with PLPF and scattering terms(PL+ex), CEFF with only PLPF(PL), and CEFF with only extra scattering terms(ex). . ...... 105 5.3 Compari son of relative difference (to OPAL model) quant ities between the CEFF with PLPF and scattering terms(PL+ex), CEFF with only PLPF(PL),CEFF with only extra scattering terms( ex), and the original CEFF(orig) .. ............................... 106 5.4 Compari son of relative difference (to OPAL model) quant ities between the original MHD (MHD), the modified MHD with PLPF and extra scattering terms(MHD(mod )),original CEFF (CEFF), the modified CEFF with PLPF and extra scattering term s(CEFF(mod)). The result of both the modified MHD and modified CEFF has been brought closer to OPAL signifi cantly, and with almost same good qualities. . . . . . . . 109 5.5 Compari son of relative difference (to OPAL model) of quant ities between the original MHD (MHD) and the modified MHD with PLPF and extra scattering terms(sc+PL). The result of the modified MHD has been brought significa ntly closer to OPAL. . . . . . . . . . . . . . . . . . . 111 5.6 Compari son of relative difference (to OPAL model) quant ities between the MHD with PLPF and scattering terms F5 - F7(MHDex), MHD with PLPF and scattering terms F5 -F8(8), MHD with PLPF and extra scattering terms F5 -F9(9), MHD with PLPF and extra scat tering terms F5 - F10(1 0) and MHD with PLPF and extra scattering terms F5 - F11 (11 ). . . . . . . . . . . . . . . . . . . . . . . . . . . . 11 4 5. 7 Compari son of relative difference (to OPAL model) quant ities between the MHD with PLPF and scattering terms F5 - F7(MHDex), MHD with PLPF and scattering terms F5 -F12(12), MHD with PLPF and extra scattering terms F5 - F13(13), MHD with PLPF and extra scat- tering terms F5 - F14(14). . ....................... 11 5 5.8 Compari son of relative difference (to OPAL model) quant ities between the CEFF with PLPF and scattering term s(CEFF), CEFF with only PLPF(PL),CEFF with only extra scattering terms(SC). . ........ 11 7 Vlll 5.9 Compari son of relative difference (to OPAL model) quant ities between the CEFF with PLPF and scattering terms F5 - F7(ceffex ), CEFF with PLPF and scattering terms F5 - F8 (8), CEFF with PLPF and extra scat tering terms F5 - F9(9), CEFF with PLPF and extra scattering terms F5 - F10(1 0) and CEFF with PLPF and extra scattering terms F5 - F11 (11) .. ................................... 11 9 5.10 Compari son of relative difference (to OPAL model) quant ities between the CEFF with PLPF and scattering terms F5 - F7(ceffex ), CEFF with PLPF and scattering terms F5 - F1 1 (11 ), CEFF with PLPF and extra scattering terms F5 - F12(12), CEFF with PLPF and extra scattering terms F5 -F13(13) and CEFF with PLPF and extra scattering terms F5 - F14 (14). . ............................. 120 5.11 Compari son of relative difference (to OPAL model) quant ities between the modified MHD (MHD(mod)), the further modified MHD with PLPF and more extra scattering term s(MHD(mr)),modified CEFF (CEFF(mod)), the further modified CEFF with PLPF and more extra scattering terms(CEFF(mr)). The result of both the MHD(more) and CEFF(more) has been brought even closer to OP AL, though not signifi cantly. . ............. 122 5.12 Compari son of relative difference (to OPAL model) quant ities between the original MHD (MHD), the further modified MHD with PLPF and more extra scattering terms(MHD(mr)),original CEFF (CEFF), the fur- ther modified CEFF with PLPF and more extra scattering terms(CEFF(mr)). The result of both the MHD(more) and CEFF(more) are now very close to OPAL compared to the original version. . ........... 124 6. 1 Angular velocity profile in the solar interior inferred from heli oseismology. [ 62] In panel(a), the rotational inversion is shown based on the subtrac- tive optimally localized averaging (SOLA) technique, while in (b ),the angular velocity is plotted as a function of radius for several selected latitudes.All inver sions are based on data from the Michel son Doppler Imager (MDI) inst rument aboard the SOHO spacecraft, averaged over 144 days. . ................................ 128 lX 6.2 Simulation box. The reverse non-e quilibrium molecular dyna mics br ings the momentum from the middle slab to the top and bottom slab by exchanging the momentum of the parti cles. Then the momentum flows back through the fluid by fri ct ion, thereby causing a gradient in the velocity in x- direction ........................ 141 6.3 The energy drift of the system with different unphysical heat flux. Nexch is the rate of swapping particles to create the unphysical flux. Nexch = 5 means the particles get swapped every 5 time steps, similar apply to Nexc h=50. We can see that the more frequent swapping the more severe drift of the total energy of the system. . . . . . . . . . . . . 144 6.4 The energy drift of the system with unphysical heat flux and the system with unphysical momentum flux. The total energy(totE -shear) of the system with momentum flux and the corresponding kinetic energy(KE); the total energy(totE -shear) of the system with heat flux and the corre sponding kinetic energy(KE) are shown . The direction of the energy shift in these two system is opposite. . . . . . . . . . . . . . . . . . . . 145 6.5 The test run of the simulation in liquid system, with the the heat flux with only one kind of species(M=M2). We see the temperature gradi- ent was developed but there is no species concentration separation. . . . 147 6.6 The result of Soret diff usion simulation, total step counts is 62000. The density profiles of different masses are still developing. . . . . . . . 149 6. 7 The result of Soret diff usion simulation, total step counts accumulated to 18 4000. The density profiles of different masses are almost satu- rated .. ................................... 150 6.8 The result of Soret diff usion simulation, total step counts accumulated to 482000. The density profiles of different masses are now stable. . .. 151 6.9 The result of Soret diff usion simulation, just compare the density pro file through different time point of the simulation. We noticed that the density profile does not change much since time step counts = 482000 .. ................................. 152 X 6.10 The result of thermal diff usion simulation, total step counts accumu- lated to 1028000. The density profiles has not clearly developed. . ... 154 6.11 The result of Soret diff usion simulation, just compare the density pro file through different time point of the simulation. We noticed that the density profile does not change much since time step counts = 424000 .. ................................. 155 6.12 The result of shear-driven Soret diff usion simulation, total step counts accumulated to 62000. The density profiles has not clearly devel- oped .. ................................... 158 6.13 The result of shear-driven Soret diff usion simulation, total step counts accumulated to 124000. The density profiles has not clearly devel- oped .. ................................... 159 6.14 The result of shear-driven Soret diff usion simulation, total step counts accumulated to 310 000. The density profiles has not clearly devel- oped .. ................................... 160 6.15 The result of shear-driven Soret diff usion simulation, total step counts accumulated to 920000. The density profiles has not clearly devel- oped .. ................................... 161 6.16 The result of shear-driven Soret diff usion simulation, total step counts accumulated to 2706000. The density profiles has not clearly devel- oped .. ................................... 162 6.17 The result of shear-driven Soret diff usion simulation, total step counts accumulated to 1261 0000. The density profiles has not clearly devel- oped .. ................................... 163 6.18 The result of shear-driven Soret diff usion simulation, total step counts accumulated to 36016 000. The density profiles has not clearly devel- oped .. ................................... 164 6.19 The result of thermal diff usion simulation, total step counts accumu- lated to 902000. The density profiles has not clearly developed. . . . . . 165 Xl 6.20 The result of thermal diff usion simulation, total step counts accumu- lated to 1814 000. The density profiles has not clearly developed .. ... 166 6.21 The result of thermal diff usion simulation, total step counts accumu- lated to 2706000. The density profiles has not clearly developed. . . . . 167 6.22 The result of thermal diff usion simulation, total step counts accumu- lated to 11 708000. The density profiles has not clearly developed. . . . 168 Xll Abstract The developments in helioseismology ensure a wealth of studies in solar physics. In par ticular, with the high preci sion of the observat ions of helioseismology, a high-q uality solar model is mandated, since even the tiny deviations between a model and the real Sun can be detected. One crucial ingredient of any solar model is the thermodyna mics of hot-dense plas mas, in parti cular the equation of state. This has motivated effo rts to develop sophisticated theoretical equations of state (EOS). It is important to realize that for the conditions of solar-interior plasmas, there are no terrestrial laboratory exper iments; the only observa tional constraints come from helioseismology. Among the most succes sful EOS is so called OPAL EOS, which is part of the Opacity Project at Livermore. It is based on an activity expansion of the quantum plasma, and realized in the so-called "physical picture ". One of its main competitor is the so called MHD EOS, which is part of the international Opacity Project (OP), a non-classified multi-country cons ortium. The approach of MHD is via the so-called "chemical picture ". Since OPAL is the most accurate equation of state so far, there has been a call for a public-domain version of it. However, the OPAL code remains proprietary, and its "emulation" makes sense. An additional reason for such a project is that the results form OPAL can only be accessed via tabl es generated by the X111 OPAL team. Their users do not have the flexibility to change the chemical composition from their end. The earlier MHD-based OPAL emulator worked well with its modifica tions of the MHD equation of state, which is the Planck-Larkin partiti on function and its cor responding scattering terms. With this modification, MHD can serve as a OPAL emula tor with all the flexibility and access ibility. However, to build a really user-fri endly OPAL emulator one should consider CEFF-based OPAL emulator. CEFF itself is already widely used practical EOS which can be easily implemented within any solar model code. In the present work we have carried the technique of the MHD-based OPAL emulator to the CEFF-based OPAL emulator and succes sfully accompli shed this goal. At the same time, we went beyond the earlier work by adding more terms. In particular, the previous MHD-based OPAL emulator was restricted to the approximation of a hydro gen plasma; our work is extended to the more realistic H-He mixture. In a separate part of the present work, we have exami ned a non-eq uilibrium effect in the solar interior. The effect is located in the zone of the sharp transition between the differential rotation in the convection zone and the solid-sphere rotation in the radiation zone beneath it. This transition was discovered by helioseismology in the 1980s, and the transition zone is called the solar tachocli ne. The tachocline is subject to strong shear and in many theori es of the solar dynamo it plays important role. Being inspired by the well-known Soret effect, which states the mass diff usion drive by a temperature gradient, we have examined if there could also be mass diff usion by a shear flow. If such an effect were to exist, it would have potential applications to the solar tachocline and dynamo. We have run a so-called reverse non-eq uilibrium molecular dynamics (RNEMD) simulation. As a test, we first confirmed the Soret effect by the simulation, then we tested for a shear driven ana logous effect. As a result, we did not see the shear-driven Soret effect in our XIV simulation. We do observe the normal Soret effect due to the temperature gradient cau sed by the numerical scheme we used. XV Chapter 1 Introduction This thesis is mainly about the Sun, our nearest star. With the advantage of being so close to it, we have the chance to study it in much detail, much more than other astronomers can study their objects. The Sun supports life on earth because it is an energy source, and its energy comes from nuclear fu sion react ions within the core. This energy produc es the outward gas pressure that balances the inward-pulling gravitational force. This hydrostatic equilibrium is stable for billions of years. It varies very slowly due to a parametric change in the pressure-density relation, the equation of state, as a result of the change in the chemical composition in the solar core. The chemical composition changes very slowly but gradually due to the conver sion of hydro gen to helium. Right now the Sun has used up about half of its hydro gen supply. It is 4.7 billion years old, and it has another 5-6 billion years to go from now. Thanks to helioseismology, we are able to get informati on about the inside of the Sun by measuring surface pulsations of the Sun. We will introduce the basics of the Sun and describe the metho ds of heli oseismology in Chapter 2. Helioseismology allows for numerical inversions of thermodynamic quant ities ac ross the Sun. The result of these inver sions are the observational data for our work which is a diagnosis of hot-dense plasmas. The most accurate inversion techniques are based on a differential ana lysis between a good solar reference model and the real Sun. The difference between the observed solar oscillation freq uencies and the analogous model frequencies is related by an integral equation to the difference between the interna l-st ructure quant ities of 1 the Sun and the corresponding model data. Present-day solar models (so-ca lled "Standard Solar Models") are so accurate that the linear approximation in the differential analysis is adequate. Due to this linearizat ion, a better the reference solar model leads to smaller mverswn errors. Therefore a good solar model is cru cial. As mentioned ab ove, the equation of state is a key ingredient from micro-phy sics. To meet the high demand helioseismology puts on its quality, various sophisticated equation of states have been proposed in the last 30 years. One is the so-called "MHD" equation of state developed by Dimitri Mihalas, David Hummer, and Werner Dappen [35, 44, 16] as part of the international "Opacity Project". The MHD equation of state postulates the existence of compound states (atoms, ions, mol ecu les) and models the interaction of the surrounding plasma on the se compo unds. It is realized as a free energy-minimization method to find the equilibrium between differ ent ionization and dissociation state s. Another sophisticated formalism is the OPAL equation of state, which is al so part of a large opacity effort, OPAL, realized at the Lawrence Liv ermore National Laboratory [56, 58, 54, 57]. Unlike MHD, OPAL does not assume the existence of any bound species but only deals with the fundamental particles electro ns and nuclei. We discuss the theory of the equation of state in more detail in Chapter 3. In the 1990s it turned out from compari sons with helioseismological observations that the OPAL equation of state is so far the best. MHD is also very good, and it is smoother and more flexible, but OPAL is more accurate. However, there are drawbacks with OP AL. Its results can only be obtai ned from tabl es which have to be generated by the OPAL team. The source of OPAL is proprietary and has not been releas ed, leading to a lacks of flexibility for outside users, in parti cular when they want to change the chemical compositions. This 2 matters since determi ning the chemical composition of the Sun has a very high astrophy s ical importance. To allow a flexible use of the OPAL equation of state, we had started a project of developing a so-called "OPAL emulator" , designed to produce high-q uality OPAL data while residing in the public domain. The first half of Chapter 4 is devoted to the introduction to previous work of my adviser, based on an OPAL emulator realized as a modification of the MHD equation of state. Although that previous work has achieved part of the desired goal, emulating the OPAL equation of state, it has not achieved the equally important flexibility, beca use MHD too is a table-based equation of state. A flexible OPAL emulator has to start from a flexible equation of state. Such a flexible equation-of-state formalism does exist in the form of the so-called CEFF equation of state, which is a descendent of a widely used modeling tool developed in 1973 by Eggleton, Faulkner and Flannery (EFF) [24]. A Coulomb pressure term was added in 1991 by Dappen and Chri stensen-Dalsga ard, the product is called CEFF, which is also discussed in the first half of Chapter 4. The presentation of my own work begins with the second half of Chapter 4. To obtain a really user-fri endly tool , we decided to develop a CEFF-based OPAL emulator. CEFF had to be modified to match the OPAL results accurately. In the second half of Chapter 4 we show these modificati ons on CEFF in detail. Technically, we proceed in two steps. First, the modification used in the prior MHD-based emulator are carri ed over to the CEFF based emulator. As a second step, we have gone beyond that by developing additional terms than tho se used in the original MHD-based OPAL emulator. We show the results of our MHD-based OPAL emulator with new terms , and the results from the CEFF-based OPAL emulator in Chapter 5. In Chapter 6, we present the results of a separate study about the ph ysics of the solar interior. This time it is about an effect 3 outside thermodynamic equilibrium. The key idea is to examine if the well-know solar differential rotation, with its shear, could induce a separation of particle according to their mass. The intuition behind this idea is based on Onsage r's symmetry relati ons, and the hope behind this study was that if such a shear-induced separation of particles exists, it might have an effect on the flow of electrons and nuclei, and thus perhaps explain the electric currents behind the solar magnetism, created by the so-called "solar dynamo". This is a high-risk investigation that I undertook in parallel with developing the OPAL emulator. Such a study can be justified by the potential high importance of its outcome. The location of the Sun where the effect - if it were to exist - would be at work is the solar tachocli ne, a transition zone of differ ential rotation due to the observed fact (from helioseismology) that the inner part of the Sun rotates like a solid but the outer part rotates differentially, with the equatorial part being faster. The transition happens to be at the bottom of the solar convection zone, at a depth of about 30% in radius (200,000 km). The well-known thermodiffusion effect is also called the Soret effect. We have thus been examining if there is an effect ana logous to the Soret effe ct, but driven by a shearing flow and not by a temperature gradient. Our method was a molecular- dynamics (MD) simulation and we have chosen a particular algorithm, reverse non-equilibrium molecular dynamics (RNEMD), which is discussed in detail in Chapter 6, along with the results of our simulations. While our results are still somewhat preliminary, it appears that there is an effect, albeit a second ary effect, in the sense that it is due to a shear- induced temperature gradient. Work remains to be done to see if it has the right quantitative amount to play a role in the solar dynamo. 4 Chapter 2 Helioseimology In this Chapter we introduce the basic facts and simple evolution model of the sun. Then we review the theory of helioseismology method and discuss how can we use our sun as a laboratory to test the equation of state theory with the helioseismology. 2.1 Basic structure of our sun Our Sun is one of the about 1 0 2 3 stars in our unive rse. It is a very normal star, but very special thanks to the exceptional location which allows us to observe it infinitely better than any other star. The Sun is a so-called "main-se quence star" , which means it still fu ses hydrog en into helium. That "burning" (as fu sion is sloppily called by astroph ysicists) is about in half time, which means the Sun is in its middle age. Its own gravity balances the outward pointing forces of gas and radiation pressure. Thu s, the Sun is basically in a hydrostatic equilibrium, which varies only very slowly with the fu sion-induced change of its chemical composition. As far as energy is concerned, the Sun is in a steady state, with energy generated in the core by fu sion continuously moving toward the surface from which it is radiated away into space. This radiation reaches all the planets it and sustains life on Earth (and wherever else in the solar system). The internal structure of the sun is strati fied, like an onion, layer by layer, as shown in Fig.2. 1. The innermost layer is called the solar core, with a density 1.622 x 105kg/ m3i and 5 Figure 2.1: The interior structure of the sun. (Figure provided by J. W. Leibacher) a temperature 1.571 x 107 K. Almost 90% of the solar mass is concentrated within the solar core. The immediate layer outside is the radiation zone, where the heat is transported by radiative diffusion and matter is mechanically stable (no convection). Beyond the radiative zone, diffusion cannot carry all energy (basically because temperature is too low here) and thus energy begins to pile up, triggering a more efficient mode of energy transport by convection. Packets of the plasma are moving up and down vigorously, carrying the heat from the bottom of convection zone to the solar surface. On the top of the convection 6 zone is the so-called photosphere, a relatively thin layer, only about 500 km thick, which is 0.07% of the Sun's radius. The light we observe from the sun is mainly from this photosphere, which justifies its name. Therefore we do not get radiation from the center of the Sun directly; it is scattered many times by the hot, dense plasma, finally reaching the photosphere after a long random walk, which takes several million years, and along which the mean photon energy has decreased from nearly one keV in the center to about 0.5 eV, distributed by a Planck distribution of 5770 K, the temperature of the photosphere. There, the photons finally become free and escape into the dark space. Thus the sunlight we can see and observe is emitted from the photosphere. We can have information about the chemical composition from the spectrum lines, but only for the pho tosphere. Fortunately, as discovered in the 1960s and 1970s, there are acoustic waves going through the entire Sun. Analogous to earthquake waves that reveal the nature of Earth 's interior these acoustic wave allow us to diagnose the Sun's interior: this method is the well-known heli oseismology. Finally, the atmosphere outside the pho tosphere consists of chromosphere and corona, they are not in thermodynamic equilibrium and become at places extremely hot( a few million degrees Kelvin), not because of fu sion but very strong magnetic activiti es. 2.2 The standard solar model The so-called 'S tandard Solar Model' (SSM) is the best physical model of our Sun so far. It requires various branches of physics in order to build a solar model, like thermodynam ics, fluid dynamics, nuclear ph ysics, radiation transfer, convection transport, mixing, ... etc. Thereby, constructing solar models is indeed a very large and complicated enterprise. 7 Here, we describe the main idea in a simplified way that nevertheless retains the basic picture. • The equation desc ribing the hydrostatic equilibrium, reflects the balance of the inward gravitational force and the outward pressure forc e. H dr . .1.b . dP Gmp y orstatlc eqm 1 num: - = -- 2 -, dr r (2. 1) where r is the radius of the Sun, P is pressure, p is density, m is the mass at r, and G the gravitational constant. The mass m within radius r can be determined by (2.2) • Nuclear fusion occurring in the core of the sun, and the energy it generated will transport via different mecha nism within different regions of the sun. In the radiation zone, the heat was carri ed by the radiation, until reaching the convection zone, the heat can be carried by the much faster convect ion. dT 3KpL Energy transport: (radiative transport) (2.3) dr 64 7rar 2 T 3 �� = ( �� ) c (convective transport) (2.4) (2.5) where Ti s the temperatur e, K is opacity, a is the Stefan- Boltzmann constant. ( �;) c is the simplest possible way to express the solar structure within convection zone, therefore one can simply use the adiabatic stratificat ion,(�� ) c = (��)s(�)(��). 8 • While the sun is in a stationary thermal equilibrium, the energy lost from the surface will be balanced by the nuclear energy generated in the core. dL 2 Thermal equilibrium: - = 4 7r r PEnuc dr L(Ro) = L0 (2.6) (2.7) where the L(r) is the luminosity at r (the energy flux going through the shell of radius r), and Enuc is the nuclear energy generation rate per unit mass. The chemical composition will change with the nuclear reaction taking place inside the core ; this drives solar evolution in a parametrical fashion. Additional ingredients have recently been added in the standard solar model, once the importance of diff usion and gravitational settling of elements was confirmed by the evidence from helioseimology. Therefore, these diff usion processes are now generally included in standard model. To be more specific, the rate of change of the hydro gen abundance X is written[9]: (2.8) here RH is the rate of change in the hydrog en abu ndance from nuclear reacti ons, D H is the diff usion coefficient and VH is the settling speed. And there are similar equations for other elements. The term DH8Xj8r tends to smooth the hydrogen distribution, while VH X will lead to separat ion, due to the different mass of elements. The lighter ones, like hydrog en, will rise towards the surface while others sinks down. 9 Since these general equations are true for all stars, what then makes every star differ ent from each other? The variety of stars is first determined by the star's mass: two stars of the same mass go roughly through a similar scenario of stellar evolution. But even with the same mass, stars can be quite different beca use of the properties of the stellar material, which are opacity, equation of state and nuclear reaction rate. They are the so called "const ituent equations", usually headed under "microphy sics", while the equations above are labeled under "macrophy sics". Among the three microphysical properties, the equation of state plays the dominant role since it is also a key ingredient of the other two quantit ies. 2.3 Helioseismology 2.3.1 History The first definitive observations of oscillations on the solar surface were made by Leighton et al., who detected periodic oscillations in the local Doppler velocity of roughly 5-min periods. These oscillati ons had quickly drawn people's attention. The breakthrough obser vation was made by Deubner (1975), in which for the first time ridges in the wavenumber frequency diagram were identi fied, thus revealing the global-mode nature of the oscilla tions. Similar results had also been reported by Rhodes et al. (1977), who then further compared the observed freq uencies with computed models, allowing them observationally to find the depth of the solar convection zone. This a ana lysis of the frequencies of these high-degree acoustic modes provided the first solid helioseismic infe rences. Compari sons of the observed and computed frequencies indicated that the solar convection zone was 10 deeper than previously assumed [5 1, 65] and showed that the frequencies were sensitive to details of the equation of state of the matter in the Sun.[64, 40, 3] 2.3.2 Observation To observe the osc illation pattern, one should measure the a series of velocity on the solar surface within a range of time. The technique employed is to measure the Doppler shift of a given line in the solar spectrum. The main idea is to focus on the two symmetric side band of a certain selected spectrum line. If the center of the spectrum shifted due to the motion of the source(Doppler eff ect) , then the intensity of the side band we were focu sed on will then change accordingly, thereby the velocity of the source can be detected. A common feature to all experiments is the need for long time series to obtain an adequate signal-to-no ise ratio and frequency resolution. Also, it is highly desirable to have few or no time gaps in the data, to avoid sidebands in the power spectra. This can be achieved from the ground by combi ning data from several observing locations (e.g. GONG, BISON), or from space (by the heli oseismological instruments MDI, VIRGO and GOLF on the SOHO spacecraf t). The snapshots taken of the Sun are actually velocity images, where each of them is in fact a superposition of thousands of osc illation modes. An amazingly high spatial reso lution can be obtained thanks to the Sun's proximity to earth and careful choice of stable wavelength detection instruments. To fully utilize the hel ioseismic data to inquire inside the Sun, what we need is the power-frequency spectrum ( Fig.2.2) for each normal mode specified by a unique set of eigenvalu es n, l, and m. In order to obtain the power spectrum, the following subsequent procedures have to be taken step by step: 11 • First, the Doppler shifts of spectrum lines are observed, which easily leads to a velocity image of the solar plate. After subtracting the Sun's self-rotating velocity, the residue in temporal domain VD(B, ¢, t) is obtained and ready to be processed. n,l,m (2.9) • Since our ultimate goal is to single out each individual normal mode, we take ad van- tage of the orthogonality property of spherical harmonic by applying a suitable weight function Wzomo ( B, ¢) to the signal VD ( B, ¢, t) and integrating. What we get is a time series of velocity Vto mo ( t) with reference degree l 0m0. Obviously we would like terms in \110 m 0 ( t) to be as close to bu0b mmo as possible. \llo mo = 1 VD(B , ¢, t)Ww mo( B , ¢)dA = L S lOmO lm A nlm ( T) COS [ Wntm ( T) t + f3n lm,l 0 m0 ( T)] (2. 10) n,l,m • Next, Fourier transform is applied to the time stream of Vto mo ( t), and a power spec trum associated with the reference l 0m0 is obta ined. Collectively doing this for various degrees yields so-called l - v diagram, which is illustrated in Fig. 2.2. Apparently in order to obtain sharp spectrum lines, we would like the time series of \ll0 m 0 ( t) to be as long and continuous as possible, which is exactly the original 12 motivation for constructing world-wide multi-site solar obsetvatories, so that the night gaps can be avoided after data are patched together. 8 6 - N :::c s 4 2 0 0 50 100 150 l 200 250 300 Figure 2.2: Power spectrum of velocity obsetvations from the SOI/MDI experiment on the SOHO spacecraft. Thelidges of power concentration correspond to separate radial orders, starting at the lowest frequency with the f mode, with n = 0. • Finally, for one l0, the frequencies and other parameters of each individual mode can be obtained from its spectrum 1/l0(v) by numerical fitting. This is usually the most difficult step and involves lots of complication due to the stochastic nature of excitation and the data's incapability of covering the whole solar surface. 13 2.3.3 Asymptotic behavior of stellar oscillations Figure 2.3: Reflection, refraction, and turning points of sound waves inside the sun. (Fig ure provided by A. G. Kosovichev) If we ignored the fact that our sun is a little ellipsoidal instead of a pe rfect sphere, the spherical harmonic is still a good approximation to represent normal mode. (True normal modes of the ellipsoid would have a main spherical analog with a little in-mixing from a few adjacent spherical modes). The three ch aracteristic wave-numbers that follow this choice would be n, which is the ra di al order; l, which is the degree; and m, which is the azimuthal order[l 0]. The displacement of a normal mode can be written as, 14 where ar,a8 and a<P are unit vectors in the r, () and ¢ direct ions. y; m is a spherical harmo nic. � r (r) is the eigenfunction that determines the radial displacement, while �h(r) is the one for displacement on the spherical surface. The horizontal wavenumber kh is related to the degree of the mode by, (2. 12) where L = jl(l + 1). Had there not been attenuations and damping, the frequency w would be real. For an object that is not spherical, w would depend on all three wave numbers Wnl m · But since the Sun can be regarded as almost perfectly spherically symmet- ric, there is no preferred axis as long as it passes through the center, so as a simplificat ion, the dependence of w on m can be relaxed. Most of the modes observed in the Sun are acoustic modes, often in high orders. In this case, we can study their asymptotic behavior (ray approximat ion,[28], Fig. 2.3) by approximating the modes locally as plane waves, etkx-twt, with the dispersion relat ion: (2. 13) where k = k r e r + kh is the wave vector, with k r and kh being its radial and horizontal components. The adiabatic sound speed c(r) describes the properti es of the modes. With the expre ssion of kh in Eq.(2.12), we can have k r written as: 15 (2. 14) this equation can be interpreted as the behavior of a ray of sound (see Fig. 2.3). While the temperature increases with decreasing radius r, the sound speed gets larger. It can be seen from Eq.(2.14) that while c incr eases kr decr eases, which means that the ray is losing its ability to go deeper. Eventually, the ray of sound will be refracted while the sound speed get too high. This happe ns at the so-called turning-point, where kr reduced to zero : w L (2. 15) L = Jl(l + 1), and we see that while l is smaller, It is bigger, therefore the low l modes can penetrate the deeper part of the sun. Modes with l = 0 reach the very center, while the highest observed mode, with about l rv 1000, are trapped near the surface within a small zone of depth of a few percent of the solar radius. The inner structure of the star can only influence the low l modes, while the outer envelope matters to modes of all degrees. We also noted that if two modes with differ ent l degree but same frequency, then they will have similar seismic behavior except near the turning point. It is beca use that the first term in Eq.(2.14) dominates over second term except near turning point. The requirement of a standing wave in the radial direction implies that the kr must have a phase shift of the multiple of 1r, while traveling from r = R to r = r t. (2. 16) substituting kr with Eq. 2. 14 leads to, 16 (2. 17) They may also be written as: (2. 18) where (2. 19) The observed frequencies of solar oscillation satisfy the simple functional relation of Eq. (2.18), which was first found by the solar observer Duva11 [19]; it is commonly referred to as the Duvall law. The importance of the Duvall law is that the absolute n can be inferred from the fact that only the correct identification of n results in a single-line plot of all values of Wnl , and such a correct identification is indeed possible.[18] 2.3.4 Adiabatic Stratification of Convection Zone The energy generated from the nuclear reaction at the solar core is first transported out wards by radiation. This transfer happe ns in a diff usion proc ess, which is relative slow. Basically it is enabled by the temperature gradient, which causes an anis otropy of local isotropic radiators that follow Stefan 's law for the flux F = a-T 4 • Over a temperature drop of dT, there is an anisotropic flux component proportional to T3dT. Since temperature decreases from the center to the surface, this anisotropic flux component points outward. 17 It is the radiative energy transport. Since this energy throughput decre ases with T 3 , at some point this mecha nism cannot carry the entire energy arriving from below. At that point, convection sets on. Now, convection can transport heat much more efficiently than radiation (it has virtually infinite conductivity ). Ty pical convective turn-over times are a month or less, which is much shorter than the cooling times for the pieces of solar matter undergoing convection. As a result, convection is largely adiabati c, and the net tempera ture profile in the convection zone if that of an adiabat aT 11 - 1 T - > --- 8P - 11 P (2.20) Tr ansition from radiative to convective happens at a location defined by the so-called "Schwarzschild criter ion", where the absolute magnitude of the radiative temperature gra- dient becomes equal to the adiabatic gradient of the solar plasma, see Fig.2.4 for a com parison. Although the actual convection proce ss in the Sun is actually very complicate d, the assumption of an adiabatic strati fication of the convection zone is mostly sufficient to yield good result s, even if it is indeed a rather crude way of treating the complexity of convective motion. Numerical simulations of convection have shown the validity of the se crude reci pes. Despite all the complicat ions that are inherent to convection, the structure of convection zone of the Sun is relatively simple. We are fairly certain that except at the top thin layer, the rest of the convection zone is almost nearly adiabati c. At adiabati c condit ions, pressure is related to density by specific heat ratio at fixed entropy, a thermodynamic quantity called adiabatic exponent and defined as 11, 18 0. 0. 0. 0. Core d(lnT) d(lnP) Radiation Zone Convection Zone /Adiabatic ·-·-·-·-·-·�·-·-·-·-·-·-·-·-·-·�;·--------- 0.2 0.4 0.6 0.8 1.0 R Figure 2.4: A rough temperature gradient �:;:� profile of the Sun. this doesn't include the complicated information of the near-surf ace. _ (8lnp) /l- 81 np s (2.21) Therefore this equation, together with the hydrostatic equilibrium equation pretty much determines the seismic structure of the adiabatic convection zone. 8p 8r Gmp r 2 8m 2 - = - 4 7lT p 8r (2.22) (2.23) 19 Since the structure of the convection zone is largely determined by the adiabatic tem perature gradient, it must be insen sitive to opacity. Also, the convection zone is too cold to allow any nuclear fu sion. Therefore, the equation of state is the one and only ingredient at the "micro-physical" end that determines the structure of convection zone. Therefore, the convection zone provides us an ideal environment to disentangle the equation of state from the complicati ons that are entangled in other parts of the Sun. We can thus study the thermodynamics of solar material under unique solar conditions of high temperature and high density, conditions that have so far been unattainable on Earth. 2.3.5 Inversion When the Sun vibrates with a normal mode, it causes periodic oscillation on the solar sur face and at the same time provides informati on from the deep inside alone the path they have traveled to the surf ace. This allows people to decipher the mysterious interior of the Sun. In principl e, if one can have all the data from those frequencies, it is possible to reconst ruct the internal structure of the sun, and it is exactly the reason that make peo ple put the effort into it and carry out so-called "solar inve rsions". In fact, some of the techniques that were applied to helioseismology were already practi ced by seismologists to probe the inside of the Earth. Here is a re-asses sment of the goal of the inve rsion: to isolate informati ons of the Sun at localized positions. ([25, 12]) Asymptotic Inversion Asym ptotic inversion was first suggested by Gough [26, 27] and then applied to solar oscillations by Chri stensen-Dals gaard [11]. In compa rison with diff erential inversions (see below), there are three thi ngs that are distinct about asymptotic inve rsions: 20 • Asym ptotic inversion is not differenti al , it doesn't need a reference model. • Asym ptotic inversion comes from the asymptotic properties of p modes, so the errors in the properties will also propagate into the inversion results. • Application to observed frequencies by Chri stensen-Da lgaard et al. [11] yields an accuracy of sound speed of better than 1 %. From previous Eq. (2.19), the left hand side, F(w;_t ) is an observed quantity, while localized sound speed c (r) is on the right hand side, so inverting this equation analytically by diffe rentiation yields information regarding local sound speed r( r ::_) = RG exp [-� 1 c /r ( w- 2 - r 2 ) - 1 / 2 dF clJJJ] ' r 7f Cs/ Rev c 2 clJJJ (2.24) The difference between real mode properties and their asymptotic approximation lim- its the accuracy of this approach. A higher accuracy of inversions can be obtai ned by numerical (diff erential) inve rsions, which at the same time serve as a good tool to test the quality of asymptotic inver sions. Differential Inversion If we have a solar model which is good enough, denoted as mod, then the difference between the model and the real sun,denoted as obs, ought to be small . With this small diffe rence, a linearized theory may be used: (2.25) 21 where the so-called "kernels" K:; ,t(r) and K�I(r) are determined from the eigenfunc tions of reference model. (In the spirit of perturbati on theory, the kernels can be used from the unperturbed reference model). An additional constraint is that the mass of the Sun and the reference model has to be same, i.e., (2.26) After applying standard inversion technique s, localized informati on of the set of the chosen variables (in the current case, density p(r) and adiabatic exponent ''(1 (r)) is extracted. This can, in turn, be used to improve the reference model after carefully study ing the discrepancies and feed the modified model back into inversion. The ultimate goal of the iteration would be a "perfect" inverted structure of the Sun. For our present study, in the fo llowing , such state-of- the-art inversions are considered the observational data which we confront with our theoretical work. As an example, Fig. 2.5 shows the difference of / l between the models and observa tion. This provides an excellent tool for selecting a better equation of state for the Sun, and at the same time guides the path towards improvements. In particular, in order to obtain a fine-tuned solar model that matches the observation, a range of different models have to be tested, with different chemical composition, different surface layer, etc. In such an analysis, an accurate and versatile equation of state is needed. No such fully adequate equation of state exists so far; its development is the main purpose of our present research. 22 6. 136 0.002 0 � !J � � r. '...... � r. 'Q -0 002 � • • • • • • • • .. .. .. G-- B 0 0 0 D D D -0 004 [', [', 6 0.8 5.983 Ml M2 M3 M4 M5 M6 M7 MB 0.8 5 Log (T) 5.777 5.436 0.95 3.762 D i I I I I / t t 1 Figure 2.5 : Intrinsic difference between 11 obtained from an inversion of heli oseismolog ical data [2]. and 11 of 4 MHD models (Ml-4: filled points) and 4 OPAL models (M5-8: empty points). respectively. All results are in the sense Sun - Model. 23 Chapter 3 Equation of state In this chapter we review the theory of the equation of state both in physical picture and in chemical picture. We then discuss the OPAL equation of state which is based on physical picture in more details, and list some of the equation of states in chemical picture, i.e. MHD EOS, CEFF EOS. The more detailed discussion about tho se equation of states in chemical picture will be in the chapter 4. 3.1 Physical Picture 3.1.1 The Activity Expansion of the Grand Canonical Partition Func- tion If a system can be treated classi cally, then the Hamiltonian for N particles has the form: N 2 � Pi � H = L....t 2m + L....t U ij i = l i <j (3 .1) where Pi is the momentum of the ith parti cle and uij = u( [ ri - rj [) is the potential energy of interaction between the ith and the jth particle. If the volume of the system is V, the canonical partiti on function is 24 QN(N, T) = N ! � 3 N J d 3 N pd 3 N r exp( �/3 L ;! � f3 I:: uij) (3.2) � z<J After integrating over momenta, it becomes 1 J 3 N � Qr(N,T) = N ! h3 N d rexp(�f3 L..t uij), t<J (3.3) where >- = J 21rn? /mkT is the termal de Broglie wavelength. Qr is called the config uration intergral. Then the canonical partition function can be written as (3.4) and in this notation the grand canonical partition function is o( NT)= � ( ef3f.k ) N Qr( V , T) ...... z, ' L....t ).,3 N ! . N =O (3.5) We then denote e;; as the activity z. The interaction function is introduced by (3.6) and hence we can write 25 = (!1 2 + 1)(!1 3 + 1) * ... = 1 + 2:: 1ij + 2:: 2:: 1 ij 1 i ' j ' + ... N ?_ i>?. 1 Therefore the configuration integral can be represented by Qr (V, T) = J d 3 r1 ... d 3 rN n (1 + f ij) z<J = J d 3 r1 ... d 3 rN[1 + (!1 2 + !1 3 + ... ) + (!12!1 3 + !12!1 4 + ... ) + ... ] To enumerate all the terms in the above expa nsion, a graph method has been developed to associate each term in the integral . Define a cluster integral b1(V , T) by 1 . b z (V, T) - 1 V (sum of all poss1ble l - clusters) , l . (3.7) where l-cluster means a cluster made with l particles, and between each particle, we use a f-bond (defined in Eq.(3.6)) to connect them. For exampl e, the 1-c luster and 2-cluster are trivial, beca use there is only one possible graph for both of them: (3.8) 26 (3.9) However, for the 3-clust er, we can see that there are actually 4 possible graphs: J k J (b) (c) (d) The type of graphs (a), (b), (c) are called singly-connected graphs and (d) is called irreducible graph. Once a f -bond is broken, a singly-connected graph will become two separate cluste rs, but an irreducible graph will stay connected. Therefore we can see four terms in the cluster integral of b3 (3.10) and the general form of b1 can be written as (3 .11) Any N-particle graph with m1 l- cluster shall follow this constraint N L lmz = N. (3.12) 1 =1 27 With this definition of cluster integra ls, one can deduce that the configuration integral is in relation to b1 by N - 'II- 1 ( )mt Qr - N. ! Vb z . 1 =1 mz . The Canonical partition function is The grand-canonical partition function is or, equivalently, one can write: from which the cluster expansion of the equation of state can be obta ined. { p "'"'00 b l kT = L.. l=l zZ (3.13) (3.14) (3.15) (3.16) (3.17) (3.18) If our system is a dilute gas, we may expand the pressure through the power of N jV, and we obtain the virial expansion also called the density expansion. Within this scope, we may take the limit of the equation of state to be 28 { where p """'00 b ' l kT = L... i=1 l z 1 """'00 lb ' l v = L... /=1 t Z The virial expansion of the equation of state is defined to be (3.19) (3.20) (3.21) (3.22) and a1(T) is called the l-th virial coef ficient. We can get the relation between virial coef ficients and cluster integra ls, by putting equation (3.19) and (3.22) all together, and require they to match in each order of z. 00 00 """'00 b ' l "'""' a ( "'""' nb ' z n ) l - 1 = L... i=1 l z L.....t l L.....t n """' oo lb ' z l 1=1 n=1 L... /=1 l This is equivalent to the condition (3.23) (b;z I 2b;z 2 I 3b;z3 I ... ) [a 1 I a,(�� nb�z n ) I a 3 (t. nb�z n ) 2 I .. · ] � b;z I 2b;z 2 13 b;z3 1 ... (3.24) By matching the coef ficients of each power of z, one can get 29 b� = 1 �b' 2 4b� 2 � 2b3 �20b� 3 + 18b�b� � 3b� In order not to be messy, we will denote b� as b1 hereafter. (3.25) (3.26) Up to this point, one can also define the irreducible integral /31 in a similar way via the cluster integrals b1 : f3l + J · · · J L IJ fijdri ... drf l.V l:Si<j:S ( l+l ) irreducible bonds (3.27) In this case, the /31 are associated with the irreducible graphs only, and the relation between the cluster integral s and the irreducible integrals can then be deduced. For example, let's take a look of b3 again: (3.28) and /3 2 is (3.29) together with 30 the relation between b3, /31 and /32 can be expressed like: and the general form can be proved to be where the set { nk} satisfies l - l L knk = l-1 . k =l At the same time the relation between a 1 and /3 1 can also be deduc ed: Now define Mayer's S function as 1 l + 1 - - l -f3t-l (l � 2) with which the equation of state can be written down as (3.30) (3.31) (3.32) (3.33) (3.34) (3.35) 31 P as -- =n+S - n kBT on (3.36) Those are just the general form of the expre ssions. Once we consider particular poten- tial, for instance, the Coulomb potential, then we will have to deal with some divergences of those coefficient s. 3.1.2 Mayer 's Regroupig and ABE 's Nodal expansion method Consider the second virial coefficient with Coulomb potentials, u(r) = zi: je 2 : (3.37) The first divergence in (3.37) can be removed by a neutralizing background, hence the leading divergence will be the term of ( u(r)jkT) 2 To follow Mayer's procedure, we should introduce u-bonds in terms off -bonds (3.38) In Fig.3.1 one can see that the interaction fi j, denoted by the f -bond, can be replaced by the summation of infinite products of two-body interaction u(rij), denoted as u-bond. 32 f-bond + 0 + 0 +0 + 0 + + + u-bond + + ... = Debye q-bond Figure 3.1: Upper: relation between f- bond and u- bond. Lower: the type of graphs contribute to Debye q- bond. Only one f -bond is allowed to connect the point pairs, however there can be arbitrary number of u-bonds. In this regrouping procedure, we only consider the u-bond interac- tions. Including all the possible distribution of cluste rs, we can regroup them according to their proto-graph, namely the form of the graph. So now the summation is not within certain cluster but among all kind of possible cluste rs. 33 Mayer demonstrated that after summing up the contributi on of the graph type in Fig. 3.1, the S function will be equivalent to Debye screened Coulomb potential: (3.39) 2:::' means the summation excludes the ring type graph, since this type of graphs is differ ent from others and can be dealt with analytically, and S · = 1 r m g 12 7r X b ' ( k : T ) rin g = n � 24�A b ' which is the same as the well-known Debye-Hiickellimit law. (3.40) (3.41) Mayer succes sfully showed that the Coulomb u-bonds can be replaced by Debye q- bonds, therefore the long-range divergence caused by the ( u(r)/kT) 2 term can then be removed. However, there is still a divergence occurring not at r = oo but at r = 0. Let's take a look of the 2-parti cle clusters in the S function: (3.42) 34 In the k-sum mati on, k vanishes beca use of the neutralizing backgroun d, and k = 2 was included in the B ring term. The terms of k � 3 are still divergent at r = 0. So what Abe proposed is to sum over all the prototypes of 3, 4, 5, ... bonds, and the result is (3.43) (3.44) In this way, the divergence can also be removed. The general expre ssion of w 3 (r) can be found in Abe's paper. We can see that the (3.43) is very similar to the 2nd virial coef ficient except for the term n 2 (2 nfJXiJ - �/3 2 AD). This happens generally to all ana logous term in the virial expansion and the expansion considered here. To summarize the results, we have : with S defined as: P as equation of state : -- = n + S-n --;::;- , kBT un From Abe's nodal expan sion, Scan be expre ssed as 00 S = B ring + L S m m= 2 (3.45) (3.46) (3.47) 35 where B ring and 8 2 are : and a 2 is 2nd virial coefficient with the Screened Coulomb potential q(r), 3.1.3 Equation of State With Activity Expansion (3.48) (3.49) (3.50) In principl e, the virial (or density) expansion and the activity expansion are both rigorous for an interacting plasma or gas. They will both give the same exact result if one can manage to include an in finite number of terms. Finite truncat ions obviously lose this rigor and the decision which of the two one should use will depend on how efficient and how accurate each method is in a given application. In the case of an interacting and reacting plasma, where bound species form and disintegrate, and where there is an equilibrium state, we can see that the activity expansion is superior to the density expansion. We might use a simple demonstration via the example of the reaction of the formation of the hydrog en molecul e: (3.51) Since here only 2-body interaction was considered, we can write down the equation of state in terms of activity z, to the second order. 36 (3.52) and the density can be duduced by: (3.53) By solving eqautions (3.52) and (3.53), we can have z H be expressed by PH (3.54) Now let's investigate Eq. (3.54) in two limits of temperatu re. First, when temperature is high, bHH --+ 0, from Eq. (3.53) we have PH = ZH, then the equation of state (3.52) can be written as: p 2 b 2 kBT = PH- P HH = PH+ P aHH (3.55) where the virial expansion is recovered. At high temperature, if the density is low enough to be ignored, then activity expansion and density expansion performs equally well. At low temperature limit, bHH is exponentially large, hence we have ZH --+ J PH/2 b HH « PH, so the Eq. (3.52) become s: (3.56) 37 The above relation tells us that at low temperature limit, hydrog en molecules form. To be able to recover this result, in terms of density expansion, we should realize that the 2nd order activity expansion coefficient should be kept. Hence from the relation between cluster expansion and density, a l b l = 1 a 2 - b 2 a 3 4b� - 2b 3 a 4 -20 b� + 18b 2 b 3 - 3b 4 (3.57) (3.58) We notice b 2 keeps appearing in each order of the virial coef ficients, so in principle one will have to consider all the virial coef ficients in order to recover the molecule formation. In practice only including a 2 will certainly not do this job. With this informati on in mind, one can see activity expansion is actually a good approach dealing with the interacting plasma system. Hereafter let us discuss the equation of state with activity expansion. As seen in the previous sections, the grand canonical partition function B(T, V, J-l) can be written in terms of activit ies: and z is activity such that: ln B L oo b l -- = zz v l=l (3.59) (3.60) 38 (2s+l) is the spin degeneracy, A is the de Broglie wavelength, A = h jv21rmkBT. From the relation between virial coefficients and irreducible integrals Eq. (3.34), we can define the similar S fu nction, however in activity expansion: S f ___f!_j_ z l+l = - f � z l+l . 1=1 l + 1 1=2 l + 1 And the equation of state can be shown in [55]: P 2:: 1 S 2:: 00 Z ( 0 ) m - 2 ( 8S ) m - - = b1z = z + + - -z - kBT m ! 8z 8z 1=1 m=2 Follow similar procedure as Abe's nodal expansion, S can be reduced to (3.61) (3.62) (3.63) (3.64) Cons equently, the equation of state of one component plasma(OCP) in activity expansion can be reordered with order of z: P 2:: 00 OXR 8C m - - =z+ xR + (C m + z - - - ) �T & & m=2 (3.65) (3.66) 39 This procedure can be generalized to the two component plasma(TCP) gas (species-a and species-b), and the equation of state is p � Za a m - 2 k T = Z a + Zb + S + L....t - 1 ( � Z a ) B m. u� m=2 � Zb a m - 2 as m as as a 2 S + L....t - 1 ( �) ( --;::;--- ) + Z a Zb � -;::;- a a + ... , m=2 m. UZb UZb UZ a UZb Z a Zb And the generalized S is (3.67) (3.68) Z a , zb, e a , eb are the activity of species a, b and the charge of species a, b, respectively. s 2 = 2: s 2 ,ij = s 2 ,aa + s 2 ,a b + s 2 , bb i,j=a, b 2 'iT 2 S 2 · · = - a 2 · · + z ·z · [ 2 1T' (/3 e ·e · )A - -( f3 e·e · ) AD] ,ZJ ,ZJ z J z J D 2 z J ' (3.69) (3.70) (3.71) Agagin, the equation of state can be reordered to have orders in z, and just shown to the order of z51 2: 40 (3.72) 3.1.4 Quantum cluster expansion This section follows the discussion in [34]. Kahn and Uhlenbeck develop a cluster expan- sion in quantum statistical mechanics. Consider N identical particles enclosed in a volume V. Then the Hamiltonian can be written as, (3.73) The canonical partition function therefore has this form: QN(V , T) = Tr e - (3 H = J d 3 N r L w:(l, ... , N) e - (3 H W a(l, ... , N) a (3.74) where {\[! a } is a complete set of orthonormal wave functions appropriate to the system considered, and the set of coordinates {r1 , ... , rN} is denoted as {1, ... , N}. The sym metry of the wave function will depend on the statistics of the particles in discussion. Let us define 41 WN (1, ... , N) N ! >.. 3 N 2::: w:(l, ... , N) e- rm wa(1, ... , N) (3.75) Then Q N (V, T) can be written in terms of W N ( 1, ... , N), 1 J 3 N QN(V, T) = N ! >..3 N d rWN(1, ... , N) We list some of the properties of the function W N ( 1, ... , N) without proof: • W N(1, ... , N) is symmetric function of its arguments. (3.76) • W N ( 1, ... , N) is invariant under a unitary transformation of the complete set of wave functions {wa}· If we define a function U 2 (1, 2) by W 2 (1, 2) = W1 (1)W 2 (2) + U 2 (1, 2), we expect that, as lr1 - r 2 l -+co , (3.77) Hence the integral of U 2 (1, 2) over r1 and r 2 should be the analog of the 2-cluster in classical statistical mechanics. 42 This definition may go on for larger cluster, and we have W1(1) = U1 (1) = 1 W 2 (1, 2) = U1(1)U1(2) + U 2 (1, 2) W 3 (1, 2,3) = U1(1)U1(2)U1(3) + U1 (1)U 2 (2, 3) where m1 is zero or a positive integer and { m1} satisfies this condition: N L lmi = N, l=l (3.78) (3.79) (3.80) (3.82) the sum over p are the distinct ways of permutation of the N particles ( 1, ... , N). We can use the relations Eq. (3.78) -Eq. (3.80) to find the inverted relations: U1 (1) = W1(1) = 1 U 2 (1, 2) = -W1(1)W1(2) + W 2 (1,2) U 3 (1, 2,3) = 2W1(1)W1(2)W1(3) - W1(1)W 2 (2, 3) (3.83) (3.84) (3.85) (3.86) 43 With these expre ssions, the l-cluster integral b1 (V, T) can be defined by (3.87) To have the partition function Therefore the canonical partition function is given by (3.89) which resembles the classical partition function. Now we want to focus on the main diffe rences between the quantum cluster integrals and the classical ones. For an ideal quantum gas, we have (ideal Bose gas) (ideal Fermi gas) 44 Unlike the classical cluster integrals, b1 of the quantum ideal gas would not vanish even when l > 1. This is not due to the particle interaction but mainly from the quantum effects depending on the nature of the particles. To calculate the second virial coefficient a 2 for any system it is suffi cient to calculate b 2 , since a 2 = - b 2 • To find b 2 we need to know W 2 (1, 2), which is a property of the two-body system. The Hamiltonian for the 2-body system in question is (3.90) with eigenfunctions Wa (l, 2), and eigenvalues E a , and the Schrodinger equation: (3.91) Let (3.92) Then We have (3.93) (3.94) 45 The quantum number a refers ti the set of quantum numbers (P ,n). The relative wave function 1/J n (r ) sat isfies the eigen value equation (3.95) which normalized to: (3.96) Use Eq.(3.93) as the wave function while calculating W 2 (1, 2), W 2 (l, 2) = 2A6 L [W a(l, 2) [ 2 e �fE a = � L L e �f p 2 / 4m e �fc (3.97) a P n In the limit as V ----+ oo the sum over P can be evaluated: (3.98) where A= J2 1rn? jmkT, the thermal wavelength. (3.99) n Repeat the same calculation for the noninteracting particles, we obtain (3.10 0) n 46 the superscript ( O ) refers to the non-interacting, (ideal) system. From Eq.(3.84), Eq.(3.87) one can get: Hence b - b ( o ) = - 1 J d3 Rd3r[W (1 2) - VV: ( 0) (1 2)] 2 2 2V 2 , 2 , = 2 V2 A3 J d3r L [11P n (r)l 2 e-En-I1P �0 )(r) l 2 e-E �,0)] n n For the noninteracting system, E � o ) forms continuum. (3.1 01) (3.10 2) (3.103) For interacting systems the spectrum of E n in general contains a discrete set of values EB, corresponding to 2-body bound sta tes, and a continuum. Let the density of states of an interacting system be g( k) and the one of a noninteracting system be g(0) ( k). Thus where EB denotes the energy of a bound state of the interacting 2-body system. 47 Let rt z (k) be the scattering phase shift of the potential v(r ) for the lth partial wave of wave number k. It can be shown that g(k) - g ( 0 ) (k) = _1-_ ""' ' (2l + 1) Brt z (k) 7r L....t z ok where the sum 2::::' extends over the values { 0, 2, 4, 6,... (bosons) l= 1,3,5,7, ... (fermions) Therefore we get the results: and after integration by parts, From Levinson's theorem, it can be proven that [23] L (2l + 1) rt z (O) = total numberof bound states = L (2l + 1) 7r l � (3.105) (3.108) 48 From above Eq. (3.10 8), one can see that the contribution from the continuum cancels the first term in the summation of the bound state terms. To repeat this procedure to one more order, it will lead to the so-called Planck-Larkin partition function. 3.1.5 The OPAL equation of state The OPAL equation of state (OPAL EOS) [56, 58, 54, 57] is part of the Opacity Project OPAL at Livermore. It was realized in so-called "physical" picture, and follows the theory exactly as we discussed in previous sections. It deals with the quantum plasma in activity expansion (ACTEX), and uses the classical regrouping procedure first, then deal with the short range quantum effect later. Quantum S function in activity expansion: (3.10 9) where a1 is the quantum-mechanical virial coefficient the dive rgence can be removed by the nodal procedure. After regrouping, we can have 00 (3.1 10) This looks like classical expre ssion, then we will need to find out the quantum mechan- ical counterpart for each terms, SR and S m . Recall the two component plasma gas with species-a and b, the cluster coefficients in classical case are with the form of: S2 ·· = -a2 ·· + z·z · 2 7r ( /3 e ·e · ) >- -- ( fJ e ·e · ) AD [ 2 7r 2 ] .�J .�J z J z J D 2 z J ' (3.111) 49 Start from here let's define and go to the quantum expression:[55] where a 2 .. >._3 o b - , t J - T ( - f3 H i j - (3 H . ') ij - - -- - - r e ') ZiZj 2 (3.1 12) (3.1 13) (3.1 14) (3.1 15) (3.1 16) (3.1 17) wij is the quantum correction term, and it will goes to 1 while taking the classical limit. Aij is the de Broglie wave length Aij = h / J2 7rP,ijkBT,J-lij is the reduced mass. The equation of state up to second order activity is [23] f3(p - Po ) = L zizj bij ij (3.1 18) where p0 stands for the pressure of the ideal gas. In the TCP, the scattering between particles might form new species. For example, in the plasma with electrons and protons, 50 the reactions between e � e or p � p pair would not produce new species since they repul se with each other, but the e � p pair has the chance to form the hydrogen atom. Hence, it will change the thermodynamical quantities, such as pressure or specific heats quite noticeably. From the previous discussion we can see the bep can be contributed by two parts, the bound state part b� P and the scattering continuum state part b� P [63, 4], (3.1 19) (3.12 0) (3.121) By doing the integration by parts twice, and apply the Levison theorem twice, one can demonstrate: [54, 38] A_ 3 b� p = ; p 2::) 2l + 1) ( e - f3 E nl � 1 + f3E nl) , nl (3.12 2) This effect involves the compensation between the scattering states and the bound states, and gives a convergent expre ssion of the internal partition function, which is often referred as the Planck-Larkin partition function: 51 PLPF = 2:: (2l + 1) (et3 En! - 1 + f3E nl) nl (3.12 4) Following this detailed discussion, one will not be mislead by the name of PLPF. The PLPF is not a real partition function, but rather a part of the second cluster coefficient, which reveals the results of the compensation between the scattering and the bound states. And b� P can be used to define the activity of new type of composite part icles. For exam ple, the bound state formed from electron-proton pair is the hydrogen atom, which is not assumed at the beginning. The activity of this new compound can be recognized as: = P _ 3 � [ - En!/kBT _ Enl] ZH _ 2zezpb ep- ZeZpA ep � (2l + 1) e 1 + k B T (3.125) In general , activity of the new bound states formed by the interaction between species a,b is defined as (3.12 6) The newly introduced activity z H has lowered the order to the first, and should be included in the original S 2 of activity { z e, zp}. S� repre sents the augmented set of activity (3.12 7) and s: P is the continuum state part of Sep (3.128) 52 This new activiti es is possible to form the higher order bound states, hence the heavier elements, such as helium. By considering all possible augmented set of activiti es m the second order, { z e, Zp,ZH.ZeH}, the equation of state is [49, 53] where i, j runs over the augmented set of activit ies: � nln 1 l 1 ZeH = D z eH nln 1 l 1 (3.12 9) (3.130) (3.131) (3.132) (3.133) (3.134) (3.135) 53 b'f j is the continuum state part of bij· For the pairs such as e - e or p - p they would not have the compensated parts. The last two terms other than b'f j in s'f j is from the Nodal expansion to remove the divergences. 3.2 Chemical picture The chemical picture is a more intuitive and heuristic way to deal with the equation of state. All kinds of atoms and their ions, as well as molecules are postulated to exist in the system. They all have their corresponding free energy terms to describe their interactions. There are reactions between species, given by ionization and molecular diss ociation (and their recombinations). In a method widely used in chemistry, conservation of the chemical potentials of all species defines the equilibrium, in which the numbers of each species can be determined. The procedure obtaining the equilibrium is known as free-energy- minimization method (FEMM). 3.2.1 Free-energy-minimization method Consider a physical system with a Hamiltonian H, confined in a box of volume V, in contact with a large heat reservoir at temperature T. The canonical partition function of the system can be expre ssed as, 54 H is the Hamiltonian of the system and f3 Helmholtz free energy is (3.137) k� T. From the partition function, the (3.138) Notice that the Ni here repre sents the number of particles within each species, like atoms, ions, and all the species are treated as "fundamental" particles on the same footing. This is ess entially the central idea of the chemical picture approach. At thermodynamic equilibrium, the free energy energy is minimized (corresponding to maximum entropy, which is the embodiment of the thermodynamic equilibriu m): (3.139) This condition will lead to the ionization equilibrium among different ionized and neutral species. This is known as the free energy minimization method (FEMM) [42, 30, 31, 32]. The equation of state can be simply derived from the free by differentiation (3.14 0) Another important assum ption in chemical picture is the modularity of the free energy [30]. In other words, the many-body interaction potentia1U(r 1, r2 , ... rN) in the Hamiltonian: 55 (3.141) does not depend on momenta, and therefore the partition function can be separated and written as the multiplication of different parts: _ 1 N Pi 1 - f3U (3 E ; !! N 2 !! N Z- h 3 N . . . dp exp(,B L 2m ) x N ! . . . dre x L g i e t=l t=l (3.14 2) Where Zt r an s is the translational free energy coming from the particle momenta: 1 N � J J N 2 Zt r ans = h 3 N . . . dp exp - (3 L 2m t=l (3.143) Confi gurational terms and internal terms are both from the contributi on of the interac tion potential. Configurational terms describes the interaction between different elements, or particles. For example, the Lenard-Jones (LJ) potential between neutral particles or the Coulomb between charg ed particles. Z confi g = � ! J ... J dr N exp[-f3U(r 1, r2, ... r"N)] (3.14 4) The internal terms counts for the internal energy levels of bound species, therefore they include the interaction between the nuclei and the surrounding electrons. N Zint = L 9 i exp(f3E i) i=l (3.145) 56 Since electrons are much lighter than all other species (which have at least the mass of the proton) they should be treated as a non-relativistic quantum ideal gas [13]. In the spirit of separation of effects, the electron gas term can be also singled out. Meanwhile, the free energy will be: F = - KBT ln(ZeZt r an s Zconfig Zint ) = Fe + Ft r an s + Fconfig + Fint (3.14 6) The modular form of the free energy makes analysis and modification easy. First, changing one term will not affect the others, therefore we can study the effect from differ ent physical contribu tions. Second, it can be seen that the largest contribution will be the translational and the degenerate gas term, in other words, the ideal gas. Then comes the correction of the inter-particle interaction, that is, the configuration term. By this reason ing we can see a hope to improve the equation of state by adding more minor, higher order contributi on from various physical effec ts. 3.2.2 The MHD equation of state The internal partition function of the hydrogen atom, which has the Coulomb potential and its discrete energy level s, is 57 Zint = l)2l + 1) exp(- E nz /kBT) nl 00 = Ln 2 exp(- E nz /kBT) . n=O (3.147) (3.148) Inserting the hydrogen energy level s of an isolated hydrogen atom - l3 � e v in Zint . we n can see that, (3.14 9) if the summation goes to infinity, then the first and second terms are divergent. Phys ically, this will never be a problem because in real world there is no such thing as iso late hydrogen atom. The atom interact with the environment and the infinite levels got destroyed and leaves the finite energy level s. The key question is how to treat the inftu- ence of the plasma environment on the internal structure of atoms. If one models these pro cesses, one obtains the appropriate internal partition function, from which the free energy ensues. There are many cut-off techniques proposed to solve this divergence problem and we just list some of them: • Saha equation: only includes the ground state energy as the internal energy. • Confined atom( CA): assumed that bound electrons are confined to a sphere deter- mined by the average volume per ion. • Statically screened Coulomb potential (SSCP): the screened Coulomb interaction limits the number of state s. 58 They are all simplified models and referred as 'hard-cutoff' which have their own flaws, like no pressure ionization in Saha equation and the disc ontinuity in CA model. And the iss ue of the SSCP has been discussed in Mihalas and Hummer 's paper. [35] The MHD model was developed by D. Mihalas, D. G. Hummer, and W. Dappen [35, 16, 44]. Unlike 'hard cutoffs ', MHD models the distribution of bound states by occupation probability, and is usually referred as a 'so ft cutoff '. Specific physical processes were considered to obtain "occupational probabili ties" Wis as continuous and monotonically decreasing functions representing the probability of finding species-s at its i-th ionization state. By introducing Wis the energy levels can disappear smoothly. Z. _ "' . . e- f3 E ;s znt - L W zs9zs · (3.150) i ,s 3.2.3 EFF & CEFF EFF model is a simple empirical model developed by Eggleton et al. [24]. It is a model that is not based on physical foundations, but rather a numerical fit forced to produce the anticipated results, ie full pressure ionization at high densities. However, despite its lack of physical foundation, it was immediately better suited for stellar modeling pur poses than earlier models that had imposed full ionization beyond an empirical determined temperature -density curve. Also, the EFF model has the advantage of being easy to realize by a computer program that can be used in-line inside large stellar structure and evolution codes. It has been fairly successful in many solar models[8]. EFF was later upgrade to incorporate the correction from the screened Coulomb potential. This upgrade is com monly referred to as CEFF. Based on its usefulne ss, our project will be based on upgrades 59 of CEFF to achieve a succes sful OPAL emulator. A more thorough discussion of the MHD and EFF/CEFF EOS is in the subsequent Chapter 4. 60 Chapter 4 MUD and CEFF EOS and their modifications In this chapter, we review the theory of the MHD EOS and CEFF EOS. Then we present the approach of modifying MHD EOS to an OPAL emulator which is done by our former members in this group. We then show our approach of modifying CEFF EOS to an OPAL emulator. 4.1 MUD As discussed in Section 3.2, we know that within the chemical picture, all thermodynam- ical quantities are derived from the free energy. The main issue is how to cut off the infinite energy level of the isolated hydrogen atom in order to yield more realistic proper ties of the hydrogen plasma. In MHD formalism, a 'so ft' cutoff procedure was introduced. Specifically, by introducing a continuous and monotonically decreasing function Wis. the occupational probability, the MHD internal partition functions are written as (4.1) i ,s with Wis defined by 61 (4.2) Wis represents the occupational probability of i-th energy level of species s. v runs over neutral particles while a runs over charged particles except electrons. rij is the radius asso ciated with state i of species j. Zj is the net charge of species j, and X is is the ionization energy of level i, and Kis is from quantum correcti on. (4.3) The free energy terms are included in MHD ([35, 44, 16]) are: (4.4) • F1 is the translational free energy term for a ideal gas without electron s. 3 F1 = -kBT L Ns( 2 lnT+ln V- lnNs+lnGs+l),Gs = (2Trm s kB/h 2 ) 3 1 2 (4.5) s fe • F2 is the internal free energy term for components that have internal degrees of freedom, such as atoms, ions and molecul es, F 2 = kBT L Ns ln Zs, Zs = L Wis g ise - f3 E is s fe (4.6) 62 • F 3 is the free energy contributi on from the partially degenerate electron, which was treated as partially degenerate quantum gas. (4.7) F n (rt) is the standard Fermi-Dirac integral of order n, rt is the degeneracy parameter, related to the number of electrons by F1 ; 2 ( rt) = tJ. v��/ 2 . • F 4 is the confi gurational free energy term, which accounts for the lowest order environmental Coulomb interactions among electrons and protons in the continuum(Sri n g term), (4.8) The equation of state can be deduced by (4.9) And from here all the thermodynamical quantities can be derived. In particular, we are interested in the adiabatic exponent 11, the strain coefficient Xt , and the isothermal compression coefficient X p · We will calculate those quantities along with pressure, and analyze the effects on them from different physical models. (8lnp) p (8p) /1 = 8 ln p ad = p 8 p ad (4. 10) 63 = (� ) = r (ap) Xt 81n T aT p p p (8lnp) p (ap) X p = aln p T = p ap T (4.11) (4. 12) Before the confirmation from the observations, we can not really say either physical or chemical picture is better. Since physical picture are rigorous however started from the elementary particles, hence to get the effect from the heavy elements will inevitably require approximation method or much more intensive calculat ions. On the other hand, the chemical picture will have no problem bringing any species if demanded however it is more heuristic and also need some basic assum ptions. Overall, the thermal quantities obtained from the OPAL EOS is closer match to the observation data than from the MHD EOS, which requires us to introduce more refined physical effect in MHD EOS to improve its accuracy. 4.2 Modified MHD 4.2.1 Planck-Larkin partition function as internal partition function The implementation of PLPF within chemical-picture formalism, (e.g. MHD, CEFF) has been rejected and doubted, and it was a subject that has been intensively debate[14]. Although PLPF is called "partition function", it is not a real partition function but rather an part of the contribution of the partition function from the physical picture(i.e. the sec ond virial coeffi cient). However, PLPF does reveal the cancelation between the high end of the bound states and low end of the scattering continuum of a bound species, hence 64 emerge a convergent, temperature dependent series to count the contribution to internal partition function, which can solve the infinite problem of hydrogen atom. To adopt the PLPF into chemical picture still lacks its physical legitimacy, but the results turn out to be well-pe rformed according to the previous and present experience. To be consistent, the higher order terms of the continuum- state-2-body Coulomb interaction of electrons and protons, are also included. in other words, a new occupation probability is introduced as [20, 21, 22] So the internal partition function will look like: and z s F 2 = - k BT L_, Ns lnZs . s ol e (4. 13) (4. 14) (4. 15) (4. 16) Remember the resulting convergent internal partition function is because of the com pensation from the scattering states terms. Those scattering states are also not originally in MHD therefore need also be included. 65 4.2.2 Scattering terms in Free Energy The first order term will involve the 2-body interaction. Since the most part of the Sun is constitute of hydrogen, we will anticipate that for most contributi ons we only need to consider hydrogen plasma. Even the second most element, helium, is about ten times smaller in number compare to hydrogen. Therefore let us discuss the terms for hydrogen plasma first. In hydrogen plasma, the basic components would be electrons e, protons p, hydrogen atom H, and the hydrogen molecules H2 • The Coulomb interaction among electrons and protons dominates over the weak interaction, like Van der Waals force between either charged or neutral species. Therefore , only 2-body interaction between charged particles is included in the scattering free energy terms. (4. 17) Recall the relation between free energy and the S function: (4. 18) with F and F0 are the free energies for the interacting real gas system and the ideal gas system, respectively. The free energy of the scattering terms therefore can be: 66 (4. 19) It leads to (4.20) (4.21) (4.22) and (4.23) where (4.24) Since the mass of proton is about 1840 times of electron's, so although we can treat protons classically but electron should be considered quantum-mechanically. bw was treated as in classical case, (4.25) 67 sPP = 27r 100 [e- f3q (r) -1]r 2 dr + 27r Ab(/3e 2 ) -�AD(B e 2 ) 2 AD is the Debye-Hi.ickel screening length, AD= JsokBTV je 2 (N e + Z'j; N P ) . (4.26) And for the partially degenerate electron gas, a semi-classical Wigner-Kirkwood expansion approach was used.[48, 36] To obtain the explicit expre ssion of bep• the perturbati on expansion in weakly coupling plasma is used. The quantum diffracti on /ep = Aep/ AD is less than 0.4 in solarplasma.[39]. The expre ssion of bep can only be worked out up to the second perturbati on, the higher terms can only be obtained at the limit of /ep ---+ 0. [17, 33] (4.31) 68 b ; P = - i({3 e 2 ) 3 (ln �; + D q ) Aep is thermal de Broglie wavelength divided by v'47[: Aep = nj )2 /-lepkBT. D q = ln3 + � ( IE - 1) = 0.8872, /E is Euler's constant.[50, 52] 4.2.3 More scattering terms (4.32) (4.33) (4.34) (4.35) If we not only consider the pure hydrogen plasma, but the H-He mixture, then obviously we will need more scattering terms related to the Helium and the the PLPF of He and He+. We know that comparing with hydrogen, helium is quite a small fraction, and it is exactly the reason that this effect hasn't been included in the previous work. However, since extending the procedure is straight forward, we would like to add more terms here and study the effects they bring. We list the terms that will be added: (4.36) 69 FHe++He++ = FpHe++ = k B T { rX> f3 (r) 2 2 2 'iT 2 2 } - 2NHe++Np V 2 'iT lo [ e q - l]r dr + 27r Av(f3Z iZj e ) - 2 Av (f3Zi Zj e ) (4.38) (4.40) FHe+He+ = 70 FHe+He++ = k B T { (X:> [ !3 (r) ] 2 2 2 7r 2 2 } - 2NHe++N He+ V 21r J o e q - 1 r dr + 2TrA. D(f3Z iZj e )- 2 A.D(f3Z iZj e ) (4.42) (4.43) 4.3 EFF and CEFF In the following section, we will discuss the main idea from the paper in 1973 by P. Eggle ton, J. Faulkner and B. Flannery [24], which leads the EFF equation of state and later the CEFF equation of state. CEFF has become the most prevailing practical equation of state among solar evolutionary model and helioseismic analysis, its compact form and fast calculation speed yet still with good performan ce are the main reasons. The following describes the components that within CEFF equation of state. 4.3.1 Approximation to Fermi-Dirac Integrals The density of a gas is given by where T is the temperature in unit of mc 2 / k and 1/J is the electron degeneracy param eter. Ae is the Compton wavelength of the electron and the energy-momentum with the 71 relation 2E + E 2 = p 2 , which coped with special relativity. !le is the atomic weight per free electron. Chandrasekhar[6] gives expansions for p* in four limits, some convergent, some asymptoti c. p* rv T 3 1 2 e .P 2:::2::: a� � e m .Prn T 3 1 2 e .P 2:::2::: a �� e m .Pr-n (T 1/J ) 3 / 2 2::: 2::: a� � 1/J - 2m (T 1/J ) n (T 1/J ) 3 2::: 2::: a �� 1/J - 2m (T 1/J t ( 1/J « -1,T « 1) ( 1/J « -1,T » 1) (1/J » 1,T 1/J « 1) (1/J » 1,T 1/J » 1) All the summations are from zero to infinity unless indicated. The authors suggested if there are two functions f(1/J) and g(1/;, T) which have the properties such that and f ( 1/J) rv e .P l:b � ) e m .P 1/; 2 2::: b � ) 1/J - 2m g ( 1/J, T) rv T 2::: b �3 ) e n,P rv 1/JT 2::: b � 4 ) 1/J - 2 n ( 1/J « - 1) ( 1/; » 1) ( 1/J « - 1) ( 1/; » 1) (4.44) (4.45) (4.46) (4.47) 72 Then the original series expansion of p* can be rewritten in terms of g and f: p* rv l / 2 J L L C��J m g n 9 3 J L L C��J m g - n l / 2 L L C ��J - m g n 9 3 L L C� �J - m g n Now note that for any M,N � 1 the expression f "'M "' N A f m n PA = _ 9 2 /3(1 + g) 2 /3 L..o L..o Pmn g - l + f (l+J)M(l+g) (! « 1,g « 1) (! « 1,g » 1) ( ! » 1,g « 1) ( ! » 1,g » 1) has the same limiting behavior in all four limits of p*. Hence the authors can choose the the form of function f ( 1/J), g ( 1/J, T) , the value of M, N and the coefficients P mn in order to fit the exact value over the whole ( ! , g) plane. Equation ( 4.44), ( 4.45) imply that 1/J rv l og( ! ) + L d� ) J m rv v12:: d �) j- m One possible choice of 1/J - f relation can be d'¢ vT+1 df f (! « 1), (!» 1) . (4.48) 73 which the solution is (4.49) Simila rly, to have the behavior of ( 4.46) and ( 4.47), one can have the relation between g and f: g=T vT+f . (4.50) One should see these choices are not unique and they merely serve the purpose of simplifying the calculat ions and still yield good approximation of Fermi-integral. Through the numerical tests by the authors, the accuracy and reliability of these choices has been reported sufficient for many astrophysical applications . With these approximations, pressure and internal energy of electron gas can be written as follow s: 74 4.3.2 An Approximate Treatment of Pressure Ionization The ionization equilibrium was determined by the Boltzmann-Saha equation in EFF. This means only the ground states was counted as the internal structure of the bound species such as atoms or ions. The equilibrium between two ad jacent states of ionizations can be written as: (4.51) where n +, n, w+, w are the abu ndances and statistical weights of two adjacent states of ionization, x being the ionization potential, T is the temperature in Kelvin. This implies that the electron sitting other than ground state are essen tially treated as free electron in this approximation. However, while the the density becomes high and the ions will interact with the neighbor ions, and the surrounding free electrons, hence the w+, w, x should be altered by the temperature and density dependent terms. The authors modifies the original Boltzmann-Saha equation into the form: (4.52) The perturbati ons to the Boltzmann-Saha equation will be included through l1f-l term. In principle, l1f-l should be function of the density of the ions of different species and the free electron s, i.e. function of n +, n, ne, where ne has the form: (4.53) 75 In the original EFF paper, ,6.1l was simplified to be depending only on the free electron density ne, with that simplification the equation ( 4.52) can have the explicit form instead of being solved recursively. And this is still fairly reasonable if the material is largely ionized, hence the effect from the free electrons will dominate. The correction of Helmholtz free energy ,6.F was postulated as ( ) O (T) ( 2 2 ) ,6.F Ne, V, T = ----:y:- N eo - N e (4.54) with specific volume V, the number Ne of free electrons, O (T) some suitable function of temperature only, and Neo is a constant, the total number of electrons (bound or free) per unit mass. Consequently, meanwhile, ,6. _ _ _ 1 8,6.F _ 20 n ll- kT ( 8Ne )v,T - kT e where neo is the total number density of elections, i.e. Neo/V. The choice of 0 is 20 = a6( kT + 20xo) where (4.55) (4.56) (4.57) 76 a o = 0.523/ Z(A) Xo = 13.5Z2(eV) Z is the mean charge per nucle us. This choice of f2 is loosely equivalent to saying that the statistical weight w is decreased relative to w+ by a factor of e-a6f r 3-a6f r 'b, where r is mean inter-electron distance and rv is the Debye radius. The choice numerical number in the definition of f2 is to ensure that when the mean inter-electron distance is of the order of a0 the has become fully pressure ionized. 4.3.3 Some issues The method described in the previous sections shows that the state of the gas in a star can be calculated explicitly, and very quickly by a single formulati on applying over a wide range of temperature and electron degeneracy. However, there are a few exceptions: • At relatively low temperature and high densities the gas can be strong imperfect that the simple approximation to pressure ionization will not be useful or it might be furthermore dominated by other effects like crystallization. • At low density and very high temperature pair creation becomes important. • At very high density and/or temperature there may be nuclear eff ects, and more complicated formulation seems inevitable. Aside with the above issues, there is another important effect that al so need to be considered, which is the contribution to the free energy due to the electrostatic interactions between charged particles, i.e., ions, nuclei , and electrons. To add this effect accounts the non-ideality of the gas and render EFF to so-called CEFF. 77 4.3.4 CEFF According to the De bye- Hi.ickel model, the free energy corresponding to the electrostatic interactions will has the form, where the sum on o: runs over all charg ed species, and ea 1( o: �e), = ( 21r 1 l 2 e 3 ) 1 F 1 ;2 (77) l:a f e Na Za ( � N z 2 e ) 1/2 X - k3 /2 Vl/2 T3 /2 F 3 /2 (TJ) l:a f e Na � a a a ' (4.58) (4.59) (4.60) (4.61) (4.62) where F n(TJ) is the standard Fermi-Dirac integral of order n, and 77 is the degeneracy parameter, which is related to the number of electrons by (4.63) This F4 term is actually transferred from the MHD package, therefore there will be subtle details of the matching of physical unit, constant, and transformation of variables. 78 Those details will be discussed more in Section 4.4. Here, the EFF equation of state was modified to include the Coulomb effect in a similar way to include the pressure ionization. + + '!!:_ = '!!!_ e - .P - x / kT + l:ifk + l:ifkc n w where the F4 comes into EFF through the relation � _ _ _l_( oF4 ) J-lc - kT o Ne V,T (4.64) By adding this Coulomb correction the resulting equation of state is so-called CEFF equation of state. CEFF EOS hence become popular owing to its effi ciency and accuracy. 4.4 Details of modified CEFF The physics we added in CEFF is essen tially the same as those has been used on modified MHD equation of state, hence we won't repeat here. However, since CEFF itself is differ ent equation of state than MHD,in many aspects, so carrying the work from one to another will need some care. In this section we will present the technical details of modifying the CEFF program. We continue the previous work done by USC group, which had modified the MHD for malism towards an early attempt to emulate the OPAL equation of state. Our strategy here is to carry those modifications over into the CEFF equation of state. 79 4.4.1 Variable transf ormation MHD and CEFF are both realized in so called "chemical picture", however they do use different set of variables to obtain the ionization equilibrium. In MHD, the number den sity of the species like hydrogen atom, hydrogen plus ion, helium atom, helium plus ion, helium plus plus ion, and free electrons(i.e. NH,NH+,NHe ,N He+,N He++,Ne) are used. In CEFF, instead of number density of each species, a set of variabl es called "reaction parameter" has been used. Reaction parameter is a measure of certain chemical reaction have taken place. In CEFF the following reactions are considered: dl : H ---+ H + + e- d2 : He ---+ H: + e- d3 : H: ---+ H: + + e- The relation between the number density and the reaction parameter is easily observed: NH = CH -dl NH+ = dl NHe = CHe- d2 NHe+ = d2 - d3 NHe++ = d3 Ne = d 1 + d2 + d3 Where CH and CH e are the number density of hydrogen nuclei and helium nucl ei, respectively. Hence we know the free energy describe in MHD is (4.65) and the free energy in CEFF will be 80 Fc E FF = Fc E FF(d1,d2, d3), and those two expression should have identical physical meaning. 4.4.2 Transf ormations in first order derivatives (4.66) The first derivatives of the free energy with respect to their own argument will have the relation: 8F 8F adi = L aN/ where i=1,2,3, j=N H,NH+, ... Let's write them down explicitly. 8F 8F 8F 8F - = - � - + -- + - 8dt 8NH 8NH+ 8Ne (4.67) (4.68) (4.69) (4.70) Now let's introduce another set of variables "ionization degree", which are defined as(number of ion/number of nuclei). Define (4.71) 81 �2 = NHe+ = d2 - d3 CHe CHe as ionization degree of NH+ ,N He+,N He++· They obviously can be inverted to the relation: d1 = CH O d2 = CHe(�2 + �3) d3 = CHe �3 and we have the relation with number density for each species: NH = CH (1- �1) NH+ =CH O NHe = CHe (1- �2 -�3 ) Therefore from we will have (4.72) (4.73) 82 To write them explicitly, 8F 8F 8dj 8�i = 2:: 8dj 8�i · J 8F _ 8F 8dl 8F 8d2 8F 8d3 _ 8F 0 0 0 _ 8F 0 80 - 8dl 80 + 8d2 80 + 8d3 80 - 8dl H + + - 8dl H 8F _ 8F 8dl 8F 8d2 8F 8d3 _ 0 8F 0 0 _ 8F 0 8�2 - 8d1 8�2 + 8d2 8�2 + 8d3 8�2 - + 8d2 He + - 8d2 He (4.75) 8F 8F 8dl 8F 8d2 8F 8d3 8F 8F 8F 8F 8�3 = 8d1 80 + 8d2 8�3 + 8d3 8�3 = O + 8d2 CHe + 8d3 CHe = CHe( 8d2 + 8d3 ) (4.76) There is another set of variabl e, the 'number fractions', defined as (4.77) 83 N He + (d2 - d3) W2 = NHe = ---:c =-- !'!..iil._ N'i£ e C He C He (d2 - d3)(NHe + NHe+ + NHe++) N'fi e 1 1 = N 2 [- d3NHe + d2(NHe + NHe+ )] + N 2 [- d3(NHe+ + NHe++) + d2NHe++] He He 1 1 = - 2 - [- d3NHe + d2(NHe + NHe+)] + - 2 [(d2 - d3)NHe++- d3NHe+] N He N He 1 1 = - 2 -[ -d3NHe + d2(NHe + NHe+)] + - 2 -[ NHe+N He++ - NHe++NHe+] N He N He 1 = - 2 - [-d3NHe + d2(NHe + NHe+ )] N He The inverse relation is thus: 84 Simila rly, from Fc E PP(dl, d2, d3) = Fc E PP(dl (Wl, W2, W3), d2(Wl, W2, W3), d3(Wl , W2, W3)) , we will have To write them explicitly, aP aP Bdj 8Wi = L 8dj 8Wi" J aP aP Bdl aP Bd2 aP Bd3 aP Nk aP Nk 8Wl = 8d1 8Wl + 8d2 8Wl + 8d38Wl = 8dl CH + 0 + 0 = 8dl CH (4.78) - - - - - 2 - 2 8F 8F 8dl 8F 8d2 8F 8d3 8F N He ( NHe++) 8F N He NHe++ 8W2 = 8d1 8W2 + 8d2 8W2 + 8d3 8W2 = O + 8d2 CHe l + NHe+ + 8d3 CHe NHe+ - - - - - 2 - 2 8F = 8F 8dl + 8F 8d2 + 8F 8d3 = O + 8F N He NHe+ + 8F N He (NHe+ + ( NHe+ ) 2 ) 8W3 8d1 8W3 8d2 8W3 8d3 8W3 8d2 CHe NHe 8d3 CHe NHe NHe Now consider the derivative respect to mass fraction of hydrogen and helium. Here we do make an approximation, which is to ignore the mass fraction of the heave elements since they have small contributi on while obtaining the ionization equilibrium. 85 Consider X+Y = l X is the mass fraction of hydrogen and Y is the mass fraction of helium. Where Mtot =AHCH + AHeCHe, with AH, AHe are the atomic weight of hydrogen and helium atoms. respectively. This relation can be inverse to Note now the derivatives of the free energy were made to the set of parameter [X, Mtot l. or [CH,CHel Hence the reaction parameter [d i] or the ionization degree [(i] or the number fraction [Wi] were held as constants. 86 4.4.3 Transf ormations in second order derivatives The mixed derivatives with respect to reaction parameter and ionization degree: aP &P &P &P 8d180 = CH( 8N H8N H 8NH8NH+ 8N H8N e &P &P &P - + + --- 8N H+8NH 8N H+8N H+ 8N H+8N e 8 2 P 8 2 P 8 2 P - + + ) 8Ne8NH 8Ne8NH+ 8Ne8N e 87 8 p 8 2 p 8 2 p 8 2 p fJdlfJ�2 = CHe( fJNHfJNHe fJNHfJNHe+ fJNHfJN e fJ 2 ft' 8 2 ft' fJ 2 ft' - + + --- - 8N H+0NHe fJN H+fJN He+ fJN H+ONe 8 2 ft' 8 2 ft' 8 2 ft' - + + ) 8Ne8NHe 8Ne8NHe+ 8N e8N e a ft' a 2 F fJdlfJ�3 = CHe( fJNH fJNHe fJ 2 ft' ----- -2 --- fJNH fJNHe+ + 8N H8N e &P &P &P - + + 2 --- - fJN H+ONHe fJN H+fJN He++ 8N H+8N e 8 2 ft' 8 2 ft' 8 2 ft' - + + 2 ) 8Ne8NHe 8Ne8NHe++ 8Ne8Ne 8 p 8 2 p 8 2 p 8 2 p --- = CH( ---- 8d280 8N He8N H 8NHe8NH+ 8N He8N e 8 2 ft' 8 2 ft' 8 2 ft' ---- + + ---- fJNHe+ON H 8NHe+8N H+ fJN He+ONe fJ 2 ft' fJ 2 ft' fJ 2 ft' --- + + ) 8N e8N H 8Ne8NH+ 8Ne8N e 88 8 p 8 2 p 8 2 p 8 2 p fJd2fJ�2 = CHe( fJNHeONHe fJNHeONHe+ fJN HeONe fJ 2 ft' fJ 2 ft' fJ 2 ft' - + + ---- fJN He+ONHe 8NHe+8N He+ fJNHe+ON e 8 2 ft' fJ 2 ft' 8 2 ft' - + + ) 8Ne8NHe 8Ne8NHe+ 8Ne8N e 89 8ft - c ( 8 2 ft 8 2 ft - 2 8 2 ft 8d38�3 - He 8NHe+8N He 8NHe+8N He++ 8N He+8N e &ft &ft &ft - + +2 ---- 8NHe++8N He 8N He++8N He++ 8NHe ++8N e 8 2 ft 8 2 ft 8 2 ft - + + 2 ) 8Ne8NHe 8Ne8NHe++ 8Ne8N e The mixed derivatives with respect to reaction parameter and number fractions: 90 - 2 oF N He ( NHe++ ) --::-- ----::- - = - - 1 + * 8dl8W2 CHe NHe+ (]2 ft (]2 ft (]2 ft [ 8NH8NHe 8NH8NHe+ 8NH8N e 8 2 ft 8 2 ft 8 2 ft - + + ---- 8N H+0NHe 8N H+0NHe+ 8N H+0Ne (]2 ft (]2 ft (]2 ft - + + l 8Ne8NHe 8Ne8NHe+ 8Ne8N e ( NAe NHe++ ) + - - * CHe NHe + (]2 ft (]2 ft (]2 ft [ 8NH8NHe+ 8NH8N He++ 8NH8N e (]2 ft (]2 ft (]2 ft - + + ---- 8N H+0NHe+ a NH+ a NHe++ 8N H+0Ne 8 2 ft 8 2 ft 8 2 ft - + + l 8Ne8NHe+ 8Ne8NHe++ 8N e8N e 91 - 2 oF =( N He NHe + ) * 8dl8W3 CHe NHe (]2 ft (]2 ft (]2 ft [ 8NH8NHe 8NH8NHe+ 8NH8N e 8 2 ft' 8 2 ft' 8 2 ft' - + + --::-- - --=-- - 8N H+0NHe 8NH+8N He+ aN H+ONe (]2 ft 8 2 ft' 8 2 ft' - + + l 8Ne8NHe 8Ne8NHe+ 8Ne8N e N'fi e ( NHe+ ( NHe+ ) 2 ) + -- --- + -- * CHe NHe NHe 8 2 ft' 8 2 ft' 8 2 ft' [ 8NH8NHe+ 8NH8N He++ 8NH8N e (]2 ft (]2 ft 8 2 ft' - + + ---- 8N H+0NHe+ 8N H+0NHe++ 8N H+0Ne 8 2 ft' (]2 ft' 8 2 ft' - + + l 8Ne8NHe+ 8Ne8NHe++ 8N e8N e 92 93 94 oF N'ke ( NHe++ ) �----=- - =-- 1 + * 8d38W2 CHe NHe+ a 2 F a 2 F a 2 F [ ---- 8NHe+0N He a NHe+ONHe+ a NHe+ONe 8 2 ft' 8 2 ft' 8 2 ft' - + + --------- 8N He++0N He a NHe++ON He+ a NHe++ON e a 2 P a 2 P a 2 P - + + l 8Ne8NHe 8Ne8NHe+ 8Ne8Ne + ( Nfie NHe++ ) * CHe NHe + a 2 F a 2 F a 2 F [ 8NHe+8N He+ 8NHe+8N He++ 8NHe+8N e 8 2 ft' 8 2 ft' 8 2 ft' - + + --------- 8N He++0N He+ a NHe++ONHe++ a NHe++ON e 8 2 ft' 8 2 ft' 8 2 ft' - + + l 8Ne8NHe+ 8Ne8N He++ 8Ne8N e 95 oF =( N 'ke NHe + ) * 8d38W3 CHe NHe 8 2 ft' (]2 ft' 8 2 ft' [ ---- 8NHe+0N He 8NHe +8N He+ 8NHe+8N e 8 2 ft' a 2 P a 2 P - + +� ---- �- 8N He++0N He 8NHe ++8N He+ a NHe++ONe 0 2 p 8 2-p 0 2 p - + + l 8Ne8NHe 8Ne8NHe+ 8Ne8Ne N'fi e ( NHe+ ( NHe+ ) 2 ) + - - ---- + ---- * CHe NHe NHe 8 2 p 0 2 -p 8 2 p [ a NHe+ a NHe+ a NHe+ a NHe++ 8NHe+8N e (]2 ft' (]2 ft' (]2 ft' - + + --------- 8N He++0N He+ a NHe++ONHe++ a NHe++ON e 8 2 ft' (]2 ft' 8 2 ft' - + + l 8Ne8NHe+ 8Ne8N He++ 8Ne8N e 96 The mixed derivatives with respect to reaction parameter and X(mass fraction of hydro- gen): aP Mtot 8d18X = CHAH * [( 97 98 4.4.4 Switching off pressure ionization The switch of pressure ionization was controlled by the parameters ca03, and caa within CEFF EOS program. Since in our modified CEFF, the PLPF will take care of pressure ionization, the original pressure ionization device of EFF is no longer needed; we are switching it off therefore by setting the relevant pair of parameters to 10 - 1 0 . In this way, EFF pressure ionization is disabled and PLPF will dominate the ionization at the high temperature end. 99 Chapter 5 Results and Discussion In this chapter we show the results from both MHD-based OPAL emulator and CEFF based OPAL emulator. The results from the MHD-based OPAL emulator, which is done by the former member of this group is also shown, in order to compare with the results from the CEFF-based OPAL emulator, which is done in this work. We then do further modificati ons on the MHD-based OPAL emulator with H-He plasma. We test how those He-related scattering terms affect the quality of the equation of state and then also carry those modifications to the CEFF -based OPAL emulator. 5.1 Hydrogen Plasma 5.1.1 PLPF and the Scattering terms From their previous research of an OPAL simulator, which was based on MHD[39], A. Liang and W. Dappen have found that the key role is played by the Planck-Larkin partition function(PLPF) : """' [ ( E nl ) E nl ] P LP F = �(2l + 1) exp - kT - 1 + k B T (5.1) 100 However, to be consistent, the second order scattering term should also be considered: Here, Aep is the thermal de Broglie wavelength divided by y'47f, Aep = h/ J2 fl ep kBT' >.. v the Debye-Hi.ickel radius JEokBTB je 2 (N e + z;N P ), and the constant Dq=l.72. The number Dq here is 1.72 instead of the original theory value 0.82. The reason for using this number is that it renders the best match between the MHD-based OPAL emulator and the OPAL EOS. To show our present results, we plot the thermodynamic quantities from the different model of equation of state, and we use the OPAL equation of states as the reference. 101 We have plotted the following physical quantities: olnp p op /1 = ( � l ) ad = � (�)ad u np p up (5.2) (5.3) (5.4) (5.5) and the relative difference w.r.t OPAL among the models we are playing with (i.e. CEFF, or MHD) was calculated as d"( Quantity) Quantity Model - QuantityOPAL Quantity Quantity OPAL 5.1.2 MHD and modified MHD (5.6) As presented in the previous work[39], it was found that the PLPF plays an important role. In that previous work, the terms involving the scattering terms were multiplied with ad justable parameter so as to achieve optimal OPAL emulation. Specifically, in [15], there is an parameter D q in Fep that was optimized in order to meet the purpose of matching the OP AL. After doing this optimizati on, the contribution from the PLPF and the scattering terms become equally dominated. The results we show here is following the results from the MHD-based OPAL emulator with the optimized parameter. 102 Fig. 5.1 shows the relative difference to OPAL model in original MHD model and the modified MHD model, with PLPF in hydrogen atom, and the extra scattering terms Fee, F PP ' and Fep· In this graph, the OPAL itself represents the horizontal zero line since the difference to itself is of course zero. Hence our goal here is to bring our models as close to zero as possible. We see the results from the modified MHD has been much improved, compared to the original MHD. XT r 0.01 0.006 0.008 MHD --+- 0.005 MHD --+- 0.006 D(sc+ PL) --- x -- MHD(sc+ PL) - - x--- 0.004 0.004 I- 0.002 0.003 t 0 f::c 0.002 !:.-. X -0.002 <.() <.() 0.001 -0.004 0 -0.006 -0.008 -0.001 -0.01 -0.002 4 4.5 5 5.5 6 6.5 7 4 4.5 5 5.5 6 6.5 7 logT logT pressure X r 0.0001 0.004 5e-05 MHD --+- - -� l�J JJ � � 0.002 0 -5e-05 0 Q: -0.0001 a. � -0.002 (l_ -0.00015 a. <.() X <.() -0.0002 -0.004 -0.00025 -0.006 -0.0003 -0.00035 -0.008 4 4.5 5 5.5 6 6.5 7 4 4.5 5 5.5 6 6.5 7 logT logT 103 Figure 5.1: Compari son of relative difference (to OPAL model) quanti ties between the original MHD (MHD)and the modified MHD with PLPF and extra scattering term s(PL+ex). The result of the modified MHD has been brought closer to OPAL sig nificantly. 5.1.3 CEFF and modified CEFF In the present project, our goal is to make the modified CEFF EOS, in a similar way as been done with MHD. The main physical idea is the same, and what has to be overcome are the practical details coming from the different platform between MHD and CEFF packa ges, as described ab ove. 104 0.015 0.01 0.005 I- t 0 X vO -0.005 -0.01 -0.015 -0.02 4 0.0003 0.0002 0.0001 a.. a: 0 vO -0.0001 -0.0002 -0.0003 4 Xr , 11 , PL +ex --+- I I ' ' )�� *x� X " * - * \ I I \ I I I � :j 'x 4.5 5 5.5 logT pressure ex --- \( -- PL -- - -)IE--- 6 6.5 X ,, PL +ex --+- /" \ ex - -- \( --- f: I PL ----)1(--- t I I I ;.W* \ ,i+: � -)( � 0 4.5 5 5.5 6 6.5 logT 0.003 0.002 0.001 f::c """' 0 vO -0.001 -0.002 -0.003 7 4 0.008 0.006 0.004 0.002 a. .:;::. a. 0 X vO -0.002 -0.004 -0.006 -0.008 7 4 � I r iii PL +ex --+- : If{' ex - - x --- : ' PL - -- -)IE--- ' :;: ' I ' , f\ "' '*x W<� t\ \ ' \ : I I I I It II !{ 4.5 l'x 1 I I " I I I I .j 5 5.5 logT X r X * ,;�'!· 4.5 5 5.5 logT 6 6.5 PL +ex --+- 6 6.5 7 7 Figure 5.2: Compari son of relative difference (to OPAL model) quantiti es between the CEFF with PLPF and scattering term s(PL+ex), CEFF with only PLPF(P L), and CEFF with only extra scattering terms(ex). 105 [b] 0. 01 5 .--------..---r---r-----,r-----r---, """" O.D1 0.005 0 � £ -0.005 -0.01 4.5 5 5.5 6 6.5 7 logT pressure 0.0003 .-------.---r---r----,.-------.----, 0.0002 0.0001 !:!: 0 � -0.0001 -0.0002 -0.0003 )\ PL ---+-- / \ sc --- x -- @!1 __ � PL +SC --- -lt'--- ::-� � *** orig --- -0.0004 L.__---'---'----'--____..JL.__---'---' 4 �5 5 5.5 6 6.5 7 logT r 0. 004 .--------..--r---r-----,r-----r---, 0.003 0.002 0.001 0 -0.001 -0.002 PL ---+-- -0.003 '-- ---" '-- -'--..L.---' '-- ---'- --' 4 4.5 5 5.5 logT 6 6.5 7 0. 008 .--------..--r---r-----,r-----r---, 0.006 0.004 PL --+- sc -- X --- a. 0.002 � a. X «l 0 -0.002 -0.004 -0.006 -0.008 '-----'---'--..L.---''-----'---' 4 4.5 5 5.5 6 65 7 logT Figure 5.3: Compari son of relative difference (to OPAL model) quantiti es between the CEFF with PLPF and scattering terms(PL+ex), CEFF with only PLPF(P L),CEFF with only extra scattering terms(ex), and the original CEFF(orig). In Fig. 5.2, we showed the the relative difference w.r.t OPAL EOS from three differ ent models based on CEFF. First, the CEFF EOS with Planck- Larkin partition function as the internal free energy for hydro gen atom. Second, the CEFF EOS with only the scatteri ng state terms between free electrons and protons. Finally, the combination of the PLPF and 106 the scatter ing- states terms. From the figure one can see that the CEFF with only PLPF or only scatter ing- state terms deviate from OPAL with roughly same order of magnitude, but with opposite direction at around ln T = 4.5 rv 5. From this results we can make some observati ons: • The effects from the PLPF and the scattering terms are equally important while calculating the thermodynamical quantit ies. • The opposite effects from PLPF and scatter ing-state terms can combine and result in an modified CEFF which is much closer to OPAL than the original CEFF. As reference, in Fig. 5.3 the original CEFF EOS has also been plotted in the same graph. One can easily see the modified CEFF with both PLPF and scatter ing- state terms has been improved quite a bit than the original CEFF. By comparing with the original CEFF, there are also some points we can make: • Note that there is a peak happened at ln T = 6.5 consist ently in all four thermody namical quantiti es with the original CEFF and the CEFF with only scatter ing- state terms model. The reason is because of the ad hoc switch of pressure ionization from the original CEFF. In the CEFF with PLPF, since PLPF will be responsible for the ionization behavior at high temperatur e, (or say, before the pressure ionization occu rs), hence the pressure ionization switch was taken off within those models with PLPF. The results from PLPF implemented models agree with OPAL EOS very well at ln T = 6.5, without any non-sm ooth bumps or peaks. • When we look at the modified CEFF carefully, we can see there is still a small feature occurring around ln T = 4.5 rv 5. This feature is not very obvious in the 107 graph of pressure, but is rather clear as a small peak in the graph of r. In both X p and xr , the feature are broader than the one in r. Since the most cancellation between the PLPF effect and scatteri ng-st ate terms effect happe ns in this range, we believe this feature is the residual due to the cancellat ion. • Other than the above features or peaks, we can also see there are tails at the high temperature end in both modified CEFF or original CEFF. We believe this is due to the exchange terms that were included in OPAL but not been implemented in CEFF or modified CEFF. 5.1.4 Comparison between modified MHD and Modified CEFF We understand that since MHD and CEFF are developed independen tly, and with differ ent approach and structure, which the former is very sophisticated and complicated while the later is much simpler and effi cient, we should expect some intr insic difference between these two equation of state. From Fig. 5.4 we can see the original MHD and original CEFF was plotted with the modified MHD and modified CEFF. We can see the MHD and the CEFF, without the modificati on, behave differently, while MHD seems to be more wiggled and CEFF has a peak at T rv 106 5 , however they still have roughly the same accuracy. Once having PLPF and extra scattering terms implemented on both of the equation of state s, the results (CEFF(mod), MHD(mod)) converge nicely toward OPAL equation of state, in terms of the shape or the magnitude. Now we can feel confident that the same approach of modifying MHD EOS also works well on modifying CEFF EOS. 108 Xr 0.005 t 0 ,s -0.005 -0.01 -0.01 5 '------'---'---L------''------'----' 4 4.5 5 5.5 6 6.5 7 logT pressure 0. 00015 .........,..--r--..--------r--r-----r---, 0.0001 5e-05 0 r 0.006 .------.....,.,.,------.---r----,r------r---, 0.005 0.004 0.003 f::c 0.002 fo 0.001 0 -0.001 '' -0.002 ·� -0.003 '-- ----'- --'---L- -----' '-- ----'- ---' 4 4.5 5 5.5 6 65 7 logT 0. 004 .-- ----.. - ---.- --r---. r-- ---,jif;- ---, 0.002 0 a.. -5e-05 J:!- 0.. a: -0.0001 0.. -0.002 ,s <-0 -0.00015 fgtl!l) . -0.0002 �. \ /'ljJ #* -0.004 -0.00025 \ , ' � -)!( -0.006 -0.0003 -0.00035 '-- ----'- ----'-' """- ----'- - -'-- -----l. ---' -0.008 '------'---='---L------''------'----' 4 4.5 5 5.5 6 6.5 7 4 4.5 5 5.5 6 6.5 7 logT logT Figure 5.4: Compari son of relative difference (to OPAL model) quanti ties between the original MHD (MHD), the modified MHD with PLPF and extra scattering terms(MHD(mod )),original CEFF (CEFF), the modified CEFF with PLPF and extra scat tering term s(CEFF(mod)). The result of both the modified MHD and modified CEFF has been brought closer to OPAL significa ntly, and with almost same good qualities. 109 5.2 Hydrogen-Helium Plasma 5.2.1 MHD and modified MHD For the hyd rogen-helium mixture, one should consider Planck-Larkin partition function for hydrog en atom, helium atom, and helium plus ion, H,He,He+ as their internal parti tion function. Fig. 5.5 shows the MHD and MHD with PLPF and three extra scattering terms, in hydrogen-helium mixture. We can see the results reach similar quality as one with pure hydro gen plasma. As pointed out in the previous work [ 43] the scattering terms involving the new species, related with helium, has not been included. We know helium is a rather small fracti on comparing to hydrogen in solar plasma. However, in order to go further, we shall include those more extra scattering terms, and study how they impact the equation of state. 110 Xr r 0.006 0.006 0.005 MHO -+--- 0.004 0.004 MHD(sc+ PL) - - x--- 0.002 0.003 1- 0.002 t 0 f::c 0.001 """' X vO 0 vO -0.002 -0.001 -0.004 -0.002 -0.003 -0.006 -0.004 4 4.5 5 5.5 6 6.5 7 4 4.5 5 5.5 6 6.5 7 logT logT pressure X r 0.0001 0.004 5e-05 0.003 0.002 0 0.001 0 :it - a.. -5e-05 a. .:;::. -0.001 a: a. -0.002 vO -0.0001 X vO -0.003 -0.00015 -0.004 -0.0002 -0.005 -0.006 -0.00025 -0.007 4 4.5 5 5.5 6 6.5 7 4 4.5 5 5.5 6 6.5 7 logT logT Figure 5.5: Compari son of relative difference (to OPAL model) of quant ities between the original MHD (MHD) and the modified MHD with PLPF and extra scattering terms(sc+PL). The result of the modified MHD has been brought signifi cantly closer to OP AL. 5.2.2 More scattering terms for modified MHD The other possible scattering pairs involving helium should also be considered are : 111 • e-He ++ • He ++ - He ++ •p + -He ++ • e-He + •p + -He + • He + -He + • He + - He ++ In the program, free energy Fep is denoted as F5, Fee as F 6 , and Fep as F7. Here, we will need to add F8 - F 14 to include the above pair interact ions. Let's list the formula of those free energy terms: 112 Notice now, Aei is the thermal de Broglie wavelength divided by v'47r, i being different species list above. Aej = nj J2 f.-lej k B T' AD the Debye-Hiickel radius J EokBTBj e 2 (Ne + ZJ; N p ) , and the constant Dq2 , Dq3 are set to be 1.72, although they are not necessa rily to be the same as Dq. Fig. 5.6 and Fig. 5.7 show that how the equation of state changes while adding more and more extra scattering terms listed ab ove. The original model we use the modified MHD with Fep, Fee, F PP ' which are F5, F 6 , F7, respectively. Starting from here, we con sider the MHD with F5 - F8, F5 - F9, etc., until F5 - F 14 , therefore we can study how those terms aff ect the equation of state. 11 3 Xr r 0.002 0.001 0.001 0.0005 0 I- -0.001 0 t f::c """' X -0.002 vO -0.0005 vO -0.003 -0.004 -0.001 -0.005 -0.0 015 4 4.5 5 5.5 6 6.5 7 4 4.5 5 5.5 6 6.5 7 logT logT pressure Xp 7e-05 0.0025 6e-05 8--+- 8--+- , 5e-05 9 - -- )( --- 0.002 9 ---)(--- 10 --- -)IE -- - 10 --- 710--- 4e-05 11 - - ·- U- 0.0015 11 ·· ···- H·-- · a.. 3e-05 mhdex a. mhdex -+-: .:;::. a: 2e-05 a. 0.001 �� · vO 1e-05 X � A 7 vO 0 0.0005 -1 e-05 \ 0 -2e-05 ·� . '\:j; - ' -1 -3e-05 -0.0005 4 4.5 5 5.5 6 6.5 7 4 4.5 5 5.5 6 6.5 7 logT logT Figure 5.6: Compari son of relative difference (to OPAL model) quantiti es between the MHD with PLPF and scattering terms F5 - F7(MHDex), MHD with PLPF and scattering terms F5 -F8(8), MHD with PLPF and extra scattering terms F5 - F9(9), MHD with PLPF and extra scattering terms F5 - F10(10) and MHD with PLPF and extra scattering terms F5 - F11 (1 1 ) . 11 4 Xr 0. 002 .------,-----r--.--------,,---.----, 0.001 0 t -0.001 ,s -0.002 -0.003 -0.004 -0.005 .___ ___._ _ __._ _ _.___ ____. .__ __._ _ _. 4 4.5 5 5.5 6 6.5 7 logT pressure 7 e-05 .------,-----r--.--------,,---.----, 6e-05 5e-05 4e-05 3e-05 2e-05 1e-05 0 -1 e-05 -2e-05 11 --+- 12 --- )( --- � 13 ----)IE-- - 14 - - ·- U- mhdex -3e-05 .___ ___._ _ __._ _ _.___ ____. .__ __._ _ _. 4 4.5 5 5.5 6 6.5 7 logT r 0. 0008 .-----,----r--.--------,,----.----, 0.0006 0.0004 0.0002 t o vO -0.0002 -0.0004 -0.0006 t -0.0008 .___ ___._ _ _._ _ _.___ ____. .__ ___._ _ _. 4 �5 5 5.5 6 6.5 7 logT 0. 0025 .-----,----r--.------,,----.--..,. 0.002 0.0015 11 --+- 12 ---)(--- ' 13 --- 710 --- 14 ·· ··-G·-- mhdex --+--' J:!- 0.. 0.. 0.001 ,s 0.0005 0 -0.0005 '------'---'----"'--'----''-----'---' 4 4.5 5 5.5 6 6.5 7 logT Figure 5.7: Compari son of relative difference (to OPAL model) quantiti es between the MHD with PLPF and scattering terms F5 - F7(MHDex), MHD with PLPF and scattering terms F5 - F12(12), MHD with PLPF and extra scattering terms F5 - F13(13), MHD with PLPF and extra scattering terms F5 - F14(14). By studying the se graph s, we observe : • The most obvious changes in equation of state is after adding F8 and F10, which are FeHe ++ and Fp He ++ ' respectively. After adding FeHe ++ the kink happe ned at T rv 11 5 10 5 in r got enlarged, and we can say it makes the result worse. But one should bare in mind those terms will result in good equation of state via their cancellation with each other therefore the ultimate purpose is to have the the combi ned contributi on. After add FpHe ++ , the kink in r shrinks and the the bumps at T rv 10 5 3 in pressure, X T and X p are all flat tened. • After adding F 12, i.e., FpHe + , the kink in r grows again, and the rest of the quantiti es has not changed too much . • Fg, Fn , F 1 3 , F 14 , which are FHe ++ He ++ ' FeHe + ' FHe + He + , FHe + He ++ seem to not contribute to solar equation of state too much. 5.2.3 CEFF and modified CEFF Fig. 5.8 shows the CEFF model with PLPF and scatter ing-state terms between electrons and protons. Since now we are dealing with H-He mixture, so the PLPF should also apply to He atom and He+ ion. Unlike He+, which is just a hydrogen-like species, He requires more complicated energy levels. In Fig .5.8 we have shown the compa rison between three models: CEFF with PLPF and three extra scattering terms ,Fee,Fep,Fpp; CEFF with only PLPF and CEFF with only extra scattering terms. From this graph we can observe some featu res which has been revealed in previous work on modified MHD: 11 6 0.01 0.005 I- t 0 X vO -0.005 -0.01 -0.015 4 0.00025 0.0002 0.00015 0.0001 a.. 5e-05 a: 0 vO -5e-05 -0.0001 -0.00015 -0.0002 -0.00025 4 � � * ' t ' ' � \ i * 4.5 5 1li ,' \ � � 4.5 5 Xr CEFF --+- 'x. * * * PL --- \( -- )l( t --- * · -- ,' Xx � . 5.5 6 6.5 logT pressure 5.5 6 6.5 logT 7 7 0.0015 0.001 0.0005 f::c 0 """' vO -0.0005 -0.001 -0.0 015 -0.002 4 0.006 0.004 0.002 a. .:;::. a. 0 X vO -0.002 -0.004 -0.006 4 � r � * * CEFF � � � · i0' -Jt; , y-- ! !\'! I k � SC - � -- I , I r \ ;If ' , I : l ' I )!; ' ' / : \ I I !J( : ' ' Of • ·- . · ; ' k"- X ' : ' "' ' ' \.! 4.5 , !'>Jti 4.5 : � 5 5.5 logT X r 5 5.5 logT � : : . f * 6 6.5 CEFt L � sc -- - f-- i_ : * ' , 6 6.5 7 7 Figure 5.8: Compari son of relative difference (to OPAL model) quantiti es between the CEFF with PLPF and scattering terms(CEFF ), CEFF with only PLPF(P L),CEFF with only extra scattering terms (SC). • The effects of PLPF and scattering terms on the equation of state have similar order of magnitude, mostly opposite sign. 11 7 • The combined result shows the cancellation between the effects from PLPF and scattering ter ms, which make it very close to zero, hence the OPAL EOS. And the some featu res that are not revealed in MHD, but mainly the intr insic behavior ofCEFF: • Notice there is a sharp peak occurs at T rv 106 5 while only including the extra scat tering terms which is ca used by the pressure ionization switch in CEFF itself. Since it is a hand-made device rather than a natural physical outcome, therefore a small discontinuity happe ned there. Once the PLPF was included , this pressure ionization device has been switched off at the same time, therefore it can be smoothly carri ed over. 5.2.4 More scattering terms for modified CEFF Similar to what has been done with the modified MHD, we also brought those scattering terms related to helium nuclei and helium ion, and try to study to see if they have similar behavior as in MHD. These results have been plotted in Fig.5.9 and Fig.5 .10: 118 Xr r 0.001 0.001 0 0.0005 -0.001 1- 0 t -0.002 f::c """' X vO -0.0005 vO -0.003 -0.004 -0.001 -0.005 -0.0 015 4 4.5 5 5.5 6 6.5 7 4 4.5 5 5.5 6 6.5 7 logT logT pressure Xp Be-05 0.003 7e-05 0.0025 6e-05 5e-05 0.002 a.. 4e-05 � \ a. 0.0015 .:;::. a: 3e-05 a. vO 2e-05 X 0.001 \ vO 1e-05 0.0005 0 0 -1 e-05 -2e-05 -0.0005 4 4.5 5 5.5 6 6.5 7 4 4.5 6 6.5 7 logT logT Figure 5.9: Compari son of relative difference (to OPAL model) quantiti es between the CEFF with PLPF and scattering terms F5 - F7(ceffex ), CEFF with PLPF and scattering terms F5 - F8 (8), CEFF with PLPF and extra scattering terms F5 - F9(9), CEFF with PLPF and extra scattering terms F5 - F10(10) and CEFF with PLPF and extra scattering terms F5 - F11 (1 1 ). 11 9 Xr r 0.001 0.0006 0 0.0004 -0.001 0.0002 ·· ···· H·· ·· I- 0 t f::c -0.002 """' X vO -0.0002 vO -0.003 -0.0004 -0.004 -0.0006 -0.005 -0.0008 4 4.5 5 5.5 6 6.5 7 4 4.5 5 5.5 6 6.5 7 logT logT pressure Xp Be-05 0.003 7e-05 0.0025 6e-05 5e-05 0.002 a.. 4e-05 a. 0.0015 .:;::. a: 3e-05 a. vO 2e-05 X 0.001 vO 1e-05 0.0005 0 0 -1 e-05 -2e-05 -0.0005 4 4.5 5 5.5 6 6.5 7 4 4.5 5 5.5 6 6.5 7 logT logT Figure 5.10 : Compari son of relative diff erence (to OPAL model) quantiti es between the CEFF with PLPF and scattering terms F5 - F7(ceffex ), CEFF with PLPF and scattering terms F5 - F1 1(1 1), CEFF with PLPF and extra scattering terms F5 - F12(12 ), CEFF with PLPF and extra scattering terms F5 - F13(13) and CEFF with PLPF and extra scattering terms F5 - F14 (14). • The overall behavior is very similar to what has been demonstrated with modified MHD. 120 • After adding F8, which is FeHe ++ ' the pressure, xr,X p are improved a little bit. But forr, the kink at T rv 10 5 got enhanced, therefore against our wish. • After adding F10, which is FpHe ++ ' the kink in r has been reduced ever smaller than the point we started. (i.e. only having Fee, F PP and Fep) In addition, the feature occurred at around T"' 10 55 in pressure, xr and X p also got smoothed. • After adding F12, the FpHe+, instead of getting further improved, the wiggle in r runs bigger. And the other features in pressure, xr and X p are not improved signifi cantly . • Again, Fg, Fn , F13, F14, which are FHe ++ He ++ ' FeHe + ' FHe + He + ' FHe + He ++ seem to not contribute to solar equation of state too much. 5.2.5 Comparison between modified CEFF and CEFF(more), and modified MHD and MHD(more) Let's call the further modified MHD and CEFF for H-He plasma, MHD(more) and CEFF(more) since they have more scattering terms than previous modified version. According to the term by term analysis in previous section, we feel it is fine to only keep F8, F9,F10 as our "more scattering terms". There are two reasons: F11 - F14 either not doing much effect or make equation of state away from OP AL. To meet our purpose here, we should just choose the subset that gives us the best result. In Fig.5 .11, we could observe the incredibly resemblance between modified CEFF and modified MHD , and CEFF(more) and MHD(more). They behaves similarly before and after adding more scattering terms, yet still some slight differe nces: 121 Xr r 0.002 0.0008 0.001 0.0006 0 0.0004 I- -0.001 0.0002 t • f::c """' 0 X -0.002 vO vO -0.0002 -0.003 -0.0004 -0.004 -0.0006 -0.005 -0.0008 4 4.5 5 5.5 6 6.5 7 4 4.5 5 5.5 6 6.5 7 logT logT pressure Xp 8e-05 0.003 7e-05 CEFF(mod) --+- 0.0025 CEFF(mo d) --+- 6e-05 MHD(mod) --- x --- MHD(m od) - --x-- - 5e-05 CEFF(mr) ----)IE-- - 0.002 CEFF(mr) --- 710 --- 4e-05 MHD(mr) - -- u- MHD(mr) a.. a. 0.0015 3e-05 .:;::. a: 2e-05 a. vO X 0.001 vO 1e-05 tj:] 0 r9 0.0005 -1 e-05 \ -2e-05 0 -3e-05 -0.0005 4 4.5 5 5.5 6 6.5 7 4 4.5 5 5.5 6 6.5 7 logT logT Figure 5.11: Compari son of relative diff erence (to OPAL model) quantiti es between the modified MHD (MHD(mod)), the further modified MHD with PLPF and more extra scat tering terms(MHD(mr)),modified CEFF (CEFF(mod)), the further modified CEFF with PLPF and more extra scattering terms(CEFF(mr)). The result of both the MHD(more) and CEFF(more) has been brought even closer to OP AL, though not significa ntly. • In the subgraph of r, there are an very interesting transition point at T "-' 10 5 5 . While T < 10 55 , CEFF(mod) and MHD(mod) agree with each other perfectly 122 while CEFF(more) and MHD(more) also form another perfect match set. How ever, when the temperature increase to T > 10 55 , there are still obv iously two sets of data between with different pairing: now CEFF(mod) agrees with CEFF(mr), tailing downward, while MHD(mod) agrees with MHD(mr), tailing upward. From this observation we want to say the modificati ons (more scattering terms) seemed dominate at rather low temperature ln T < 5.5, so the equation of state will depend on either having those terms or not. But while temperature is higher, ln T > 5.5, the original behavior of MHD and CEFF dominate so the equation of state behave depend on either it is CEFF-based or MHD-based. • In the subgraph of pressure, we can see CEFF(mod)-CEFF(mr) and MHD(mod) MHD(mr) form two pairs, which have similar shape but with magnitude shifted. It looks like the magnitude still comes from the intr insic properties within CEFF or MHD, but the shape and effect of more scattering terms on the equation of state still behave similarly. To close this section, let's show the final results of the modified MHD and the modified CEFF with respect to the original MHD and original CEFF in Fig. 5.12 . We see that for both MHD EOS and CEFF EOS, the strategy of improving its quality by introducing the PLPF term and scattering terms works equally well. We now have a good tool in the form of a modified CEFF EOS, CEFF(more) , which can be applied in solar or stellar modeling. 123 [b] XT r 0.008 0.006 0.006 0.005 CEFF --+-- 0.004 MHO - -- x--- 0.004 0.003 CEFF(mr) ··· iii · - · """" 0.002 0.002 1'- x._ MHD(mr) ·· ·- o ·-- t '""" ";/.. >x 0 i: 0.001 / X X <0 0 <0 -0.002 -0.004 -0.001 -0.002 -0.006 -0.003 -0.008 -0.004 4 4.5 5 5.5 6 6.5 7 4 4.5 5 5.5 6 6.5 7 logT logT pressure Xp 0.00015 0.006 0.0001 0.004 5e-05 0.002 0 a. 0 !:!: � 0.. -5e-05 a. <0 X -0.002 -0.0001 <0 -0.00015 -0.004 -0.0002 -0.006 -0.00025 -0.008 4 4.5 5 5.5 6 6.5 7 4 4.5 5 5.5 6 6.5 7 logT logT Figure 5.12 : Compari son of relative diff erence (to OPAL model) quantiti es between the original MHD (MHD), the further modified MHD with PLPF and more extra scattering terms(MHD(mr)), original CEFF (CEFF), the further modified CEFF with PLPF and more extra scattering terms(CEFF(mr)). The result of both the MHD(more) and CEFF(more) are now very close to OPAL compared to the original version. 124 5.3 Outlook We feel confident about the performance of our CEFF-based OPAL emulator. We have not only succes sfully brought the experience from the modified MHD to the modified CEFF, beyond the prior H-only work, we have further developed and added more scattering terms necessa ry for the H-He plasma, both for the MHD and CEFF equations of state. This way, we have succeeded in matching OPAL the closest ever. We have accompli shed the ultimate goal of this project, and we plan to make our software availabl e on the website of the Journal of Com putational Physics. 125 Chapter 6 The Shear-driven Soret Ef fect In this chapter we discuss the solar tachocline problem and give the introduction to the Soret effect. We then review the reverse non-eq uilibrium molecular dynamics (RNEMD) simulation method, and introduce the setting of our simulation box. We then show our test run and the real simulation results. 6.1 The solar tachocline Thanks to hel ioseismology, we have a probe to the invisible internal structure of the Sun. As with human X-ray or ultrasound diagno sis, it is always fascinating to probe something inside that cannot be seen, beca use it can reveal unexpected th ings. At the same time it creates new probl ems for the curi ous human nature. The solar tachocline is exactly such a probl em. The internal rotational curve of the sun, inferred by heli oseismology, is shown in Fig. 6.1. Differential rotation appears in the convection zone, with a period at the equator of about 25 days and at the poles of more than 30 da ys. Hence the rotation rate decreas es monotonically from the equator toward the poles by about 30 %. It is only thanks to helioseismology that one can infer the internal rotati on. Before, the differ ential rotation was only observed on the surface velocity map. We can see that within the convection zone, the rotation curve is not too much different from the surface differential rotation 126 pattern, especially in the mid-latitudes area. Near the bottom of the convection zone, on the other hand, appear to have a sharp transition between the differential rotati on within the convection zone, and the nearly uniform rotation in the radiation zone. This rigid body-like rotation of the radiation zone and this sharp transition between radiation zone and convection zone had not been expected before helioseismology. This sharp transition in rotational velocity is called the solar tachocli ne. 127 450 425 N' 400 � 400 " 375 N ' a 350 350 325 b 300 �� 300 0.50 0.60 0.70 0.80 0.90 1.00 r/R Figure 6.1: Angular velocity profile in the solar interior inferred from helio seismology.[ 62] In panel( a), the rotational inversion is shown based on the subtractive optimally localized averaging (SOLA) technique, while in (b),the angular velocity is plotted as a function of radius for several selected latitudes.All inversions are based on data from the Michelson Doppler Imager (MDI) instrument aboard the SOHO spacecraft, averaged over 144 da ys. Recent seismic estimates [7, 1] of the tachocline location is about Tt rv 0.693 ± 0.003R0 near the equator which is below the base of the convection zone rb rv 0.713 ± 0.003R0, but may lie within the deeper convective overshoot region. While moving to higher latit udes, the location of the ta chocline also shif ts upward, to about Tt rv 128 0717 ± 0.003R0 at latitude 60°. Therefore the tachocline is obviously prolate. The shape of the tachocline is not like the shape of the base of the convection zone, which has not yet been detected with significant latitude dependence. The estimation of the thickne ss of the tachocline is highly depend on the way of defin ing the tachocline. Charbonneau et al. [7] characterize the thacocline transition in terms of an error function: 1 { [2(r-r t )]} f(r : rt , .6.t ) = 2 = 1 + e r f .6.t . (6.1) Their result yields the thick ness of tachocline is .6.t/ R0 = 0.0039±0.013 at the equator and .6.t/ � = 0.0042 ± 0.013 at 60° latitude. This suggests the tachocline broaden a bit on the higher latitude. On the contrary, Basu and Antia[1] suggested a more obvious widening of the tachocline with latitude. Started from .6.t rv 0.016R0 at the equator to .6.t rv 0.038R0 at 60° latitude. They also argued that this widening may not be smooth but with a sharp transition at mid latitude. Throughout most of the tach ocline, the vertical shear rate is about a order of magnitude larger that the latitudinal shear rate: dv¢,/dr rv ±1. 5 X 10 - 6s - 1 and r-1d vrpjd() rv 2 X 10 - 7 s-1. The solar tachocline is important beca use some solar-cycle dynamo theor ies include it as a crucial part that generate the solar dynamo and maintain this 11 -year solar cycle. First, note that there is no generally accepted theory of the solar cycle yet. Second, the question why the tachocline is stable and keeps itself within a narrow width that separate the convection zone from the radiation zone is also still an interesting and open question. There are several tachocline models that have been proposed: Spiegel & Zahn [61] argued that intense turb ulence could be generated by a shear inst ability of the differential rotat ion, 129 which if model with an anisotropic viscosity with enhanced horizontal components, then the associated angular momentum transport could prevent the radiative spreading of the shear layer, and then keep it thin. On the other hand, Gough & Mclnty re[29] proposed a weak meridional circulation driven by stresses from the overlying convection zone, and a core magnetic field prevents the shear layer from spreading into the deep interior. In order to confirm these model s, large-s cale simulation can be used. However, one has to deal with a system of highly-coupled, non-linear equations, and various numerical probl ems cau sed by the wide range of timescales inherent in the phy sics of the system, and all sorts of instabilities, physical and numerical. In our Sun, which is extremely com plicated, many possible effects are competing with each other, which ca uses this type of research to require lots of effo rt. However, up to this point, instead of jumping into this complicate enormous system, we tried to examine an important partial and iso lated prob lem, in the hope that in the future its solution would become part of general dynamo models. To illustrate, let's take a step back. We know that within the tachocli ne, the constitute material is the plasma with high density, high temperature, subject to the shearing velocity profile. Hence we would like to restrict ourselves to the question how this shearing plasma is functioning from a microscopic point of view. At the atomic level , the plasma is consisting of particles of different species, such as elect rons, protons, ions ... etc. In principle, they are subject to the Coulomb forc e, and their trajectories or behavior should depend on their most important properties, namely their charg e. However, since there are many different kinds of species, with differ ent masses, especially as far as electrons are concerned, compared to ions and nuclei, we deal with mass difference of the order of 10 3 . We can not stop wondering if can 130 this varied intrinsic properties of these particle ca use some effect in the shearing flow of the plasma? The most intri guing example is the phenomenon called Soret-Ludwig effect, also known as thermal diff usion. It is an effect purely cau sed by mass diffe rences of parti cles, where isotopes can separate from each other if subject to temperature gradient. The theoretical foundation of the effect was given by Onsager's work on weakly irrever sible thermodynamic systems. 6.2 Soret effect The Soret effect was first found by Ludwig in 1856[41], and it got further studied by Char les Soret in 1879[60]. He found that while the salt solution contained in a tub e, with temperature gradient was applied by holding the two end of the tube in differ ent temperature, then the salt will be more concentrated on the cold end than the hot end. This phenomenon of the salt preferentially diff use against the temperature gradient is then the well-known "Ludwig-So ret effect" , or thermal diff usion: the mass diff usion drive by a temperature gradient. Thermal diff usion has been propo sed to explain many differ ent phenomena of a wide range,e g., migration and concentration of DNA along tempera ture gradient[5], the isotopes separation in lava [37], or the constraining of the timing of the abrupt climate change of Younger Dryas from thermally fracti onated gases in polar ice[59]. Thermal diff usion also shows great potential in application of nanotechnology and biotechnology. 131 When thermal diff usion happe ns together with normal diff usion, which obeys Pick's law, in a binary system, a net mass flux is caused by a linear combination of the concen- trati on gradient and the temperature gradient:[ 47] (6.2) where 11 is the flux of species 1, p is average mass density, x1 the mole fraction of species 1, and mutual diff usion coefficient, D12 and thermal diff usion coefficient Dr. 8x1j8z is the concentration gradient of the species 1, and 8T / 8z is the temperature. While the non-eq uilibrium steady state was establi shed, the mass flux vani shes, then we have (6.3) (6.4) where Sr was defined as Soret coef ficient. Notice that Eq. (6.2) is the phenomenological equation, the determination of D12 and Dr will rely on the experiments or simulations. There have been many studies on the Soret effect using molecular dyna mics (MD) simulations. The ultimate goal for these types of research is to try to study the quantitative explanation of Soret effect of some specific system. Our attempt here, however, is differ ent. We did not dig into the details of the study on the Soret effect, but rather try to ask a question: if there is such an effect, similar to the Soret effe ct, but instead of coupling 132 between mass flux and heat flux, it has the coupling between the momentum flux and mass flux? In other words, the Eq.6.2 then can be modified as: (6.5) where Dm is then the cor responding coef ficient. We wish to see if this effect does exist by running the MD simulation, since it is a mature method for studying the Soret effect. Let's continue our discussion on the method of MD simulation in section 6.3. 6.3 MD simulation 6.3.1 Molecular dynamics Computer simulation acts as a bridge between microscopic length and time scale and mac roscopic real world. By having a large number of particles and the interaction between them, one can study their structure and properti es as an ensemble, via the theory of statis tics mecha nics. From this approach we hope to learn something new from the simulation as from the experiment in laboratory. The Molecular dyna mics (MD) simulation, in par ticular, can simulate the dynamical properties, like transport coefficient s, time-dependent responses to perturbat ions, rheological properties and spectra, which make it more ad van- tageous. The basic ideas of the MD: 133 • The Newton 's equation of motion: the molecular dyna mics simulation consists of the numerical integration of the classical equation of motion, i.e. the Newton 's sec- ond law. In a simple system, we just need to solve the equation: (6.6) • The pair-potential. The three-body or higher order interaction are usually neglected. therefore the U in the above equation can usually be written as: U(rN) = LLu(ri, rj) (6.7) i j > i 6.3.2 Integration method There are two simple integration schemes have been used in MD simulation, namely leapfrog and Verl et. They are equivalent algebraically, and both are accurate up to ,6.t 3 • In the fo llowing , we show the derivation of the Verlet formula. We start out from the Tayl or expansion of the coordinate variable: h 2 x(t +h ) = x(t) + hx(t) + ( - , )x(t) + O(h 3 ) 2. (6.8) where h ,6.t , x(t) being the velocity, x(t) being the acceleration which can be obtain from (6.6). Simila rly, we can also write down, h 2 x(t - h) = x(t) - hx(t) + ( - , )x(t) + O(h 3 ) 2. (6.9) 134 Use Eq.(6.8) - Eq.(6.9), and rearrange, we can have, x(t +h ) = 2x(t) - x(t - h) + h 2 i(t) + O(h 4 ) and the velocity can be obtai ned: x(t) = [x(t +h ) - x(t - h)] + O(h 2 ) 2h For the leapfrog meth od, we rewrite the Tayl or expansion as: x(t +h ) = x(t) + h[x(t) + (h/2)i(t)] + O(h 3 ) :i: ( t+h/ 2 ) If consider (6.10) (6.11) (6.12) x(t - h/2) = x(t) - (h/2)x(t) , x(t + h/2) = x(t) + (h/2)x(t) (6.13) one can have the relat ion: x(t + h/2) = x(t - h/2) + hx(t) Hence, the leapfrog integration formulae are: x(t + h/2) = x(t - h/2) + hx(t) x(t +h ) = x(t) + hx(t + h/2) (6.14) (6.15) (6.16) 135 6.3.3 Periodic Boundary Conditions If there is a simulation box with physical wall constructing by the parti cles with fixed position s, the ratio of the number of the particles on the surface to the internal particles will rv N'/; -(3 . In real world the surface effect can be neglected since typically Nm rv 10 2 3 and the surface atom are merely rv 10 - 8 of Nm. However for a MD number, Nm = 1000, then almost 500 particles will be on the surfac e. If we are interested in a infinite system and the surface effect is not our concern, then the periodic boundary condition should be applied. To introduce the periodic boundary condition is like to have infinite, space-filling array of the replicas of the original simulation box. The effe cts from the periodic boundary are essentially twofol d: • While a parti cle leaves the simulation box there will be a particle entering from the opposite surface immediately. • While calculating the pair potential, a particle interacting with the particles from the other side of the boundary will be considered as they are adjacent-a wraparound effect. In this way it is obvious that there is no physical boundaries and is possible for one to model a homogeneous system. To implement a periodic boundary, for exampl e, a simulation box with size Lx in x direction, and the x-coo rdinate of particle i, rix· are allowed to be between - Lx/2 and Lx /2. The test of the in integrating the equation of motion is • if rix � Lx/2, replace it by rix - Lx • if rix :s; - Lx/2 replace it by rix + Lx 136 The test of the calculating the pair interaction is very similar: • if T'ijx 2' Lx / 2 , replace it by T'ijx - Lx • if T'ijx < - Lx/2 replace it by T'ijx + Lx 6.3.4 Equation of motion The Lennar d-Jones(LJ) potential, originally proposed for the liquid argon, has two prin- cipal features of an interatomic force. One of them is the resistance of the compre ssion, hence the interaction repels at close range. The other is the induced dipole interaction which will bind the atom into solid or liquid state, therefore the atom pair will have long- range attracti on. For a pair of atom i and j located at ri and rj, the potential energy is { 4E [ (_!!___) 1 2 - (_!!___) 6 ] r· · < r ( ) r;j r;j ' t} c u lij = 0, lij > T' c where T'ij h - rj[· The parameter E governs the strength of the interact ions and O" defines the length scale. The strong repulsion in short range comes from the quantum effect that the cores of atoms repel each other; while the long-range attractive tail repre sents the van der Waal s interaction due to electron correlat ions. The potential can be further simplified by ignoring the attractive tail and have + E, r-- < r = 2 1 1 6 0" t} c 137 The choice of rc makes u (rc) = 0. The force corresponding to u (r) is f = -grad u (r) So the force of atom j act on atom i is while rij :s; r c• or otherwi se zero. The equation of motion follows form Newton 's second law, 6.3.5 Dimensionless unit Nm m ri = 1 i = 2:::: 1 ij j=l,jfi (6.17) (6.18) (6.19) For MD simulation based on LJ form. (Eq.(6.1 8)), the most suitable unit is by choosing O", m and E as the unit of length, mass and energy. Here we need to make the replacement: length: r --+ rO" energy: e --+ eE Under this new definitions of uni ts, equation of motion can be written as: (6.20) (6.21) (6.22) (6.23) 138 h = 48 2::: ( r ij 14 - � r ij B ) r i j , j( f i ) and the dimension less kinetic and potential ener gies, per atom, are where vi is the velocity. 6.4 Reverse Non-equilibrium Molecular dynamics (6.24) (6.25) (6.26) Reverse non-e quilibrium molecular dynamics (RNEMD) is a branch of NEMD developed by Miiller-Plat he.[45, 46, 47] The name "rever se" comes from the reversal of the experi- mental ca use- and-effect picture while calculating the transport coef ficient. What that mean is, for example, in calculating viscosities, the effect : the momentum flux or stress, is imposed, while the cause, the velocity gradient or shear rate, is measured from the simu- lation. The advantage of this method is that the total energy and linear moment conserve, therefore no external thermostatting is needed. The resulting raw data is robu st and rapidly conve rging. This method has been tested on the calculation of the shear viscosi ty, the ther mal conductivity and the Soret coefficient (thermal diff usion). 139 6.4.1 Method:calculating shear viscosity as an example Consider a system under shear field �, which is the gradient of one component of the average fluid velocity, say, the x direct ion, with respect to the perpendicular direct ion, say the z direct ion. Then the x component of the momentum Px will transport in the z direction, denoted as the momentum flux jz(p x)· The connection between jz (px) and a;;; is the defined as the viscosity coef ficient. . ( ) Bvx J z Px = - rt Bz (6.27) In RNEMD, the momentum flux was imposed on the system in a unphysical way. Suppo sed we would like to have the shear look like on the top of the simulation box, z = Lz/2 and the bottom z = -Lz/2 with positive x- component velocity, +vx, and the negative velocity -vx in the middle, at z = 0. The way to obtain this velocity profile is by imposing the unphysical momentum transfer via swapping the x- component of the momentum of the two selected parti cles. One should define the slabs on top, bottom and middle of the simulation. Then, by scanning the particles in the slab, one can always find the particles with the most undesired Vx· For instance, since we hope the top slab will have +vx, then we will select the particle with the most negative Vx within the slab. The bottom slab is simply the replica of the top slab via the periodic boundary condition, so we don't have to select any particles from it. Then we choose the most positive Vx in the middle since we hope it can move toward - x direction. Then these two selected particles trade their Px, to complete this unphy sical momentum transfer. 140 By imposing this swapping of Px, the system will eventually develop the physical momentum transfer via the interaction between particles. This proce ss is summarize in Fig.6 .2 . Lz Momentum Flux 1-------1 (physical) Figure 6.2: Simulation box. The reverse non-e quilibrium molecular dynamics bri ngs the momentum from the middle slab to the top and bottom slab by exchanging the momentum of the parti cles. Then the momentum flows back through the fluid by fri ction, thereby ca using a gradient in the velocity in x-direction 141 6.4.2 Methods for thermal conduction and thermal diff usion With the same idea from the previous example of finding viscosi ty, the RNEMD can be easily generated to study thermal conduction and thermal diff usion • For thermal conduction, we should define the slabs with high temperature and low temperature. By finding the particle with most undesired temperature within each slab, and then switch their velocity. (It is like finding the coldest atom in the hot slab and exchange it with the hottest atom in the cold slab, if the former is still hotter than the later.) • For Ludwig-Soret effect, or thermal diff usion, the procedure is very similar to the thermal conduction case. However, since the system now is mixture of two kinds of particles with differ ent mas ses, hence the excha nge of the velocity should only occur between species with same mass. 6.5 Our RNEMD for sheared mixture system We intended to see if the shear-driven Soret effect would happen in the MD simulation, since it was never been reported in experiment. Our strategy is to construct a system and condition which can observe the normal Soret effect, and then replace the temperature gradient with a velocity gradient (as subj ected to a shear force) in a way was described in the RNEMD algorithm. Since this is a rather fundamental research, so we are not intending to apply the realistic parameter as in solar tac hocline, since after all the solar plasma would be more complicated Coulomb system. In this simulation we only dealt with the Lennard- Jones potential, so it will be like the system of noble gas and its isotope. 142 Our simulation box is exactly as depicted in Fig.6.2. Hence the shear was applied both on top and bottom, also in the middle of the simulation box. To analyze our data, we simply slice our simulation box into 50 thin slabs in z- directi on. We calculate the major quantities within each slab, like temperature, velocity and density, and then plot them against z-axis. 6.6 Program check To make sure our program will work, we had to run some test simulations to check if it can produce the anticipated results. Therefore we would have liked to reproduce the well know Soret effect in our simulation, and see the cases with and without the mass difference between species. We do realize that the energy will be drif ted, and this is inevitable since in this algo rithm there is no thermostat applied. We will need to choose the appropriate frequency of swapping the particle for unphysical heat flux, denoted as Nexch.(The parti cle will exchange every Nexch step coun ts.) We see it is obvious that the system with more fre quent particle change will be more unstable. The energy drift in Nexch = 5 case is much rapid that the Nexch =50 case. Therefore we will use Nexch =50 in our simulation. 143 to the algorithm we only exchange the x- component of the velocity in the sheared case, but have to change all (x-,y-,z-) components of the velocity in the system with heat flux. 200000 400000 600000 800000 # of step counts 1e+06 KE X 1.2e+06 * D 1.4e+06 Figure 6.4: The energy drift of the system with unphysical heat flux and the system with unphysical momentum flux. The total energy(totE-shear) of the system with momentum flux and the corresponding kinetic energy(KE); the total energy(totE-shear) of the system with heat flux and the corresponding kinetic energy(KE) are shown . The direction of the energy shift in these two system is opposite. 145 In Fig.6.5 we have the result from a test run, in the condition that are the same as in the real simulation, but with the mass of two species set the same. This is equivalent to say there is just one species in the system. From the graph of temperature profile, we see the temperature gradient was developed during this test run but there is no separation in the density profile. Also, the velocity profile of the velocity in x- component is randomly distributed around zero, which is what we expected as there is no translation of the system as a whole. This confirms that this system is reliable, since it doesn't create any fake separation of the density profile due to the numerical effects. Hence, if we observed the separation of the density profile in the later simulation, we should conclude it is from the mass difference effect. 146 • density = 0.8 • temperature= 0.8 In Fig.6.6 - Fig.6.8 we present the results for the simulation on step counts of 62000, 184000 and 482000 in velocity profile, temperature profile and particle density pro file. In these graph, the green line(with subscript "2", like T2, V2, M 2) stands for the lighter particles while the red line stands for the heavier particles. Their mass ratio is Mlight/ Mheavy = 1 1 0• In Fig.6.6 the density profile has still been being developing, but the temperature gradient has already been developed quite clear. It seems that the density profile obtains it final shape around step counts = 482000. Fig.6.9 shows the process of the density profile development. The Soret effect is confirmed existing in our set-ups. We will just repeat the same test for another condition in the next section. 148 T-profile 9 --�---..--,--�-r�--r-.-�-. 8.5 8 7.5 7 I- 6.5 6 5.5 5 4.5 4 ��-L�--��--�� 0 0.1 0.2 0.3 0.4 0.5 0.6 0. 7 0.8 0.9 z-direction V-profile 0.06 .--�---..-----.--r---.------.--,--r-.----, 0.04 0.02 0 -0.02 -0.04 '"*"' '""'"'""" -0.06 -0.08 -0.1 -0.12 -0.14 -0. 1 6 L--'-------'------'--..__--'----'------'--.J...._--'--_J 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 z-direction density-profile 1 .3 .----.-------.-----,.-�---..--,--�-r�-----, 1.2 1.1 0.9 0.8 0.7 0.6 '----'----'-----'��-L�--��-----' 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 z-direction Figure 6.6: The result of Soret diffusion simulation, total step counts is 62000. The density profiles of different masses are still developing. 149 T-profile density-profile 9.5 1.4 9 1.3 8.5 1.2 8 1.1 7.5 I- 7 :::2: 6.5 0.9 6 0.8 5.5 0.7 5 0.6 4.5 0.5 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 z-direction z-direction V-profile 0.04 0.02 0 -0.02 5 -0.04 -0.06 -0.08 -0.1 -0.12 -0.14 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 z-direction Figure 6.7: The result of Soret diffusion simulation, total step counts accumulated to 184000. The density profiles of different masses are almost saturated. 150 I- T-profile 11 T -+- 10 7 9 8 7 6 5 ��--��--��--�� 0 0.1 0.2 0.3 0.4 0.5 0.6 0. 7 0.8 0.9 1 z-direction V-profile 0. 06 .----.----,---.----.----,-----,.--.-----,----,-,,-----, 0.04 0.02 0 -0.02 -0.04 -0.06 -0.08 -0.1 -0. 1 2 L__ .:.L_ __J_ __ L__ _L___J_ ___.J L__L_ ----'- ___.J L___j 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 z-direction density-profile 1 . 5 .------.-----,---r-----.----.------.---.-----.-----,----, 1.4 1.3 1.2 1.1 0.9 0.8 0.7 M -+ M2 - - x --- 0.6 '----'------'----'---'-�---'------'--�'--�� 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 z-direction Figure 6.8: The result of Soret diffusion simulation, total step counts accumulated to 482000. The density profiles of different masses are now stable. 151 >- � c Q) "0 density-profile 1.5 .-----r-----r-----�----�----,-----,-----�-----.-----r-----, 1.2 1.1 0.9 0.8 0.7 0.6 0.5 L__ ____ ...__ ____ ..__ ____ .J...._ ____ ..J...._ ____ _._ ____ __._ ____ _J_ ____ __J_ ____ _�.._ ____ _J 0 0.1 0.2 0.3 0.4 0.5 z-direction 0.6 0.7 0.8 0.9 Figure 6.9: The result of Soret diffusion simulation, just compare the density profile through different time point of the simulation. We noticed that the density profile does not change much since time step counts= 482000. 6.8 So ret effect in liquid phase We have run this simulation in liquid phase. The corresponding parameters are: • � T= 0.002 152 • density = 0.8 • temperature = 1.2 In Fig.6.1 0 we present the results for the simulation of 1028000 step counts in velocity profile, temperature profile and particle density profile. In this graph, the green line stands for the lighter particles while the red line stands for the heavier particles. Their mass ratio is Mlight/ Mheavy = 1 1 0• The thermal diffusion again was demonstrated in the liquid phase. The change of the density profile through the simulation is shown in Fig.6.11. We see very similar process as the example in gas phase. 153 12 11 10 I- 9 8 7 0.05 0 -0.05 5 -0.1 -0.15 -0.2 T-profile 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 z-direction V-profile 1 n x 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 z-direction 1 density-profile 1 .5 .----.-------.-----,.--.------.------.--.----.------.------, 1.4 1.3 1.2 1.1 0.9 0.8 0.7 0. 6 <=---'----'-----''----'----'-----'--'---'-----'------' 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 z-direction Figure 6.10: The result of thermal diffusion simulation, total step counts accumulated to 1028000. The density profiles has not clearly developed. 154 >- � c Q) '0 density-profile 1.5.----..----,----�----�----�-----r-----r----�----�----· 1.4 1.3 1.2 1.1 0.9 0.8 0.7 M-122 --+-- •• .'� . . M2-122 - - x--- M-424 _ __ * __ _ ;r· '+i. e ll : M2-424 ·-··- H---.. • : . ' if_ 'e .. - M-726 ....,. i�· ' :;{ M2-726 .:1\ ,lili M-1 028 -- - •- - �- * \. M2-1028 - -A � \ / � � /!-U· El · t;J i A A ·� � _.0 I \_ &__� · · · Gl · .,; / · . .<J .. .. • - ,' I .EJ . .'1 '!\ Y.. . ,1 ' - · . . '' .· � r· : ' e I" ' J 'If, x, \ x / � \ i'J 1 \� · · . · !I' � i11 >< ' x GJ\ / r: � f \\ ·. •� ip.- -i D-· ' \ ";,. .0 \ '-rl q \ : ' " "" · .' / x. ':f X X ;t 'X' -� ' I \ ' (J I ' , ::�. � - '* � ,fi.i ! "' , , • G: , ll\ ,_ •• Y"' X -� r- x-;J:_ k .X \ i il � . ..: - � \ : . . Wi / a· 1 X :>( , -\ - ·· · • ';1:. '- t; , ! ·. I i ll' - �� X.. ' \ \ [j . IZi j. ,. :. - z!, · · · s '-'-fi! , \ �--- -x , /'"! / r : ' · ' , :t / . .... .. li' tJ " · I ' X fi \ ! I � .· . ' 1, '11( ; \ ': , \ \ ;/ \ , x-''\ ;P] ·¥- .JJ .,_ • \ .. . . "'J ,.. ', �- "- ! ; '- I '. ,e . ' · - : · \ · ·� I I\ 0 ! & jt .. * . •• . �\¥.� -- \. . . . / � A l'l\ ! \ / d, � ';II ''<' )/;, )./ .. ·_· .} • \ ':.; " ' · . ' ' ·, !ii-.r ' . ' . � ,.. .. B- el \. · --g: � \ / �-· \ • • .- tJ / \ � · �- • l'l . I lli , \ / . � •:: 0.6�--��--�----�----�----�-----L-----L-----L----�----� 0 0.1 0.2 0.3 0.4 0.5 z-direction 0.6 0.7 0.8 0.9 Figure 6.11: The result of Soret diffusion simulation, just compare the density profile through different time point of the simulation. We noticed that the density profile does not change much since time step counts= 424000. 6.9 Shear-driven So ret effect in gas phase Now we come to the simulation that we were really interested in. We have replaced the temperature flux with the momentum flux in x- direction, with the system of mixture of 155 two species with different masses. We would like to see, if the shear-driven Soret effect can be observed in such kind of simulation? This simulation is under this condition: • .6. T = 0.001 • density = 0.8 • temperature= 0.8 Fig.6.12 to Fig.6.18 show the results of the simulation at the time steps at 62000, 124000,310000,920000,2706000, 12610000 and 36016000, respectively. We should say, however, compared to the heat-driven Soret effect, the density profiles was not obviously separated. However, while simulation runs, not only the velocity gradient was established, but also a peculiar temperature profile emerging. This temperature looks like a "M" shape in our graph. Note that it is symmetric about z- direction, unlike the velocity profile. In other words, the average of the x- velocity is positive or negative doesn't matter. The temperature peaks at the average velocity closed to zero, the then slightly drops while the magnitude of the velocity grows. The temperature gradient is actually expected in this kind of simulation [ 46], and the reason for it is due to the momentum exchange. In this algorithm, it always picks the particle with the largest peculiar velocity ui =vi � ( v) against the flow. After exchange, this particle will be following the flow, i.e., the peculiar velocity become smaller, hence the temperature drops since the definition of the temperature is through the peculiar velocity. This mechanism also serves as an internal thermostat that drains the viscous heating in the sheared flow. We should bear in mind this temperature gradient is merely due to the numerical scheme but not related to any real physical mechanism. 156 Although the temperature gradient is artificial however the separation in the density profile of two species with different masses does reveal this temperature gradient. In order to see if there is real shear-driven density profile separation, we should exclude this features first. As a result, we should report we did not detect a pure shear-driven Soret effect. The pattern of the separation in the density profile is small and consistent with the temperature gradient, hence due to the normal Soret effect. We did not see the definite separation of the density profile that reveals the the definite velocity gradient, as in the case in the normal Soret effect. 157 I- 6.4 6.2 6 5.8 5.6 5.4 5.2 5 4.8 1 0.8 0.6 0.4 0.2 5 0 -0.2 -0.4 -0.6 -0.8 -1 T-profile /{ I I I \ I I I I t I 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 z-direction V-profile 0 0.1 0.2 0.3 0.4 0.5 0.6 0. 7 0.8 0.9 1 z-direction density-profile 1 . 15 ,----y--,------.-.-----.------r-----r--.------.----, 1.1 1.05 0.95 0.9 0.85 0.8 .__......__.____.._..__......__._____.._..___.____, 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 z-direction Figure 6.12: The result of shear-driven Soret diffusion simulation, total step counts accu mulated to 62000. The density profiles has not clearly developed. 158 T-profile 6.4 6.2 6 5.8 I- 5.6 5.4 5.2 5 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 z-direction V-profile 1 0.8 v --+-- ,· 0.6 V2 ---X-- 0.4 0.2 5 0 -0.2 -0.4 -0.6 -0.8 -1 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 z-direction :::2: 1 density-profile 1 .2 ,----y--,------.-.--------r-----r--.------.----, 1.15 1.1 1.05 0.95 0.9 0.85 M--+- M2 -- x --- 0.8 .__......___.____._..__......__._____.._..___.____, 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 z-direction Figure 6.13: The result of shear-driven Soret diffusion simulation, total step counts accu mulated to 124000. The density profiles has not clearly developed. 159 T-profile 6.2 6 5.8 5.6 I- 5.4 5.2 5 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 z-direction V-profile 1 0.8 0.6 0.4 0.2 5 0 -0.2 -0.4 -0.6 -0.8 -1 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 z-direction :::2: 1 density-profile 1 .2 ,----y--,------.-.-----.------r-----r--.------.----, 1.15 1.1 1.05 0.95 0.9 M -t- M2 -- x --- 0. 85 .__......___.___.._..____,___.____.._..___.____, 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 z-direction Figure 6.14: The result of shear-driven Soret diffusion simulation, total step counts accu mulated to 310000. The density profiles has not clearly developed. 160 T-profile 6.4 6.2 6 5.8 I- 5.6 5.4 5.2 5 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 z-direction V-profile 1 0.8 0.6 0.4 0.2 5 0 -0.2 -0.4 -0.6 -0.8 -1 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 z-direction 1.3 1.25 1.2 1.15 1.1 :::2: 1.05 1 0.95 0.9 0.85 0.8 0.75 1 0 density-profile M -t- M2 -- x --- 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 z-direction Figure 6.15: The result of shear-driven Soret diffusion simulation, total step counts accu mulated to 920000. The density profiles has not clearly developed. 161 T-profile 6.4 6.2 6 5.8 5.6 I- 5.2 5 4.8 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 z-direction V-profile 1 0.8 0.6 0.4 0.2 5 0 -0.2 -0.4 -0.6 -0.8 -1 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 z-direction 1.25 1.2 1.15 1.1 1.05 :::2: 1 0.95 0.9 0.85 0.8 0.75 1 0 density-profile 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 z-direction Figure 6.16: The result of shear-driven Soret diffusion simulation, total step counts accu mulated to 2706000. The density profiles has not clearly developed. 162 T-profile density-profile 6.2 1.3 6 M-+-- 5.8 1.2 M2 -- x --- 5.6 1.1 5.4 :::2: 1- 5.2 5 0.9 4.8 0.8 4.6 0.7 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 z-direction z-direction V-profile 1 0.8 0.6 0.4 0.2 5 0 -0.2 -0.4 -0.6 -0.8 -1 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 z-direction Figure 6.17: The result of shear-driven Soret diffusion simulation, total step counts accu mulated to 12610000. The density profiles has not clearly developed. 163 T-profile 5.4 5.2 5 4.8 I- 4.6 4.4 4.2 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 z-direction V-profile 1 0.8 v --+-- 0.6 7 0.4 0.2 5 0 -0.2 -0.4 -0.6 -0.8 -1 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 z-direction :::2: 1 1 density-profile 1 .2 .---..------.------,---r---.------r----r-----,.-----.,---, 1.15 1.1 1.05 0.95 0.9 0.85 0.8 ,__...___.___.___.____.____.____.____.,.___...______. 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 z-direction Figure 6.18: The result of shear-driven Soret diffusion simulation, total step counts accu mulated to 36016000. The density profiles has not clearly developed. 6.10 Shear-driven Soret effect in liquid phase We have also simulated the sheared system in different system, and try to see if similar behavior may be retained. This simulation is under this condition: • � T= 0.001 164 • density = 0.8 • temperature = 1.2 We repeated the simulation as in section 6.9. The similar behavior also be retained. See Fig.6.19 - Fig.6.22 T-profile density-profile 9 1.25 8.8 1.2 8.6 1.15 8.4 1.1 8.2 1.05 1- 8 :::2: 1 7.8 7.6 0.95 ' 7.4 � 0.9 7.2 0.85 7 0.8 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 z-direction z-direction V-profile 1.5 0.5 X 0 > -0.5 -1 -1.5 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 z-direction Figure 6.19: The result of thermal diffusion simulation, total step counts accumulated to 902000. The density profiles has not clearly developed. 165 T-profile 9.5 9 8.5 I- 8 7.5 7 6.5 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 z-direction V-profile 1.5 0.5 5 0 / -0.5 -1 -1.5 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 z-direction 1.25 1.2 1.15 1.1 1.05 :::2: 1 0.95 0.9 0.85 0.8 0.75 1 ' 0 density-profile M -t- M2 -- x --- 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 z-direction Figure 6.20: The result of thermal diffusion simulation, total step counts accumulated to 1814000. The density profiles has not clearly developed. 166 T-profile density-profile 9 1.25 1.2 M -+-- 8.5 1.15 M2 -- x --- 1.1 8 1.05 1- :::2: 7.5 0.95 7 0.9 0.85 0.8 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 0 z-direction z-direction V-profile 1.5 0.5 5 0 -0.5 -1 -1.5 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 z-direction Figure 6.21: The result of thermal diffusion simulation, total step counts accumulated to 2706000. The density profiles has not clearly developed. 167 T-profile 8.8 8.6 8.4 8.2 8 I- 7.8 7.6 7.4 { 7.2 7 6.8 6.6 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 z-direction V-profile 1.5 0.5 " 5 0 -0.5 . -1 -1.5 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 z-direction 1 1 density-profile 1 . 25 .---,----.----.,.-----,-----.-.---,----.---,----, 1.2 1.15 1.1 1.05 0.95 �fx.f 0.9 If 1/ \ ,,l * M -t- M2 --x --- 0.85 * 0.8 ' --- ' - - -- '- ----'- ....:.....L. - -'-' - - - ' - - -- '-----'- --- -' 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 z-direction Figure 6.22: The result of thermal diffusion simulation, total step counts accumulated to 11708000. The density profiles has not clearly developed. 6.11 Outlook The Shear-driven Soret effect does not exist as we had assumed in Eq.6.5. However, a shear-induced temperature gradient has been found in the simulation, and this temperature gradient does cause the separation between species with different masses. However, this 168 temperature gradient might have well been caused by the numerical scheme adopted, and we have therefore not yet found conclusive evidence of a real physical effect. 169 Chapter 7 Conclusion In this thesis we have been dealing with two different types of projects, both related to the Sun, but with very different physical regimes. As the title of the thesis implies, one of these projects is related to the physics in thermal equilibrium, more precisely to models of the equation of state (EOS). The other project is related to the transport phenomenon of the steady state of fluid, which is not in thermal equilibrium. At first, such different regimes might appear contradictory since they refer to the same object, our Sun. However, in practice, in astrophysics the concept of the so-called local thermal equilibrium (LTE) is widely accepted. Depending on the degree of accuracy, it is adequate. Under the right conditions, within certain regions, the distribution function of the velocity of particles will be Maxwellian to a very good degree. In that case one can locally assume thermal equilibrium, even if there is a temperature gradient on larger scales. In the part dedicated to thermal equilibrium, we have continued previous work on the MHD EOS. Features that were missing in the MHD EOS have been borrowed from its major competing formalism, the OPAL EOS. The OPAL EOS is known to be the best EOS so far when compared with solar observations, however its source has never been released and its results can only be obtained from large-scale tables with fixed chemical composition, as released by the OPAL group. The prior developed modified MHD does perform as an OPAL emulator, and provide the community the flexibility to test their theories. However, although MHD is free to the public, by itself it is not a user-friendly 170 tool. In our present project we aimed at providing a user-friendly tool that can be easily implemented and given any kind of chemical composition freely and easily. We carried the experiences from the modified MHD EOS to a more flexible formalism, CEFF EOS. CEFF EOS is itself already a good EOS and has been widely used. This makes upgrading it to an OPAL emulator especially desirable. We have succeeded in successfully updating CEFF to be as good as the prior tool, the modified MHD, but this time we have also added more physics, because the current work is valid for H and He, in contrast to the prior H only formalism. We therefore had to develop new terms related to helium. We have added those missing terms, tested them first on the MHD EOS and then carried them to CEFF. We have shown that our effort has brought the emulator closer to OPAL than ever before. In the part dedicated to a regime outside of thermal equilibrium, we have used a MD simulation to test if there exists a shear-induced analogy of the Soret effect. The relevance of such an effect would be for the solar tachocline, which is defined as the sharp transition between the differential rotation in the convection zone and the solid-sphere rotation in the radiative zone beneath. We have adopted a "Reverse-NEMD" method for this simulation, and we have tested both for a normal and a shear-driven Soret effect. As a result, we have not seen a shear-driven Soret effect in analogy to the normal Soret effect. We do observe a temperature gradient, and the density profile of different mass does separate under this temperature gradient, but the temperature gradient could have been caused by the numerical scheme and not by physics. Therefore, as this point, we cannot positively report on a shear-induced effect of separation of particles according to their mass. 171 Bibliography [1] S. Basu and H. M. Antia. Changes in Solar Dynamics from 1995 to 2002. Apl., 585:553-565, March 2003. [2] S. Basu, W. Dappen, and A. Nayfonov. Helioseismic Analysis of the Hydrogen Partition Function in the Solar Interior. Apl., 518:985-993, June 1999. [3] G. Berthomieu, J. Provost, A. Rocca, A. J. Cooper, D. 0. Gough, and Y. Osaki. Sensitivity of five minute eigenfrequencies to the structure of the sun. In H. A. Hill & W. A. Dziembowski, editor, Nonradial and Nonlinear Stellar Pulsation, volume 125 of Lecture Notes in Physics, Berlin Springer Verlag, pages 307-312, 1980. [4] E. Beth and G. Uhlenbeck. The quantum theory of the non-ideal gas. II. Behaviour at low temperatures. Physica, 4:915-924, October 1937. [5] D. Braun and A. Libchaber. PERSPECTIVE: Thermal force approach to molecular evolution. Physical Biology, 1:1P, March 2004. [6] S. Chandrasekhar. An introduction to the study of stellar structure. 1957. [7] P. Charbonneau, M. Dikpati, and P. A. Gilman. Stability of the Solar Latitudinal Dif ferential Rotation Inferred from Helioseismic Data. Apl., 526:523-537, November 1999. [8] J. Christensen-Dalsgaard. On solar models and their periods of oscillation. Monthly Notices of the Royal Astronomical Society, 199:735- 761, May 1982. [9] J. Christensen-Dalsgaard. Helioseismology. Rev. Mod. Phys., 74:1073, March 2003. 172 [10] J. Christensen-Dalsgaar d and W. Daeppen. Solar oscillations and the equation of state. Astron. Astrophys. Rev., 4:267-361, 1992. [11] J. Christensen-Dalsgaard, D. Gough, and J. Toomre. Seismology of the sun. Science, 229:923-93 1, September 1985. [12] J. Christensen-Dalsgaard, J. Schou, and M. J. Thompson. A comparison of methods for inverting helioseismic data. Mon. Not. R. astr. Soc., 242:353-369, February 1990. [13] W. Daeppen. An analytical version of the free-energy-minimization method for the equation of state of stellar plasmas. A. & A., 91:212-220, November 1980. [14] W. Dappen, L. Anderson, and D. Mihalas. Statistical mechanics of partially ion ized stellar plasma - The Planck-Larkin partition function, polarization shifts, and simulations of optical spectra. ApJ., 319:195-206, August 1987. [15] W. Dappen and D. Mao. A smooth equation of state for solar and stellar abundance determinations. Journal of Physics A Mathematical General, 42(21):2 14006, May 2009. [16] W. Dappen, D. Mihalas, D. G. Hummer, and B. W. Mihalas. The equation of state for stellar envelopes. III - Thermodynamic quantities. ApJ., 332:261-2 70, September 1988. [17] H. E. Dewitt. Statistical Mechanics of High- Temperature Quantum Plasmas Beyond the Ring Approximation. Journal of Mathematical Physics, 7:616-626, April 1966. [18] T. L. Duvall, Jr. A dispersion law for solar oscillations. Nature, 300:242, November 1982. [19] T. L. Duvall, Jr. and J. W. Harvey. Rotational frequency splitting of solar oscillations. Nature, 310:19-22, July 1984. [20] W. Ebeling. Equation of state and Saha equation of partially ionized plasmas. Phys ica, 38:378-388, May 1968. [21] W. Ebeling. Coulomb interaction and ionization equilibrium in partially ionized plas mas. Physica, 43:293-306, July 1969. 173 [22] W. Ebeling. Statistical derivation of the mass-action law for interacting gases and plasmas. Physica, 73:573-584, May 1974. [23] W. Ebeling, W.D. Kraeft, and D. Kremp. Theory of bound states and ionization equilibrium in plasmas and solids. April 1979. [24] P. P. Eggleton, J. Faulkner, and B. P. Flannery. An Approximate Equation of State for Stellar Material. A. & A., 23:325, March 1973. [25] D. Gough. Inverting helioseismic data. Solar Phys., 100:65-99, October 1985. [26] D. 0. Gough. On the rotation of the sun. Royal Society of London Philosophical Transactions Series A, 313:27-38, November 1984. [27] D. 0. Gough. Asymptotic sound-speed inversions. In D. 0. Gough, editor, NATO ASIC Proc. 169: Seismology of the Sun and the Distant Stars, pages 125- 140, 1986. [28] D. 0. Gough. Linear adiabatic stellar pulsation. In J.-P. Zahn & J. Zinn-Justin, editor, Astrophysical Fluid Dynamics - Les Houches 1987, pages 399-560, 1993. [29] D. 0. Gough and M. E. Mcintyre. Inevitability of a magnetic field in the Sun's radiative interior. Nature, 394:755-757, August 1998. [30] H. C. Graboske, D. J. Harwood, and F. J. Rogers. Thermodynamic Properties of Nonideal Gases. I. Free-Energy Minimization Method. Physical Review, 186:210- 225, October 1969. [31] G. M. Harris. Attractive Two-Body Interactions in Partially Ionized Plasmas. Physi cal Review, 125:11311 140, February 1962. [32] G. M. Harris, J. E. Roberts, and J. G. Trulio. Equilibrium Properties of a Partially Ionized Plasma. Physical Review, 119:1832-1841, September 1960. [33] H. J. Hoffmann and W. Ebeling. On the equation of state of fully ionized quantum plasmas. Physica, 39:593-598, September 1968. [34] K. Huang. Statistical Mechanics, 2nd Edition. April 1987. 174 [35] D. G. Hummer and D. Mihalas. The equation of state for stellar envelopes. I- an occupation probability formalism for the truncation of internal partition functions. ApJ., 33 1:794--- 814, August 1988. [36] T. Kihara, Y. Midzuno, and T. Shizume. Virial Coefficients and Intermolecular Poten tial of Helium. Journal of the Physical Society of Japan, 10:249, April 1955. [37] D. J. Lacks, G. Goel, C. J. Bopp, IV, J. A. van Orman, C. E. Lesher, and C. C. Lundstrom. Isotope Fractionation by Thermal Diffusion in Silicate Melts. Physical Review Letters, 108(6):065901, February 2012. [38] A. I. Larkin. Thermodynamic functions of a low temperature plasma. Sov. Phys. JETP, 11:1363, October 1960. [39] A. Liang and W. Dappen. Emulating the Opal Equation of State. In D. Danesy, editor, SOHO 14 Helio- and Asteroseismology: Towards a Golden Future, volume 559 of ESA Special Publication, page 548, October 2004. [40] S. H. Lubow, E. J. Rhodes, Jr., and R. K. Ulrich. Five minute oscillations as a probe of the solar interior. In H. A. Hill & W. A. Dziembowski, editor, Nonradial and Nonlinear Stellar Pulsation, volume 125 of Lecture Notes in Physics, Berlin Springer Verlag, pages 300-306, 1980. [41] C. Ludwig. Diffusion zwischen ungleich erwarmten Orten gleich zusammengesetzter Losungen. Sitz. Math. Naturwiss. Classe Kaiserlichen Akad. Wiss., 20:539, August 1856. [42] G. Maki Harris. Equilibrium Properties of a Multicomponent Ionized Gas. J. Chern. Phys., 31:1211-1 220, November 1959. [43] D. Mao. IMPR OVED THERMODYNAMICS OF THE DENSE SOLAR PLASMA AND MOLECULAR -DYNAMICS SIMULATIONS OF THE NUCLEAR-REACTION RATES. May 2008. [44] D. Mihalas, W. Dappen, and D. G. Hummer. The equation of state for stellar envelopes. II - Algorithm and selected results. Apl., 331:815-825, August 1988. 175 [45] F. Miiller-Plathe. A simple nonequilibrium molecular dynamics method for calcu lating the thermal conductivity. Journal of Chemical Physics, 106:6082-6085, April 1997. [46] F. Miiller-P lathe. Reversing the perturbation in nonequilibrium molecular dynamics: An easy way to calculate the shear viscosity of fluids. Phys. Rev. E, 59:4894-4898, May 1999. [47] D. Reith and F. Miiller-Plathe. On the nature of thermal diffusion in binary Lennard Janes liquids. Journal of Chemical Physics, 112:2436-- 2443, February 2000. [ 48] F. J. Rogers. Phase Shifts of the Static Screened Coulomb Potential. Phys. Rev. A, 4:1145-1155, September 1971. [ 49] F. J. Rogers. Statistical mechanics of Coulomb gases of arbitrary charge. Phys. Rev. A, 10:2441-2456, December 1974. [50] F. J. Rogers. On the compensation of bound and scattering state contributions to the partition function. Physics Letters A, 61:358-360, June 1977. [51] F. J. Rogers. The energy balance and hydrodynamics of the solar chromosphere and corona. Proc. !AU Colloq., 36:3, March 1977. [52] F. J. Rogers. Formation of composites in equilibrium plasmas. Phys. Rev. A, 19:375- 388, January 1979. [53] F. J. Rogers. Integral-equation method for partially ionized plasmas. Phys. Rev. A, 29:868-879, February 1984. [54] F. J. Rogers. Occupation numbers for reacting plasmas - The role of the Planck Larkin partition function. Apl., 310:723-728, November 1986. [55] F. J. Rogers and H. E. Dewitt. Statistical Mechanics of Reacting Coulomb Gases. Phys. Rev. A, 8:1061- 1076, August 1973. [56] F. J. Rogers and C. A. Iglesias. Rosseland mean opacities for variable compositions. ApJ., 401:361-366, December 1992. 176 [57] F. J. Rogers and A. Nayfonov. Updated and Expanded OPAL Equation-of-State Tables: Implications for Helioseismology. Apl., 576:1064- 1074, September 2002. [58] F. J. Rogers, F. J. Swenson, and C. A. Iglesias. OPAL Equation-of-State Tables for Astrophysical Applications. Apl., 456:902, January 1996. [59] J. P. Severinghaus, T. Sowers, E. J. Brook, R. B. Alley, and M. L. Bender. Timing of abrupt climate change at the end of the Younger Dryas interval from thermally fractionated gases in polar ice. Nature, 391:141- 146, January 1998. [60] C. Soret. Concentrations differentes d'une dissolution dont deux parties sont a' des temperatures differentes. Arch. Sci. Phys, Nat., 2:48, August 1879. [61] E. A. Spiegel and J.-P. Zahn. The solar tachocline. A & A, 265:106-- 114, November 1992. [62] M. J. Thompson, J. Christensen-Dalsgaard, M. S. Miesch, and J. Toomre. The Inter nal Rotation of the Sun. Annual Review of Astronomy &Astrophysics, 41:599-643, 2003. [63] G. Uhlenbeck and E. Beth. The quantum theory of the non-ideal gas. I. deviations from the classical theory. Physica, 3:729, October 1936. [64] R. K. Ulrich. The influence of partial ionization and scattering states on the solar interior structure. Apl., 258:404- 413, July 1982. [65] R. K. Ulrich and E. J. Rhodes, Jr. The sensitivity of nonradial P mode eigenfrequen cies to solar envelope structure. Apl, 218:521-529, December 1977. 177

Abstract (if available)

Abstract The developments in helioseismology ensure a wealth of studies in solar physics. In par- ticular, with the high precision of the observations of helioseismology, a high-quality solar model is mandated, since even the tiny deviations between a model and the real Sun can be detected. ❧ One crucial ingredient of any solar model is the thermodynamics of hot-dense plasmas, in particular the equation of state. This has motivated efforts to develop sophisticated theoretical equations of state (EOS). It is important to realize that for the conditions of solar-interior plasmas, there are no terrestrial laboratory experiments

Linked assets

University of Southern California Dissertations and Theses

Conceptually similar

PDF

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

PDF

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

PDF

On the Feynman path into the sun

PDF

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

PDF

The physics of pulsed streamer discharge in high pressure air and applications to engine technologies

PDF

Advanced simulation models and concepts for positron and electron acceleration in plasma wakefield accelerators

PDF

Improvements in the equation of state for the partially ionized plasmas of the solar interior

PDF

Finite size effect and Friedel oscillations for a Friedel-Anderson impurity by FAIR method

PDF

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

PDF

Modeling graphene: magnetic, transport and optical properties

PDF

Transient pulsed plasma for pollution remediation and energy conversion applications

PDF

Particle-in-cell simulations of kinetic-scale turbulence in the solar wind

PDF

Modeling nanodevices: from semiconductor heterostructures to Josephson junction arrays

PDF

The cathode plasma simulation

PDF

Three-dimensional kinetic simulations of whistler turbulence in solar wind on parallel supercomputers

PDF

Light management in nanostructures: nanowire solar cells and optical epitaxial growth

PDF

Deterministic evolutionary game theory models for the nonlinear dynamics and control of chemotherapeutic resistance

PDF

An experimental investigation by optical methods of the physics and chemistry of transient plasma ignition

PDF

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

PDF

Hot electron driven photocatalysis and plasmonic sensing using metallic gratings

Asset Metadata

Creator Lin, Hsiao-Hsuan (author)

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

School College of Letters, Arts and Sciences

Degree Doctor of Philosophy

Degree Program Physics

Publication Date 05/03/2012

Defense Date 03/22/2012

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

Tag equation of state,helioseismology,OAI-PMH Harvest,plasma physics,solar physics,stellar evolution

Language English

Contributor Electronically uploaded by the author (provenance)

Advisor Daeppen, Werner (committee chair), Dappen, Werner (committee chair), Chang, Tu-nan (committee member), Haas, Stephan W. (committee member), Kunc, Joseph (committee member), Nakano, Aiichiro (committee member)

Creator Email hsiaohsl@usc.edu,sanabel18@gmail.com

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

Unique identifier UC11290206

Identifier usctheses-c3-27826 (legacy record id)

Legacy Identifier etd-LinHsiaoHs-734-1.pdf

Dmrecord 27826

Document Type Dissertation

Rights Lin, Hsiao-Hsuan

Type texts

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

Access Conditions The author retains rights to his/her dissertation, thesis or other graduate work according to U.S. copyright law. Electronic access is being provided by the USC Libraries in agreement with the a...

Repository Name University of Southern California Digital Library

Repository Location USC Digital Library, University of Southern California, University Park Campus MC 2810, 3434 South Grand Avenue, 2nd Floor, Los Angeles, California 90089-2810, USA