4:. {1.5.1 .1) .. n .u .5..,.y.n....‘... w) ‘ 511:4. .,.;;..,.2z. ..,. ..:.....,...a_ x. 31y _ L211}? 2 3,33. : as”? [.3 . . ii .L 1:12. m. .11? .3. i, :L, ._ t. is K L .a..:?:. z: . ‘ 631:3. , ‘ 1:2,: :7 1 .. .9. “its; V 2.4.1.1 This is to certify that the dissertation entitled INGREDIENTS FOR ACCURATE SIMULATIONS OF CONVECTION IN STELLAR ENVELOPES presented by Regner Trampedar‘h has been accepted towards fulfillment of the requirements for Eh. I) degree in Astrophysics difi/{W Major professor Date '2 Dec 20” ? MSU i: an Affirmative Action/Equal Opportunity [Institution 0-1277. PLACE IN RETURN BOX (0 r emo TO AVOID FINES retu MAY BE RECALLED with ve this checkout from your record. m on or before date due. earlier due date if requested. DATE DUE DATE DUE DATE DUE 6/01 c:/C|RC/DateDue.p65-p.15 INGREDIENTS FOR ACCURATE SIIV‘IULAII‘IONS OF CONVECTION IN STELLAR ENVELOPES b\ Regner Trampedach A DISSERTATION To be submitted to Michigan State University in partial fulfillment of the requirements for the degree of DOCTOR OF PHILOSOPHY Department of Physics and Astronomy 2004 ABSTRACT INGREDIENTS FOR ACCURATE SINIULENTIONS OF CONVECTION IN STELLAR ENVELOPES by Regner Trampedach I present the ingredients for high precision, 3D hydrodynamical sin‘iulations of convection in stellar atmospheres, as well as a number of applications. I have devel— oped a new scheme for evaluating radiative transfer, an improved equation of state and I have investigated a number of directions for improving the numerical stability of the convection simulations. The equation of state (EOS) used "for the simulations. is updated by including post-Holtsmark micro—field distrilmtions and relativistic electron—degeneracy as pre- viously published. I have further included quantum effects, higher-order Coulomb interactions and improved treatment of extended particles. These processes (except relativistic degeneracy) have a significant effect in the solar convection zone, and most of them peak at a depth of only 10 Mm. I also include a range of astrophysically sig— nificant molecules, besides Hg and the I’I'zi—ion. This EOS will be used directly in the convection simulations, providing the thermodynamic state of the plasma, and as a foundation for a. new calculation of opacities for stellar atmospheres and interiors. A new scheme for evaluating radiative transfer in dynamic and multi—dimensional stellar atmosphere calculations is developed. The idea being, that if carefully chosen, very few wavelengths can reproduce the full radiative transfer solution. This method is based on a calibration against a full solution to a. ID reference atmosphere, and is therefore not relevant for static 1D stellar atmosphere modeling. The first tests of the method are very promising, and reveal that the new method is an improven‘ient over the former opacity binning technique. The range of con\i'ective fluctuations is spanned more accurately and not only the radiative heating, but also the first three angular moments of the specific intensity, can be evaluated reliably. Work on implementing the method in the convection—code, is in progress. These developments will be employed in the future for a number of detailed simu— lations of primary targets for the upcoming, space-based. astero—seismology missions, and will include oz Cen A and B, nBoo, Procyon and ,13Hyi. \‘Vork on a 10 Mm deep solar simulation was severely hampered by numerical instabilities, but investigating the issue has revealed a number of potential solutions that will be tested in the near future. The work on individual stars will soon be superseded by an effort to compute a grid of convection simulations in Terr, logg and metallicity. [Fe/H]. in the spirit of present—day, grids of conventional atmosphere models. iii ACKNOWLEDGMENTS This dissertation would not have been possible without the help of a great num— ber of great people. I would in particular like to thank my advisor, Dr. Robert F. Stein, for generous hospitality during my visits to MSU before I started in the graduate program, for inviting me to pursue a PhD. with him at MSU. and for a great number of helpful insights during all these years. I would like to thank Dr. W'erner Dappen, at University of Southern California (USC) for generous hospitality during my numerous visits to Los Angeles, and for guiding me through the intricacies of the equation of state, as well as China Town and Old Town, LA. I would also like to thank Martin Asplund for hospitality during my visits in Upssala, Sweden, and, lately, Canberra, Australia Doing research entails a lot of reading-up on what other people have been doing on the subject you are working on; To learn about the subject—learn how to do, and how not to do it, find a new angle to a problem, etc. As you can see from a quick glance at the bibliography, quite a lot of reading has gone in to this research, and if I had, had to visit the library and make a photocopy of each of those papers, I would barely have made it through the second page of references. This dissertation and the included papers, have benefited tremendously from NASAs Astrophysics Data System (ADS) and I send the people who power the abstract server, my deep—felt gratitude. It has been delightful to share office with friend and long-time convection— colleague Dali Georgobianiv-and very convenient, both in terms of coffee—breaks and discussions about work. I would also like to thank Dr. Horace Smith for initiating the excellent institution of “bad, fifties sci—fi nights”, during which I got acquainted with blissfully obscure productions from an emerging Danish film industry. The mere mentioning of titles such as, “Reptilicus” and “Journey to the 7”] Planet” now evokes a strange mix of patriotism and, well—embarrassment. Last, but by no means least; This dissertation has been powered by oatmeal— raisin cookies, lovingly provided by my wife, Charlotte l\r’Iudar. Her patience, support and love has enabled me to keep my sanity through this endeavor. TABLE OF CONTENTS LIST OF TABLES ................................ LIST OF FIGURES ............................... 1 Introduction .................................. 2 Hydrodynamics ................................ 2.1 The Navier—Stokes Equations ....................... 2.2 Diffusion .................................. 2.2.1 Second Order Diffusion ...................... 2.2.2 Fourth Order Diffusion ...................... 2.2.3 Quenched Diffusion ........................ 2.2.4 A Gallery of Diffusion Coefficients ................ 2.2.4.1 The Original Mixed 2+4—Order Version ........ 2.2.4.2 The Original Form used with Quenched Diffusion . . 2.3 Alternative solutions ........................... 2.3.1 Subtracting a Smooth Average ................. 2.3.2 Monotonic Interpolation Schemes ................ 3 Input Physics ................................. 3.1 Equation of State ............................. 3.2 “A Synoptic Comparison of the MHD and the OPAL Equations of State” ................ 3.3 Pressure Ionization in the MHD Equation of State ........... 3.3.1 Introduction ............................ 3.3.2 Pressure Ionization in MHD ................... 3.3.3 A brief review of the MHD equation of state .......... 3.3.4 Examination of the influence of the \11 term ........... 3.3.5 Conclusion ............................. 3.4 MHD 2000 ................................. 3.4.1 Introduction ............................ 3.4.2 Chemical Compositions ...................... 3.4.3 Relativistic Electrons ....................... 3.4.4 The Free Energy of the Coulomb Interactions ......... 3.4.4.1 Quantum Diffraction .................. 3.4.4.2 The Free Energy of the Exchange Interaction vi 3.4.4.3 Coulomb Interactions Beyond Debye—Hiickel ..... 84 3.4.4.4 Shielding by Bound Electrons ............. 9 3.4.4.5 The New [74 ....................... 9' 3.4.5 Improved Micro-field Distribution ................ 9‘ 3.4.6 Including More Molecules .................... 9' 3.4.7 Energy Levels of the HQ— and Iii-molecules ........... 10 3.4.8 Conclusion ............................. 10? 3.5 Opacity .................................. 10: 4 Radiative Transfer .............................. 10! 4.1 Opacity sampling for 3D convection simulations ............ 10! 4.1.1 Introduction ............................ 10! 4.1.2 Primer on Opacity Binning ................... 11f 4.1.3 Primer on Opacity Sampling ................... Ill 4.1.4 SOS ................................ 11' 4.1.5 1D radiative transfer ....................... 121 4.1.6 Spanning convection ....................... 12f 4.1.7 Conclusion ............................. 12: 5 Improvements to stellar structure models, based on 3D convec- tion simulations ............................... l2t 5.1 T-T relations from convection simulations ................ 12? 5.1.1 Introduction ............................ 12! 5.1.2 The Basis for T-T relations .................... 13f 5.1.2.1 Grey Radiative Transfer ................ 133 5.1.2.2 Average Radiative Transfer .............. 133 5.1.3 The 1D-envelopes ......................... 13( 5.1.3.1 Implementation of the T—T relation .......... 132 5.1.4 The 3D—sin‘1ulations ........................ 13! 5.1.5 Averaging procedures ....................... 14: 5.1.6 The solar '[l—T relation ...................... 14" 5.1.7 Fitting formulas .......................... 14! 5.1.8 Results and Discussion ...................... 15f 5.1.9 The depth of outer convection zones .............. 161 5.1.10 Conclusion ............................. 16..l 5.2 Calibrating the mixing—length ...................... 162 5.2.1 Introduction ............................ 162 5.2.2 Mixing—length vs. 3D convection ................. 171 5.2.3 The simulations .......................... 17‘ 5.2.4 The envelope models ....................... 17' 5.2.5 Matching to envelopes ...................... 17! 5.2.6 Results and discussion ...................... 18‘ 5.2.6.1 Depth of the Solar convection zone .......... 18. 5.2.7 Conclusion ............................. 18' vii 6 Conclusions .................................. 190 6.1 Summary ................................. 190 6.2 Future Work ................................ 193 6.2.0.1 The convection simulations .............. 193 6.2.0.2 Applications ...................... 194 Bibliography ................................... 196 viii 9 10 LIST OF TABLES Chemical mixtures in [j\7,-/NH] (see text) ................ 6 Coefficients, a”, for Eq. (83). ...................... 8 Coefficients, 1),], for Eq. (84). ...................... 8 Coefficients, CU, for Eq. (85). ...................... 8 Coefficients for Jex, Eq. (103) ...................... 8 Coefficients for g(A). Eq. (106) ..................... 8' Parameters and derived convection zone depths for the seven simula- tions ..................................... 141 Coefficients for Eq. (171). ........................ 1.5? Response to changes in (11. lllh‘ and a ................... 16: Parameters for the simulations. and derived as and convection zone depths, dcz. ................................ 17‘ LIST OF FIGURES Schematic of a pure two—zone oscillation. The thin curve is an example of a smooth interpolating function, illustrating the zero derivatives at the grid—points, which are marked with X. ............... 11 Schematic of the response to a. step function (solid line with dia- monds). The dashed lined shows the corresponding linear second- order difference, and the solid line is the linear fourth—order difference. 13 Cubic—spline interpolation (dashed line) across a step—function (dia- monds). The solid line shows a cubic interpolation where the deriva— tives have been forced to be zero (not a spline). ............ 20 This figure shows contours of the difference in pressure between in— cluding and ignoring the \Il—term, in the sense loglolAP/PI. The dashed lines going from lower left to upper right are tracks of stellar structure for stars in the mass range 0.6-1.511/19, as indicated above each track. The dotted contour lines mark the hydrogen and second helium ionization zones respectively, with contour lines in steps of 10%. 59 The relative difference in mean—electron—mass, he, as a function of temperature in a solar stratification. The solid line is for the GN96 mixture truncated from 30 to 11 elements, and the dashed line is the optimized 11 element composition, as listed in Tab. 1 .......... 66 The two degeneracy factors, 0,, Eq. (64), and 9,, Eq. (78). The solid line shows the non—relativistic case, for which 68(7), 0) 2 (96(7), 0). The dashed lines display the relativistic 1),, for /3 up to 1.0 in steps of 0.1, and the clash—dotted lines show the same for Oe. The behavior in ,3 is monotonic in both cases ......................... 74 The quantum diffraction term, as evaluated from Eq. (65) and (66), compared with the Padé—formula given by Riemann et a1. (1995) and with analytical expansions for small and large 7 ............. 75 10 11 12 13 14- The RMS—deviations of the ? evaluated from the fifth-order binomial fits, from the exact expression, Eq. (81). The dotted contours are spaced 0.0125% apart and solid contours are 0.05% apart. ...... The exchange integral, Jex(7}, L3). as function of degeneracy parameter 7}. The upper solid line shows the numerical integration of [ex(l]), indistinguishable from the fit of Eq. (93). The dot—dashed line (right hand axis) is the relative deviation between the numerical integration and the fit for [ex. The other lines show the behavior of JPx with increasing L3 (in steps of 0.1) going from top to bottom. The fit of Eq. 103 is shown with dashed lines, but only discernible above 7} c: 30. The effect of higher order correlations on the free energy. The solid line is my fit to g 2 F/FDH. Eq. (105) and the dashed line is the corresponding gr; : If/UDH. The long-dashed line diverging at. -0.5. is the Abe cluster expansion (Abe 1959). and the dot—dashed line is the fit to the fluid—phase free-energy presented in Stringfcllow cl al. (1990). The vertical dashed line shows the location of the transition to the crystalline phase. ......................... Ratios of molecular number-densities between the present EOS and those obtained with the EOS in the Uppsala-package. which is part of the MARCS stellar atmosphere-models. The thin, horizontal. dashed line marks a ratio of one .......................... Heating q, flux H, J, B, 31f and Rosseland opacity for four of the test cases (horizontally). In the three upper panels. the dashed lines (right axis) are the differences in the sense (SOS — monochromatic). In the lower panel the dashed line is the Rosseland opacity from my method. The horizontal dotted lines are the zero—points for the differences, and the vertical ones show the location of T : 1. From left to right, the columns represent a red giant, a cool solar-type star, a metal-poor “Sun” and Procyon———a sub-giant. .................... Comparison of the heating from the opacity binning (bin, thin lines) and my new method (SOS, thick lines), with the “monochromatic” case (see text). A (. . .) denotes horizontal averages over a 1.5><6 Mm, 82X200 point vertical slice from a solar simulation (solid). RMS de— notes the corresponding horizontal root—mean—square averages (dashed lines). ................................... The effect of various averaging techniques applied to a. solar simula— tion (see text for details). Notice how cases a) and 1)) follow closely in the atmosphere down to (T) : Terr, at log10 7' 2 —0.2. ....... J.) k 9! 10 12 00 Comparison of the T—T relation from the simulation. with some often used solar T-T relations. a) The difference in Rosseland T-T relations between the simulation and an Atlas9 atmosphere model (Kurucz 1993; Castelli et at. 1997). b) Differences. on a mow—scale. between the simulation and an Atlas9 model, and four semi—empirical atn‘iospl‘ierc models: The model by Holweger & .\'Iiiller (1974). the classical fit. by Krishna Swamy (1966a: 19661)), the HSRA model (Cingcrich (i ai. 1971) and the VAL model (\I’crnezza el al. 1981). ........... 145 Plot of the differences between the actual T—T relations and the global fit (solid line) and the individual fits (dashed line). The horizontal solid and dashed lines are the corresponding RMS deviations from the fits. The shaded areas show the RMS of the temporal variation of the temperature ............................. 154 Plot comparing the global (triangles) and individually (filled circles with error—bars) fit parameters for the T-T relation-fits. The line segments show the loglOT gradient of the global fits. The error bars are the Rh’lS—scatter from T-T relation fits to single time-steps of the simulations. ................................ 156 The change of the T-T relation with stellar mass. on the zero-age main—sequence. The horizontal dotted line indicates the effective tem— perature, and the vertical dotted line, optical depth unity. The dashed U! line shows the grey atmosphere. ..................... 1 The change of the T-T relation with gravity, for fixed. solar Teff : 5777 I\'. The horizontal dotted line indicates the effective temperature. and the vertical dotted line, optical depth unity. The dashed line shows the grey atmosphere. ....................... 157 Scaled T-T relations for the seven simulations, relative to the solar simulation in the sense: T. X (Tome/Tea.) —— 7i?- The more vigorous the convection, the flatter the T-T relation ................ 165 The position of the seven simulations in the HR—diagram. The size of the symbols, reflects the diameter of the star. I have also plotted evolutionary tracks, with masses as indicated, to put the simulations in context .................................. 176 A plot of the 0’s found from the matching procedure (stars), com- pared with the linear fit from Eqs. (183) and (184) (triangles). The (lower) diamonds show the 2D calibration by Ludwig 61 al. (1999), which I have also multiplied by 1.14 to agree with my result for the Sun (upper diamonds). The line—segments show the local LogT- derivative of the fitting expressions .................... 181 xii This figure shows a’s behavior with Teff and gsurf. The surface is the fit given in Eq. (183). The dashed lines shows evolutionary tracks for stellar masses as indicated. The stars show the 0’s as listed in Tab. 10 and the little circles, connected by lines, show the projection onto the fitting-plane. ............................. 183 xiii Chapter 1 Introduction Understanding stars is a very multi—faceted endeavor. Stars are basically self—gravitating “Great Balls of Fire” (Lewis 1957) with their cores having temperatures and densities high enough for nuclear fusion to take place, and supply the star with energy. The stability of the star is provided by the pressure of the gas counter-acting the inward force of gravity, establishing a hydro-static equilibrium. The equation of state (EOS), supplying the thermodynamical relations between temperature, density and pressure, must be a key ingredient in our model of a star. The energy produced in the core travels out through the star, in one of three ways; either by conduction, radiation or convection. Simple conduction is only important for very dense objects, where the gas begins to act as a metal. This happens for, e.g., neutron stars and white dwarfs, but also in the centers of gaseous planets. Radiative transfer of energy takes place in all stars and is also a key ingredient for any stellar model. The opaqueness, or opacity, of the gas determines the ease with which light can travel through the star. With a higher opacity, the energy cannot escape as freely and the interior will heat up, until the energy that escapes matches the energy that is produced inside, setting up a. radiative equilibrium. The opacity is composed of a number of interactions between light and matter. Cross—sections for the various processes are evaluated from detailed knowledge of the quantrim—mechanical wave- functions of all possible states of the particles involved. Calculation of these cross- sections is a field of research in itself, and the vast calculations involved have only been possible within the last two decades with the progress of computing-power. To calculate the opacity, relevant for a stellar model, we need to multiply the cross—sections of individual processes by the number-densities of the interacting par- ticles, and add-up all the processes. The number—densities of the various particles, and the population of the energy—levels of each kind of particle, are also results of the equation of state, which therefore has a dual rdle. The equation of state is also a major atomic physics project, as the interactions between the particles constituting the gas, have to be accounted for. Although the gas is electrically neutral on a whole, most of the particles will be charged, being ionized due to the high temperatures. The particles therefore feel the presence each others by their charges, especially at higher densities where close encounters will be more frequent. In close encounters quantum effects will have an effect. The particles will no longer be point—like, but “smeared out” according to Heisenberg’s uncertainty principle (Heisenberg 1927), and encounters between identical particles will be governed by Pauli’s exclusion principle (Pauli 1925). All of these effects (and many more) determine the thermodynamics of the gas, as well as the ionization- balances and the population of energy—levels, needed for opacity calculations. The last mode of energy-transport is convection, which consists of macroscopic bubbles or flows of warmer than average gas. moving out towards lower temperatures. shedding their excess heat and sinking down to higher temperatures to be recycled. If the opacity is very high and the energy has a hard time escaping via radiation, the gas will become unstable towards convection and this will take over the transport of energy. Convection is inherently a dynamic process. involving a lot of non—linear phe- nomena. As such, it is not amenable to simplified theories, and a good and relatively simple description of convection suitable for stellar modeling has eluded us. The so-called mixing-length-formulation of convection is simple and very suitable for use in stellar structure codes. It is, however, fraught by it’s assumptions breaking down where a theory of convection is most needed, and by free parameters that have to be calibrated from outside the theory. For envelope convection, i. (2., where the convection zone reaches all the way to the surface, but does not necessarily reach the center, further complications arise at the surface. Here radiative transport begins to become important again, as the distance to the surface (where light can escape) becomes comparable with the mean-free-path of photons. This region is called the stellar atmosphere, the modeling of which forms its own field of research. As the gas becomes optically thin, and photons can begin to escape freely, evaluation of radiative transfer becomes much more complicated. The radiative transport of energy can no longer be described by a single opacity, but has to be evaluated at upwards of 104 wavelengths. Combine this transition between optically thick and optically thin, with the transition between convective and radiative transport of energy, and you have a very i1‘1tractable problem. The best way of studying this region, is through three-dimensional (3D) hydro dynamical simulations, for which a very small number of approximations have to b< introduced. They are based on the same equation of state and opacities also user in conventional (e.g., 1D) stellar models, but the hydrodynamics of the problem i: solved explicitly, based on the fundamental laws of conservation of mass, momentun and energy. Such simulations can be used for the same purposes as conventiona stellar atmosphere models, e.g., interpretation of observations including abundancc analysis, and as upper boundary conditions for stellar structure and evolution models The above sketch of stellar structure hopefully gives an impression of the inter- relatedness of all the phenomena participating in making a star a star, and serves a: the motivation behind the present dissertation. In chapter 2, the hydrodynamical convection simulations are introduced, and z stability problem is discussed. In chapter 3 I discuss the atomic physics entering botl the convection simulations and the stellar structure models, presenting the curren‘ state of astrophysical equations of state in Sect. 3.2 and 3.3, and proposing some improvements to one particular EOS-project, in Sect. 3.4. In Sect. 3.5 I give a short overview of recent developments in opacity calculations, and expectations for the near future. Chapter 4, deals with the evaluation of radiative transfer in the convection sim- ulations, and how to bring the problem down to tractable dimensions, and yet only lose a little in accuracy compared to solutions including 2; 105 wavelength-points. Chapter 5 presents a few applications of convection simulations in conjunctior with stellar structure calculations. In Sect. 5.1. the convection simulations are used as boundary conditions for stellar structure models, which are used in Sect. 5.2 for a calibration of the main parameter, a, of the mixing—length formulation. A summary and an outlook for the future. is offered in chapter 6. Chapter 2 Hydrodynamics The code for the convection simulations used in this work, was written by Nordlund & Stein and is further described in Nordlund (1982), Stein (1989). Nordlund 8: Stein (1990) and Stein & Nordlund (2003). The simulations are of compressible convection, subjected to realistic radiative transfer including the effect of line—blanketing. The boundaries are transmitting, with the entropy of the in—flow being adjusted, so as to result in the desired flux for the given simulation. A state of the art equation of state, presented in Sect. 3.2, is used, together with equally high quality opacities. The properties and morphology of the convection in the simulations, is described by, 6.9., Stein & Nordlund (1989) and Nordlund & Stein (1991). 2.1 The Navier-Stokes Equations The conservative or divergence form of the Navier—Stokes equations is a 7% + V-(Q’u) =0 (1 091‘ + \7 (ouu—a) :fe (.3 0t aLoftot ,. at + V‘fQCtotU—O'U-l—(I):f,mu (3 where Q is the density, elm. = 5 + é-u'z is the total specific internal energy, It is tlr velocity field, q is the heat-flux vector and fe is the external force per unit volume. The stress—tensor, 0’, can be written az—PgI+T, (4 where Pg is the gas pressure, I is the unity tensor with components I,-,- = 5,,- and C is the viscous stress tensor with components Tij- Expanding Eqs. (l)—(3), eliminating occurrences of Eq. (1) and substituting Eq (4), we arrive at the following three equations: the equation of mass conservation often referred to as the continuity equation, 8? -:—-V—V-. 5 a, uee’w ( xI the equation of momentum conservation, also called the equation of motion, a Qfiz—pu-Vu+f0—Vpg+VT. (. and the equation of energy conservation (with E = 95). 0E_ —5t—_—\7-(uE)+(T—P)V-u+V-q. ( Since the densities and pressures change by 4—5 orders of magnitude from the t( to the bottom of the simulations, it is advantageous to precondition the equatior for use in the highly stratified case of stellar atmospheres. This is done by rewritii (5) — (7) in terms of logarithmic density, lng, velocity, u, and internal energy per un mass, 5: (‘31 81:9 : —u - Vlng — V - u. (i 8 P 1 —u:—u-Vu+g———ngnP+—VT (1 at e e and 8 T—P 1 ,—€—=—u~Ve+ v-u+—v.q, (II ()t p Q where, for the external force, we use fe/Q :g, the gravitational acceleration, whi< is assumed to be constant with depth. These equations are much more well behave than Eqs. (5) — (7) and the various derivatives are also smoother. As the gradient of a flux is the time—derivative of the quantity transported l: the flux, we define de E 9—1V - q, where Qrad is the specific (radiative) heating. O O 2.2 Diffusmn In the energy equation the term in the viscous stress—tt—nisor. T is the dissil')ati« term, which we will re-write in a manner similar to Qrad: (2,45,. E p'ltl'V - u. The T-term in the momentum equation, is the (velocity-) diffusion. We w re-write this in component form, as 0v,- _ l ('3 (3&1, _ 1 0'1},- 9-z _ ‘20 ”1'99 ‘_Z‘) S (1 ( diff Q j .1] (1] Q j (.lj (1 Bu,- Comparing Eq. (9) and (12) we see that the diffusion coefficients, V]. have dimensio [L2/T]. With the dissipation tensor, T, so defined. the dissipation is T l (311., (711,- 011, 0a,- V'QCZ—V'UZ“ Tit—Z V'.—.-———- V' .—— 1 Q I” p p ,2]: J (3:17]- :1: J (3r, (3.1“,- {21: J (Jr,- ( 2.2.1 Second Order Diffusion Similar to Eq. (11), the second-order diffusion for a per—unit-mass quantity f, ad a contribution of the form f)f_ 1 (2 . 0 at T + anrj lug-”90:13 The summation is over directions, j = (1?, y and If the order of the numerical differencing is higher than 1, the derivatives wil be oblivious to two-zone oscillations, as depicted in Fig. 1. It can easily be show that the derivative of a smooth interpolating function will be zero at the peaks of: two—zone oscillation. For higher than first—order differencing schemes, it is therefore necessary to ex pand Eq. (14) into 8f Off 0122,- alng 33f .. E—+;[V2Jal—j+ (001] +V2J (21] UJIJ' . (I) We are using a first-order differencing scheme. but the expansion will nonetheles: prove useful in the subsequent discussion. The linear finite difference version of the second derivative is .A2 . .. . . (Afffmfhn=Avsuvv+n—Anrsu=uw4w¢nv+uwar(w If the quantity is a pure oscillatory signal with amplitude, d, and a wavelengtl of two grid-points (a so-called two-zone oscillation), as shown in Fig. 1, then we haw fi—i = fi+1 = —f,' = id, and consequently Nt:_M; 07 A17,- Ar, 7 for a uniform grid. Adding to f, a term proportional to Affi/Arcj, will therefort dampen a two~zone oscillation. Similar to the derivation of Eq. (17), we find that the first—order first-derivatiw 10 Fig. 1. Schematic of a pure two—zone oscillation. The thin curve is an example of a smooth interpolating function, illustrating the zero derivatives at the grid—points, which are marked with x. is twice the function, phase—shifted by half a grid-point. The second—derivative con- tribution to the second—order diffusion will therefore in general dominate, over the first-order derivative term. 2.2.2 Fourth Order Diffusion In order to minimize unduly damping of physical features, we can go to higher order schemes, which will have a higher ratio of two—zone—oscillations to smoother features. Fourth—order diffusion acids a contribution of the form 8/ 1 a2 021* the sign of which will be justified below. 11 The corresponding dissipation is . '2 Q (3211,; (19 .4 I 2: V4 T_:— . use ,- . J 032 2,] J The linear, finite differences version of the third— and fourth—derivatives on uniform grid, are Aye E) Z A3ff'i — %) = A2f(i) _ A'Zf(,_1) (20 (A103 [43f] (,- _1 :-4c—v+va—n—wm+an> (a and 4A4f. 4- 3~1 3-1 - mo Afl oz _—_ f(z — 2) — 4f(z — 1) + 6f(z) — 4f(i + 1) + f(2. + 2). (23 In the case of a pure two—zone oscillation, we therefore have that A3f, = 8 f, am A4f, : 16f,, or in general AN], : 2Nf,. This is true for the linear difference scheme and for third—order differences it can be shown that ANf, : 12N/2f,-. To destroy a two-zone oscillation with the pure fourth-order term, we therefor need to subtract a term proportional to AU}. Using a linear difference scheme, two zone oscillation will therefore be four times more prominent in fourth—order deriva tives than in second—order derivatives. But what about physical features like shar] gradients? 12 I ' r ' T ' l 0 2 4 6 8 10 Fig. 2. Schematic of the response to a step function (solid line with diamonds). The dashed lined shows the corresponding linear second-order difference. and the solid line is the linear fourth—order difference. Fig. 2 shows the response of second- and fourth—order derivatives to a step function — the ultimate in large gradients. From the figure, we see how second- order diffusion simply smooths out the large gradient, spreading it over more zones. Fourth—order diffusion does the same, and more efficiently, but it also introduces some extra peaks on either side of the gradient. These extra peaks, coupled with the (global) cubic-spline interpolations and differentiations used elsewhere in the simula— tions, quickly and efficiently excite two—zone oscillations. This phenomenon is called ringing, and it is excited by sharp edges, as well as localized peaks. We want to keep the gradients as large as possible, and capture as much detail as possible, with a finite resolution of the simulations. We therefore need a scheme with a much larger response to two—zone-oscillations than to step—functions. For linear second—order derivatives, this ratio is 4 / 1, and for fourth-order derivatives it is 13 16/3253. So we do gain a little going from second— to fourth-order derivatives, bu at the expense of extra excitation of two—zone oscillations. 2.2.3 Quenched Diffusion One way around the problem of two-zone oscillations, is the so—called quenched diffu sion, introduced by Stein & Nordlund in the early ‘903 (private communications). Th quenching consists of using an expression for 1/4‘, that involves a (signed) first-orde difference, sealed with the ratio of (unsigned) third—order differences to (unsigned first-order differences, as _ A; MI Ar,- _ A1“,- max3 IAfI 7 (~24 centered on i — 3' The maxg—operator takes the maximum value in a three-poin neighborhood of i, that is, i — 1, i, i + 1. The third—order difference is defined by Ec (20). By scaling the first—order difference in this way one obtains a diffusive flux tha has the sign of the first—order difference, but the order of magnitude of the third~orde difference. In smooth parts of the solution the diffusive flux will thus be quencher (hence the name), and the effect of the diffusion operator will be of similar magnitud as that of a fourth—order diffusion operator. In steep parts of the solution, on th- other hand, the ratio of the third- to first-order difference will be of the order unit; or larger, and the diffusion will be similar to normal second-order diffusion. The diffusion contributions are combined into where the first part is a seconcl—order diffusion that is active mainly in shocks. Th second part is qualitatively similar to a fourth-order diffusion, but has the advantag that sharp edges and localized peaks do not give rise to the ringing that tends to de velop with a. normal fourth-order operator. In the combined second- and fourtl'l-orde diffusion, enough second—order diffusion must be present to counteract the tendenc; for ringing from the fourth—order operator. Since this tendency is not present whei (25) is used, one can remove the second—order diffusion altogether, except in shocks where significant local diffusion is always needed. With the quenched diffusion, only the shock capturing a2 term is kept in thl pure second-order part of the diffusion, and all three terms: advection-, shock— anr sound~wave-terms, are included in the quenched part (see Sect. 2.2.4.2). The dissipation that corresponds to (25) is 811,,- Bu,- Qvisc : Z 8'— VQ‘jf + V4jA3U2‘ - (26 23.7 .1] 1,] The diffusion at the top and the bottom are set to zero in both diffusion scheme to avoid boundary effects. All of the non-local operators take advantage of th: periodic properties in the horizontal directions and are skewed at the top and botton boundary. The quenched diffusion works very well, as long as the two—zone oscillations ii the first derivative are larger than the broader features in the first derivative — tha‘ is, when the two-zone oscillations are centered around zero. Otherwise it fails to pick up the two-zone oscillations. The latter is often the case deeper in the simulations 15 and was of less concern for the shallow simulations pursued so far, but. during worl on a. 10 Mm deep solar simulation, this proved detrimental. These symptoms hat also been noticed higher in the atmosphere of a number of different simulations introducing noise in the radiative transfer calculations and two-zone oscillations i1 temperature. 2.2.4 A Gallery of Diffusion Coefficients The various forms of the diffusion coefficients all contain combinations of three term with the following motivation: The (L1 -term is proportional to the fluid velocity to prevent ringing at sharp change: in advected quantities. a2 -term is proportional to a finite difference velocity convergence (when positive) and is necessary to prevent excessive steepening of shocks. a3 -term is proportional to the sound speed, and is necessary to stabilize sounc waves and weak shocks. Also needed to stabilize 2—zone oscillations. 2.2.4.1 The Original Mixed 2+4-Order Version /_\.:rr _. 122,- : (allujl + (igAj'uj)—Ei , (27‘ where the velocity convergence is defined as k/2 Ail-(My) = X [15,241 — 1th} > 0 . (28 iZI—k/2 16 If the fluid in two adjacent volume-elements are approaching, Arc-(ivy) is equal to tilt relative velocity, in k volume-elemei‘its around the two. The fourth-order diffusion coefficient is I/4j:(~/1(L1(ttjl+ (I3CS)# where 05 is the adiabatic sound-speed. 2.2.4.2 The Original Form used with Quenched Diffusion We use a diffusion suppression factor to suppress unnecessary shock—diffusion in the deeper layers where there are no shocl The second—order diffusion coefficient is and the fourth—order diffusion coefficient is As:- 1/4J- : —1—2]-[a1|uj|+ a-ng AIM”) + (L3CS:| . (3 2.3 Alternative solutions Neither of the two diffusion-schemes outlined above, the second plus fourth orde diffusion, or the quenched diffusion, are optimal for the convection simulations. Two zone oscillations are excited by the cubic—spline interpolation el‘nployed, and these twl diffusion schemes do not dampen them efficiently without also smoothing physica features. We need another scheme for either dampening or avoiding these two—zonn oscillations all together. 2.3.1 Subtracting a Smooth Average One straight forward solution to the above problem. is simply to subtract out tlr physical variation, to leave the undesired two-zone oscillations exposed: i + N f.’ = /1 — Z 10'ij . (33 j: Z- — TV where the weights 10,- define the functional form of the averaging procedure an< 2N +1 is the number of points involved in, the averaging. We have chosen a Gaussial averaging kernel , apt-(170)?) ZEif-‘iw expt—(j/OV) ’ (34. w, I where o is the width of the Gaussian and the denominator takes care of the nor malization. Using a full—width—half—maximum of W : 30\/ln2 of three points, ha: proven optimal for getting a large signal from the two-zone oscillations. The facto; of three in the expression for 14/ stems from the counting of grid points. Normall: 18 there would be factor of two to result in the width and not just the position. it, c half-maximum. In the discrete case. however. we also have to count the point in th middle. In choosing N there is obviously a trade—off between the computational expens and the “smoothness” of the average. Rather surprisingly a three point average i very close to a five point average, and a seven point average is indistinguishable fron the five point average. These are the results for real case scenarios, simply goini through a snapshot of a simulation in both horizontal and vertical directions, an< comparing f’ for 2N + 1 :3, 5 and 7. l 1 Usin W" = 3 and A" = I therefore results in the wei )‘hts u‘ :— , .1. ‘1‘ . and f) 1 2 4 _1, z 41.-. +1- 1 <35 which has first-order differences (36 From Eq. (20), we recognize this as the third—order difference of f. The procedurc is commutative, so whether subtracting the Gaussian average before or after thi differencing, gives the same result. By subtracting out a Gaussian average, we merel; change the first—order differences into third—order differences, immediately bringim us back to the problems of Sect. 2.2.2. The extra peaks introduced around stee] 19 gradients (see Fig. 2), do not depend on the choice of N. This idea has therefore been abandoned. 2.3.2 Monotonic Interpolation Schemes In Fig. 3 we show how a cubic—spline interpolation across a step-function, generates oscillations in the derivatives on either sides of the step. A cubic spline, is a se— ries of piece—wise cubic polynomials, connected by the requirement. that the second derivatives be continuous (de Boor 1978). We see from Fig. 3 that the interpolating function is not monotonic between grid—points, and that this is the reason for the excitation of two—zone oscillations around sharp edges. 1.2 1.0- 0.8— 0.6- 0.4- Fig. 3. Cubic-spline interpolation (dashed line) across a step—function (diamonds). The solid line shows a cubic interpolation where the derivatives have been forced to be zero (not a spline). There are a number of schemes for calculating monotonic interpolating functions, 20 of which some are local and compact and others are global in nature. The princi] of the requirement of monotonicity, is illustrated by the solid curve in Fig. 3. This a piece-wise cubic interpolation where the derivatives on either side of the step ha been set to zero. This results in a more visually pleasing behavior. This is admittec a rather subjective criterion for the quality of an interpolant, but a valuable first te The deciding test is how the scheme fares when applied to advection—, shock—tut. and wave—tests with analytical solutions. There are a great number of monotonic interpolation schemes on the marl (e.g., Fritsch & Butland 1984, and references therein). All of these schemes WOL‘ qualitatively result in the same interpolating function in the case shown in Fig. 3, b behave differently when presented with triangular functions, Gaussians, parabol: multiple step functions and similarly challenging tests. One monotonic scheme, which is rather pleasing for its simplicity, uses piece-w cubic—splines, disconnecting two parts of the spline where the solution is not mor tonic between grid—points. This obviously raises the question about the bounda conditions on the two parts of the spline. The solid curve in Fig. 3 is an example this scheme. The total variation diminishing (TVD) schemes, as introduced by Harten (198: use a measure of the overall amount of oscillation of a quantity v. N TVW‘U] = E; 11-61410) * with ~. (3 and require this quantity not to increase with time, i. A good presentation. of TV 21 was recently written by Trac & Pen (2003). In general, the order of the interpolation goes down across a discontinuity, but this is a small price to pay for increased stability and a more physical solution in the smoother parts of the simulation. There are several promising ways of solving the problems of instabilities in the convection simulations, and future investigations of these will reveal which is the best solution to the problem at hand. In parallel with this work, a new version of the code, specifically optimized for massive parallelization, is being written by Nordlund, Stein, Carlsson & Hansteen (private communications). This new version employs a staggered grid, with energy and density defined at cell—centers and velocities on cell-boundaries, which will most likely be stable against two-zone oscillations. Other problems are introduced with this formulation, though, in particular regarding the stability of the upper and lower boundaries of the simulation domain. 22 Chapter 3 Input Physics The atomic physics underlying a model of an astrophysical phenomena is critical for the quality of that model. In Sect. 3.2 I present an analysis and comparison between two of the leading equation of state projects; the so—called Mihalas—Hummer—Dappen (MHD) equation of state (Hummer & Mihalas 1988; Mihalas ei (Ii. 1988; Dappen et a1. 1988) and the OPAL equation of state pursued at Lawrence Livermore National Laboratory (Rogers 1986; Iglesias & Rogers 1995; Rogers et al. 1996, and references therein). An artificial term, \11, in the MHD EOS that ensures pressure ionization in cold and dense plasmas, has been suspected of contaminatii'ig this EOS under solar con— ditions, and is examined in Sect. 3.3. The close scrutiny of the MHD EOS in Sects. 3.2 and 3.3 paves the way for the improvements, aimed at broadening the range of ap}_)licability and increasing the accuracy of the MHD EOS, presented in Sect. 3.4. This work is equally aimed at the convection simulations and the general stellar structure and evolution problem. All 23 of the improvements are important for the latter, but for the convection simulations the addition of many more molecules will have the largest effect. The improved MHD EOS, presented in Sect. 3.4 will in the future be used as basis for a new opacity calculation, and in Sect. 3.5 I give a short. review of recent progress in the calculation of absorption coefficients. 3.1 Equation of State The equation of state is essential in most astrophysical contexts, both for the ther- modynamic relations it provides, and as the foundation for opacity calculations. Exploring the differences between the main EOS projects, I found that the most contested part of the solar EOS, the most uncertain and the most extreme plasma conditions, occur very close to the surface~less than 10 Mm below the photosphere. The deep 10 Mm solar convection simulation will get the full impact of the improve- ments presented in this chapter, but even the shallow 3 Mm solar simulations will be affected to an extent that can be measured by helioseismology (Christensen-Dalsgaard et al. 1988; Basu & Christensen—Dalsgaard 1997; Di Mauro et (Li. 2002). 24 3.2 “A Synoptic Comparison of the MHD and the OPAL Equations of State” by Trampedach, R., Dappen, W. & Baturin, V A., 2004, ApJ, (accepted) The currently most widely used equations of state in the astrophysical com- munity are the OPAL equation of state, pursued at Lawrence Livermore National Laboratory (Rogers 1986; Rogers et (Li. 1996), and the MHD equation of state (Hum— mer 85 Mihalas 1988; Mihalas 61‘. (1.1. 1988; Dappen ci (1i. 1988), which is part of the international Opacity Project (OP) (Seaton 1995', Berrington 1997). The present paper (Trampedach ct (ti. 2004c) compares these two equations of state to set the stage for further developing the MHD EOS in Sect. 3.4. 25 A synoptic comparison of the MHD and the OPAL equations of state R. Tra m pedachl Department of Physics and Astronomy. .A/fic/zi_qa72 Slate (11111111111111, [Casi Lansing, Ml 48824, USA trampedachOpa.msu.edu W. Dappen1 Department of Physics and .1’Ist7‘onomy, NSC, Los Angeles, CA 90089-1342, USA dappenQusc.edu and V. A. Batiurin Sle'l‘nberg Astronomical Ins/flute, Universitelskg Prosper-1‘. [3, Moscow 119899, Russia vabOsai.msu.su ABSTRACT A detailed comparison is carried out between two popular equations of state (FOS), the MihaIas-I-Iummer-Dappen (MHD) and the OPAL equations ofstate, which have found widespread use in solar and stellar modeling during the past two decades. They are parts oftwo independent efforts to recalculate stellar opacities: the international Opacity Project (OP) and the Livermore— basd OPAL project. We examine the difference between the two equations of state in a broad sense, over the whole applicable ,9 — T range, and for three different chemical mixtures. Such a global comparison highlights both their differences and their similarities. We. find that omitting a questionable hard—sphere correction, 7'. to the Coulomb interaction, greatly improves the agreement. between the MHD and OPAL EOS. We also find signs ofdiffer- ences that could stem from quantum effects not yet included in the MHD EOS, and differences in the ionization zones that are probably caused by differences in the micro-field distributions employed. Our analysis does not only give a clearer perception oft'he limitationsofeach equation ofstate for astrophysical applications, but also serves as guidance for future work on the. physical issues behind the. differences. The outcome should be an improvement ofboth equations ofstate. Subject headings: Atomic processes—lilqimtion of state—Plasmas—Sun: interior 1. Introduction Stellar modeling, and in particular helio- and asteroseismoiogy, require an equation of state and corresponding thermodynamic quantities that are smooth, consistent, valid over a large range oft'em- lTeoretisk Astrofysik Center, Danmarks Grundforsk— mngsfond, Institut for Fysik 0g Astronomi, Aarhus Uni- versrtet, DK-8000 Aarhus C, Den mark 26 peratures and densities, and that incorporate the most important chemical elements ofastrophysical relevance Christensen-Dalsgaard & Dappen (for a review see 1992). In astrophysics, the equation ofstate plays two basic roles. On the. one hand, it supplies the ther- modynamics neccessary for describing gaseous ob- jects such as stars and gas-planets. On the other hand it also provides the foundation for opac- ity calculations, in the form of ionization equi— libria and level populations. Thanks to helioseis— mology. the Sun has broadened this perspective. The remarkable precision by which we have now peered into the Sun, puts strong demands on any physics going into a solar model. This, to such a degree, that we can turn around the argument and use the Sun as an astrophysical laboratory to study Coulomb systems under conditions not yet achieved on Earth. Although the solar plasma is only slightly non- ideal, the tight observational constraints prompts the use of methods normally reserved for studies of more stronglymoupled plasmas. in this way the solar experiment addresses a much broader range of plasmas, e.g., Jovian planets, brown dwarfs and low-mass stars, as well as white dwarfs ((fauble et al. 1998). The two equation—of-state efforts we compare in this paper are associated with the two leading opacity calculations of the eighties and nineties. The MHD EOS (Hummer & Mihalas 1988; Mi— halas et al. 1988; Dappen et al. 1988) was devel— oped for the international Opacity Project (OP) described and summarized in the two volumes by Seaton (1995) and Berrington (1997). and the OPAL EOS is the equation ofstate underlying the OPAL opacity project at Livermore (Rogers 1986: Rogers et' al. 1996, and references therein). Another highly successful EOS address the ex— treme conditions in low-mass stars and giant plan- ets and include the transition to the fluid phase (Saumon Xx. Chabrier 1989; Saumon et al. 1995). la a trade-off between accuracy and range of va— lidity, this EOS has so far only been computed for H/He-mixtures, rendering it less suitable for he— lioseismic investigations. Comparisons with this EOS should, however, be an essential part of cf- forts to further develop precise stellar EOS. The OP and OPAL projects are based on two rather different philosophies; the chemical picture and the physical picture, respectively, as detailed in Sect. 2.1.] and 2.1.2. The effect of Coulomb interactions is reviewed in Sect. 2.2, and a cor— rection, 1', to these, that seems to account for a substantial part of the differences between the two formalisms, is explored in Sect. 2.2.2. Detailed comparisons between the MHD and OPAL EOS have proved very useful for discov- 27 ering the importance and consequences of several physical effects (Dappen et al. 1990: Dapper] 1992, 199(5). ln Sect. 3. we extend these comparisons to a systematic search in the entire 7‘—/) plane. and in Sect. 4 we take a closer look at the EOS under solar circumstances. The consensus of the last few years has been that in helioseismic comparisons the OPAL EOS is closer to the Sun than the MHD EOS (Christensen— l)alsgaard et al. l996) although both are remark- ably better than earlier theories. However. recent helioseismic inversions for the adiabatic exponent 7] : ((9liip/(9liig)ad (Basu et al. l999: Di Mauro et al. 2002) shows that the M Hf) EOS fares better than OPAL in the upper 7% of the sun includ- ing the ionization zones of hydrogen and helium. These new developments again highlights the im— portance of competing equation-of-state efforts and systematic comparisons such as the present. 2. Beyond Ideal Plasmas The simplest model of a plasma is a non- ionizing mixture of nuclei and electrons. obeying the classical perfect gas law. However. an ideal gas can be more general than a perfect gas. ldeal only refers to the interactions between particles in the gas. The interactions in any gas redistribute energy and momentum between the particles, giv- ing rise to statistical equilibrium. In an ideal gas these interactions do not contribute to the energy of the gas, implying that they are point interac— tions. Since the Coulomb potential is not a (5— function. real plasmas cannot be ideal. Deviations from the perfect gas law, such as ion- izat ion. internal degrees of freedom (excited st ates, spin). radiation and Fermi—Dirac statistics of elec— trons are all in the ideal regime. And the parti- cles forming the gas can be classical or quantum, material or photonic; as long as their interactions have no range. the gas is still ideal. All such ideal effects can be calculated as exactly as desired. The ideal picture, is however, not adequate even for the solar case. At the solar center, an ideal- gas calculation leaves about. 20% of the gas un— ionized. On the other hand, the mere size of the neutral (unperturbed) atoms, do not permit more than 7% of the hydrogen to be unionized at these densities, provided the atoms stay in the ground state and are closely packed. At. the temperature at the center of the Sun neither of these assump- tions can possibly hold and the mere introduction of size and packed immediately imply interactions between the constituents of the plasma and it is therefore no longer ideal. In a plasma of charges, Z, with average inter- particle distance (r), we define the coupling pa- rameter, F, as the ratio of average potential bind— ing energy over mean kinetic energy hBT deg P _ lmT(r) I (I) Plasmas with r >> 1 are strongly coupled. eg. the interior of white dwarfs, where coupling can be- come so strong as to force crystallization. Those with P << 1 are weakly coupled. as in stars more massive than slightly sub-solar. As one can suspect, P is the dimensionless cou— pling parameter according to which one can clas- sify theories. VVeakly-coupled plasmas lend to sys- tematic perturbat'ive ideas (eg. in powers of P). strongly coupled plasma need more creative treat— ments. Improvements in the equation ofstate be- yond the model of a mixture of ideal gases are dif- ficult, both for conceptual and technical reasons. 2.1. Chemical and Physical Picture The present comparison is not merely between two EOS-projects, but also between two funda- mentally different. approaches to the problem. The chemical picture is named for its foundation in the notion of a chemical equilibrium between a set of pre—defined molecules, atoms and ions. ln the physical picture only the “elementary" particles ofthe problem are assumed from the out— set — that is, nuclei and electrons. Composite. particles appear from the formulation. 2.1.]. Chemical Picture: MHD EOS Most realistic equations of state that have ap— peared in the last. 30 years belong to the chemical picture and are based on the free—energy minimiza- tion method. This method uses approximate sta- tistical mechanical models (for example the non— relativistic electron gas, Debye-Hiickel theory for ionic species, hard—core atoms to simulate pressure ionization via configurational terms, quantum me- chanical models ofatoms in perturbed fields, etc.) From these models a macroscopic free energy is 28 constructed as a function of temperature T. vol ume V, and the particle numbers N] ..... Nm 0 the in components of the plasma. At given T am i". this free energy is minimized subject to tlu st'oichiomet‘ric constraints. The solution of thi: minimum problem then gives both the equilibrinn concentrations and. if inserted in the free energ] and its derivatives. the equation of state and the thermodynamic quantities. Obviously. this procedure automatically guar antees thermodynamic consistency. As an exam ple. when the ("oulomb pressure correction (set Sect. 2.2) to the ideal-gas contribution originate: from the free energy (and not merely as a correc tion to the pressure). there will be corresponding terms in all the other thermodynamic variables, a: well as changes to the equilibrium concentrations One major advantage of using the chemical pic ture lies in the possibility to model complicate< plasmas. and to obtain numerically smooth am consistent thermodynamical quantities. In the chemical picture. perturbed atoms mus be introduced on a more-or-less ad-hoc basis t< avoid the familiar divergence of internal partitior functions (see (.9. Ebeling et al. 1976). in othe words. the approximation of unperturbed atom precludes the application of standard statistica mechanics. i.e. the attribution of a Boltzmann factor to each atomic state. The conventional rem edy is to modify the atomic states. 63.9. by cutting off the highly excited states in function ofdensitj and temperature. The MHD equation-ofstate is based on an oc cupation probability formalism (Hummer 8.: Mi halas l988), where the internal partition function Zlm ofspecies s are weighted stuns Eis [\fBT int _ . , Z, _E 111,-...gisexp — i (2 Here. is label state i ofspecies s, and F]... is th- energy and (Ii,q the statistical weight ofthat state The coefficients 112,, are the occupation probabil ities that take into account charged and neut'ra surrounding particles. In physical terms, in” give the fraction ofall particles ofspecies s that can ex ist in state i with an electron bound to the ator or ion. and l — this gives the fraction of thos that. are so heavily perturbed by nearby neighbor that: their states are effectively destroyed. Per turbations by neutral particles are based on at excluded-volume treatment. and perturbations by charges are calculated from a fit to a quantum- mechanical Stark—ionization theory (for details see Hummer 81. Mihalas 1988). The Opacity Project and, with it, the MHD equation—of-state restricts itselfto the case ofstel- lar envelopes, where density is sufficiently low that the concept. of atoms makes sense. This was the mainjustification for realizing the Opacity-Project in the chemical picture and base it on the Mihalas. Hummer. Dapper] equation of state (Hummer 8/, Mihalas I988; Mihalas et al. 1988; Dappen et al. 1988, hereinafter MHD). The Opacity Project is mainly an effort. to compute accurate atomic data. and to use these in opacity calculations. Plasma effects on occupation numbers are ofsecondary in- terest. 2..12. Physical Picture: OPA l. EOS The chemical pictures heuristic separation of the atomic-physics from the statistical mechan- ics is avoided in the physical picture. It starts out, from the grand canonical ensemble of a sys- tem of electrons and nuclei interacting through the Coulombpotential (Rogers 1981b. 1986, 1994). Round clusters ofnuclei and electrons, correspond- ing to ions, atoms and molecules are sampled in this ensemble. Any effects of the plasma environ- ment on the internal states are obtained directly from the statistical mechanical analysis, rather than by assertion as in the chemical picture. There is an impressive body ofliterature on the physical picture. important sources of infortna— tion with many references are the books by Ebel— ing et al. (1976), Kraeft et al. (1986). and Ebeling et. al. (1991). However, the majority ofwork on the physical picture was not dedicated to the problem of obtaining a high—precision equation of state for stellar interiors. Such an attempt was made for the first time by the OPA L—t.eam at Lawrence Liv- ermore National Laboratory (Rogers 1986; lgle- sias & Rogers 1995; Rogers et al. 1996. and refer- ences therein), and used as a. foundation for the OPAL opacities (Iglesias et al. 1987, 1992: lglesias 8r. Rogers 1991; lglesias 8!. Rogers 1996; Rogers 8/. lglesias 1992). The OPAL approach avoids the ad-hoc cutoff procedures necessary in free energy minimization schemes. The method also provides a systematic 29 procedure for including plasma effects in the l ton absorption coefficients. An effective po tial method is used to generate atomic data w have an accuracy similar to single configura Hartree-Fock calculations (Rogers 1981a). In contrast to the chemical picture. plague: divergent partition functions, the physical pic has the power to avoid them altogether. Parti functions ofbound clusters ofpart icles (e.g. at and ions) are divergent in the Saba approach, has a compensating divergent: scattering state in the physical picture (Ebeling et al. 19761Ro 1977). A major advantage of the physical pic is that it incorporates this compensation at outset. A further advantage is that no assu tions about energy-level shifts have to be mad follows from the formalism that there are non As a result. the Boltzmann sutn appearin the atomic (ionic) free energy is replaced by tlu called Planck—Larkin partition function (PLl given by (Ebeling et al. 1976: Kraeft et al. 1 Rogers 1986) Eh PLPP : :9” [exp (— ART) - l + is His erT The PLPF is convergent without additional cu criteria as are required in the chemical picture. stress, however, that despite its name the PLP not a partition function. but merely an auxil term in a virial coefficient (see, e._q., Dappen e 1987). The major disadvantage ofthe physical pict is its formulation in density expansions. EX] sions that first of all are very cumbersome to c out, which means that only terms up till 3 in . sity have been evaluated (Alastuey 8:. Perez 1 Alastuey et al. 1994, 1995). Second, the slow vergence ofthe problem, means that even this traordinary accomplishment has a rather lin‘ range of validity. The chemical picture, on other band, do not need to rely on expansions. complicated expressions, possibly with the cor asyn‘iptotic behavior, can be used freely. 2.2. The Coulomb correction The Coulomb correction, that is, the co quence of an overall attractive binding force neutral plasma deserves close attention, bec. it. describes the main truly non—ideal effect. der conditions as found in the interior of nor- mal stars. Already in a number of early papers (e.g. Berthomieu et al. 1980; Ulrich 1982; Ulrich & Rhodes 1983; Shibahashi et al. 1983, 1984) it was suggested that improvements in the equation of state, especially the inclusion of a Coulomb cor- rection, could reduce discrepancies between com- puted and observed p—modes in the Sun. Respond- ing to this, Christensen—Dalsgaard et al. (1988), showed that the MHD equation ofstate indeed im- proved the agreement with helioseismology. That the largest change was caused by the Coulomb cor- rection was not immediately clear, since the Mle equation of state also incorporates other improve— ments over previous work. From early comparisons between the MRI) and OPAL equations ofstate (Dappen et al. 1990), it turned out, rather surprisingly, that. the net ef— fect of the other major improvement, the influ— ence ofhydrogen and helium bound states on ther— modynamic quantities, became to a large degree eclipsed beneath the influence of the Coulomb— term. In the solar hydrogen and helium ionization zones the Coulomlyterm is the dominant correc— tion to the ionizing perfect gas. This discovery led to an upgrade of the simple, but astrophysically useful Eggleton et. al. (1973) (EFF) equation of state through the inclusion of the Coulomb inter- action term (CEFF) (see Christensen—Dalsgaard 1991; Christensen—Dalsgaard 8x. Dappen 1992). The leading—order Coulomb correction is given by the Debye-Hiickel (DI—I) theory, which replaces the long-range Coulomb potential with a screened potential, as outlined below. 2.2).]. The Debye & Hiickel (1923) theory of elec— trolytes, describes polarization in liquid solutions of electrons and positive ions. This description also applies to ionizing gases. Assuming the parti- cles can move freely, the electrons will congregate around the ions, and the ions will repel each other due to their charges. With their smaller mass and higher speeds, the paths of electrons are deflected by the ions increasing the chance of finding an electron closer to an ion. This screening by the electrons decreases the repulsion between the ions, acting as an overall attractive force in the plasma. The Debye- Hilckrel approximation The fundamental assumption of Debye and I-Iiickel is that ofstatistical equilibrium, according to which the local density of particles of type j (including electrons) immersed in a potential 2,0 around an ion. 2', can be expressed as "1(7'1')=(mlexpl-Zzifiwtil/WT) , (3) where Zje and (DJ) are the charge and mean den— sity of the particles and 7'2_,-(r) are the perturbed densities. e’=(r) is the plasma-potential or the ef— fective (screened) inter—particle potential. Over—all charge neutrality dictates that Zomz-zo <:> (7)9)::(nj)7;j. (4) j #e W'ith these perturbed densities, the corresponding charge density is PM = ZMum—W(“t/W+zteatr.) j = Z Zjetni) [w _ ewe-WT] + 2 jic resulting in the Poisson equation VQMH) = -47T€ ZZj("tile—Z’M'S‘VWT + 21-502) .7' (6) And now comes the most critical of Debye"s ap— proximations: To make Eq. (6) more tractable, the exponential is expanded in a power series, and only terms up to first order are retained. The zero— order term is the net-charge, Eq. (4). Solving Eq. (6) with the remaining first-order terms results in a screened Coulomb potential—the Debye-Huckel potential Ze r “W” -+ (7) hr) = where /\DH is the Debye-length The approximation of disregarding higher order terms affects the low temperature and high den— sity region where the inter-particle interactions be— comes too large to be described by just. the first or- der term. This is a manifestation of the problems with the classical, long-range part. of the Coulomb field in a plasma. Investigations taking the physical picture point ofview indicate that the original potential defined in (6), is a good choice for a plasma potential (Rogers 1981b), and only the truncation ofthe ex— ponential resulting in the Debye-Hiickel potential is oflimited validity (Rogers 1994) At high densities the effect is in fact over- estimated by using the I)ebye-I-Iiickel potential (7). The relative Coulomb pressure in the Debye- Hiickel theory, expressed in terms of the coupling parameter. pDH/kBT : —f3/2/\/1—2, is a negative contribution to the pressure. At very high densi— ties, the over-estimation of the Coulomb pressure can be so severe as to result in a negative total pressure. The negative pressure differences seen in the comparison plots in Sects. 3 and 4. sug~ gests that. the amplitude of the Coulomb pressure is larger in OPAL than in MHD. This statement is true when the r-correction, mentioned below, is applied to the MHD EOS. To get a feeling for the behavior ofthe Coulomb pressure, we use the perfect. gas law to obtain the approximate expression r a R1/3/,l-]/3<22>1/3, (9) where [1 is the mean-molecular weight. This leads us to anticipate differences between OPAL and MHD, stemming from different treatments of the plasma interactions, to increase with R, and that such differences will be somewhat reduced when we mix in helium and metals. 000 41- 4:- 41- The r correction in DH theory As they were investigating electrolytic solutions of molecules under terrestrial conditions, it was natural for Debye and I-liickel to consider elec‘ trolytes made up of hard spheres. Assuming there is a distance of closest approach. rm,n to the ion, Eq. (7) is modified to ZP. e—(T*Tmt'.)/Ar)tt — l 'l’ rmin/ADH 7' for r 2 rm," and constant, d;(rmm), inside, remov— ing the short range divergence. To obtain the free energy, we apply the so—called recharging proce— dure detailed in Fowler 81. Guggenheim (1956) to Eq. (10), and get the result. without rmin, multi- plied by the factor we) (10) T(.’L') : 3[ln(1+ gr) — .13 + gram—3 , (11) I31 where .L' : rmm/ADH. In short. the rechargin, procedure consists of varying all charges in th potential and integrating from zero to full charge Equation (11) is the analytical result of this inte gration and is based on the assumption that rm, is independent. of the charge of any particles. Th r-factor goes from one to zero as .1: increase. re ducing the Coulomb pressure which was overesti mated before. \Vith the T correction we can avoit the negative pressures mentioned above. Craboske et al. (1969) proposed to use . -1 l'.'i/2l7tel] . (l'2 Ft/Qlue) T'mm Z (7)!“J [ART for stellar plasmas, and it was later used in tlr MHD EOS but not in OPAL. This choise of rmin i merely the distance ofequipartition between ther mal and potential energy of electrons approach ing ions. Since the charges are opposite therl are, however, no classical limits to their approach Also notice that since this choice of rm-m depend. explicitly on charge, the recharging procedure wil result in a different form of T. A thorough and critical review of the Debye Hiickel theory can be found in Fowler & Cuggen heim (1956). Chp. IX. and a very clear presen tation is found in Kippenhahn (K: \Veigert (1992) though the latter does not mention 7'. 9.9.9. Other higher-order Coulomb corrections Obviously, the r correction is just one particu lar higher-order Coulomb correction. We can us it as a model for developing more general expres sions, by allowing some liberty in the choice c rmin. Let us begin by asking about the distanc of closest approach for quantum—mechanical elec trons. Heisenberg"s uncertainty relation puts firr limits on how localized particle can be — it i smeared out over a volume the size ofa de Brogli wavelength A : h/p. This de-localization elim nates the infinite charge densities associated wit classical point—particles, and hence the short—rang divergence of the Coulomb potential. Based on that, we can tentatively suggest a dis tance of closest approach which is the combine radii of the electron and ion: -.1;/\e + %<’\ion>' Th diffraction parameter, Vii: between two particle i and j, emerging from a more careful quantun' mechanical analysis implies the use of the De Broglie wavelength in relative coordinates rm... = A.,- = (212/2i1.,i.-Br)"/? oc r-l/Q. (13) where [1,,- is the reduced mass. Comparing the Tefunction with the quantum diffraction modifica— tion in Fig. 5 of (Rogers 1994), we see a similarity in the functional form. The asymptotic behavior differs though: ‘r(.ir) —> 13—1 for .r —) oo in the hard-sphere model, whereas quantum diffraction goes as 171/2. The two functions are very close up to .’L‘ 2 1, though, suggesting that preliminary investigations ofquantum diffraction effects in the MI-If) EOS could be carried out by means of the r-function and a new rm,“ as given by Eq. (13). Dividing Eq. (13) by ADI—I. we find that the cor- rection is now a function of 9 only. That is. going from a hard-sphere model of interactions. to in- cluding quantum diffraction, the factor alleviating the short-range divergence of the Coulomb poten- tial becomes a function of 9 instead of R. Abandoning the hard—sphere ion correction for the benefit of quantum diffraction, still leaves us with only the first. term of the Coulomb interac— tions. Could the higher order terms be represented by r in some form? It turns out that 7' would only fit in a very limited range, and it would be more fruitful to use proper expressions. The present analysis however, shows that the effect of includ— ing higher-order Coulomb terms, is much smaller than has previously been estimated by the MHD EOS. It therefore might be a fair approximation to leave them out: for at least the solar case. 3. The EOS landscape in 0 and T ‘- For this comparison, we have computed MHD EOS tables with exactly the same p/T—grid points as the OPAL-tables (Rogers et. al. 1996). to en- sure that the equation-of-state comparison is not influenced by interpolation errors. We do actually use the respective interpolation routines to access the table-values, but by interpolating on the ex— act gridpoints for identical mixtures, we should not lose precision in the process. We compare tables with three different chem— ical mixtures, successively adding more elements to the plasma: Mix 1 is pure hydrogen, mix 2 a hydrogen-helium mixture and mix 3 is a 6-element mixture that, besides hydrogen and helium, also includes C, N, O and Ne. In Table 1 we list 32 the exact mixtures. both by mass abundance, X,- of chemical element, 2', and as logarithmic num ber fractions relative to hydrogen [NM/V11]. The choice of mixtures is that ofthe currently available OPAL-tables, to avoid interpolations in X and Z In the comparisons of this section. we have omit ted the radiative contributions. The MHD equation of state now includes rel- ativistically degenerate electrons. (Gong et al 20011)) as do the new version of OPAL (Rogers 85. Nayfonov 2002). This. of course. is significant for stellar modeling and important for helioseismit investigations ofthe Solar radiative zone (IClliot k Ix'osovichev 1998.). For the present paper, however it is irrelevant due to the lack of controversy or the subject. and we will therefore limit ourselves to dealing with non-relativistic electrons. All plots ofdifferences in this paper present ab- solute differences. Since the absolute quantities span less than an order of magnitude and as they have quite complicated behaviors. we found that normalizing the differences would confuse more than illuminate. The solar track (also presented in Sect. 4) overlaid on the surface plots is not hid- den behind the surface. so as to give an idea ofthe behavior in otherwise obscured regions. While the MHD tables and the pure—hydrogen OPAL table have the same resolution. Mix-2 anc Mix-3 OPAL tables have three times higher reso- lution both in T and g. This can only affect the comparisons ofthe solar Mix 2 and 3 cases, Sect. 4 where it might introduce some extra interpolation- wiggles in the OPAL—MHD differences. The ta- ble comparisons are all done on the low resolutior grid. For the case of pure hydrogen (Mix 1) we plot the logarithmic absolute pressure, but for the other mixtures we plot the logarithm of a reducec pressure, P/(gT). to make it easier to identify non- ideal effects and the location of ionization zones This choice will ofcourse not affect the differences of the logarithms. Apart from the actual pressure we also investi‘ gate the three derivatives (l91nP dlnP tl / _ __ . / —' ‘8] A/ : M 81119 T' Mi alnT 0' ( ’1 (14 where 71 is the adiabatic derivative often calle< f]. These three derivatives form a complete se‘ TABLE 1 CHEMICAL MIXTURES '2 AND 3 (SEE TEXT) element Ari (7(3) [.Ni / .-'V H] X,‘ (9%) [.N,‘ / Nle H 80.00 0.00000 8000 0.00000 He 20.00 —] .20098 I600 - I .29789 C 0.00 — 0.762643 —3.09693 N 0.00 — 0.223398 -3.69693 0 0.00 — 2.171950 -2.76693 Ne 0.00 —- 0.842053 -3.27923 and fully describe the equation ofstate. 3.1 . Pure hydrogen. we start with the simplest. mixture. that is, pure hydrogen (Mix 1). The case of hydrogen is, however, far from simple, not the least because of its negative ion and molecular species. All in all five species of hydrogen: H, H+, H‘, Hg and H; are included in both EOS. The number of negative hydrogen ions does never exceed a few parts in a thousand compared to the other hydrogen species. Already at mod— erate temperature, they dissociate into hydrogen atoms. Despite its low abundance, H— does have an impact on the electron balance since it. is the only (significant) electron sink. The heavy ele— ments with their low abundances are most affected by this. Apart. from this indirect effect on the heavy elements, the most important feature of the H'-ion is of course its bound—free and free—free opacity, which is the primary source of opacity in atmospheres of G, K and M stars. The positive and neutral hydrogen molecules can be seen in the low-temperature—high—density corner of the tables, where their abundance reaches up to 28% of hydrogen, by mass. At, slightly lower densities, which is of greater as- trophysical interest, these molecules only become important at. temperatures below those considered here. The most important feature in the hydrogen- EOS landscapes of Figs. l—4 is, by far, the ioniza- tion (from atom to positive ion), seen as a curved rift. in all the derivatives. It is hardly visible in the surface-plot of the full pressure (Fig. 'l), but 33 becomes obvious in those of the reduced pressur (Figs. 5 and 9). The OPAL—team introduced the quantity RzTé’g‘]. (15. where T6 : T/l06, as a convenient quantity t describe the approximate g — T—stratification ( many stars. This is clearly seen in Fig. l wher we also plotted three iso-R tracks. bracketing th solar track. In the lower panel of Fig. 1 these isc R tracks also highlight. the main feature in th differences: A sharply rising ridge. bell-shaped i logT and centered around logT : 5.5, lying alon logR 2 0. This ridge is a signature ofdifferences i the pressure ionization. The sign of AP in Pig. tells us that MHD has a more abrupt pressure ior ization than the softer OPAL. The reason for thi difference is still not completely clear. It migl be related to differences in the treatment of th short—range suppression of the Coulomb forces, a mentioned in Sect. 2.2.2 and 2.2.3, or it could b a result of differences in the mechanism of pre.c sure ionization (lglesias (Kr. Rogers 1995; Basu et a I999; Gong et. al. 2001a). We now turn to the logarithmic pressure deriva tives, X0 and X7“: displayed in Fig. 2 and 3, respec tively. In both figures, the ionization zone is easil recognized as the canyon or ridge starting in th low-temperature-low-density corner, slowly bend ing over to follow the solar track and disappear a about logT : 6. ln quite a. large region of the g — T plane bot derivatives are equal to one reflecting that. the ga is a perfect. gas. In this region the differences ar very small (ie. less than 0.03%), confirming tha both the chemical and the physical picture con— verges appropriately to the perfect gas case. At low temperatures X0 and XT are dominated by temperature ionization, which is about an order of magnitude more prominent in XT than in X9. This region is a fairly well known regime and here we can directly compare the two pictures. The dif— ferences are indeed small in this region, less than 1% and less than a tenth of the differences in the high—R ridge. The rise of XT in the low-T—high—g corner is due to Hg—molecules. About 28% by mass, of the Fig. ].— Comparison of long in the two pure hydrogen tables. The upper panel shows the ab— solute value from the MHD EOS and the lower panel shows the difference; OPAL minus MHD. The strange boundaries of the surface simply re— flects the shape of the tables. We also overlay the solar track from Sect. 4 for comparison. On this plot alone we also show iso—R tracks (dotted lines) for long = —2,—1,0, going from low to high densities. 34 hydrogen atoms are bound in molecules in this gion. but at higher densities they quickly press dissociate. The fact that the differences incre while the absolute value decreases indicates t MHD is pressure dissociating faster than OPA The differences are again dominated by sharp ridge at high H, but in contrast to press (Fig. 1), the differences in yo and XT return zero for high temperatures and densities. As pressure, the solar track falls over or climbs H—edge in the middle of the ionization zone, a is traversing the iso—R at log]? t: 0. These high]? differences occur in a reg where there is competition between the Coulo terms and electron degeneracy. This makes interpretation much more difficult. Two possi reasons are the previously mentioned short—rat L or» ‘ .7 ~= -‘.-. Fig. 2.— The logarithmic pressure derivative w respect to density X0 = (8|nP/81ng)T for p' hydrogen in the upper panel, and its differen (OPAL minus MHD) in the lower panel. part of the Coulomb interactions and the changes induced in the internal atomic states by the dense, perturbing surroundings. In the MHD EOS, all energy levels of internal states are assumed to be unaltered by the plasma environment. That is, the effect of the perturba- tion by surrounding neutral and charged particles on the internal state is restricted to a lowering of the occupation probability ofthe given state only. In the OPAL EOS, the net result looks similar, but there the relative stability of energy levels to per- turbations is not merely postulated but the result of in-situ calculations of the Schrodinger or Dirac equation for each configuration of nuclei and elec« trons, based on parameterized Yukawa potentials (Rogers 1981a), as mentioned in Sect. 2.1.2. Looking at 71 in Fig. 4 we immediately no— M“? Fig. 3.— The logarithmic pressure derivative with respect to temperature XT : ((91n P/aln T)(, for pure hydrogen in the upper panel, and its differ— ences (OPAL minus MHD) in the. lower panel. 35 tice how well this quantity displays the ionization zones while leaving out everything else. This prop erty is also reflected in the differences. which here are of about the same magnitude in the ionization zone as in the high—R ridge. The high}? differ- ences have also changed characteristics, changing sign periodically, while retaining the overall bell~ shape in logTofthe amplitude. We mention, how— ever, that at least some of this behavior might be due to the numerical differentiation scheme used in the OPAL EOS (see Sect. 5 and Fig. 25). 3.2. Hydrogen and helium mixture The effect of helium in the thermodynamical quantities is revealed by the addition of 20% he- lium and comparison with the pure—hydrogen case. Fig. 4.— The adiabatic logarithmic pres— sure derivative with respect to density '7] : (fllnP/Bln g); for pure hydrogen in the upper panel, and its differences (OPAL minus MHD) in the lower panel. The first thing we notice from Fig. 5 is how well the reduced pressure P/(gT) reveals all the dis— sociation and ionization zones (except H‘); The Hg—formation in the low—T-high-g corner and the prominent ionization of hydrogen together with the two He ionization zones, the first eventually merging with the H ionization. The effect of de— generate electrons is evident in the high—T—high—g corner. We also notice another thing: while the pure hydrogen OPAL—table was cutting the high—g, low- T corner, leaving a little less for the comparison, the mixture OPAL-tables allow a full comparison since they have the same boundaries as the MHD tables. The slightly larger table reveals a new fea— ture in the differences. For pure hydrogen. the pressure difference drops suddenly in the high—g, low-T corner, due to faster molecule formation in OPAL as compared to MHD. But in the slightly larger tables used for the remainder ofthis section, this difference suddenly goes to zero before it falls down the high-R edge. In the pressure differences, one can just barely identify the first helium ionization zone, whereas the second is too faint to be seen here. The high—R differences are a little smaller than for pure hydro— gen, as anticipated from Eq. (9) and the discussion following it. This can be most clearly seen by com- paring the dip in the hydrogen ionization zone. The addition of helium is also evident in the logarithmic pressure derivatives in Fig. 6 and 7. First we see the deep rift (ridges in Fig. 7) of the hydrogen ionization zone. Then comes a small groove from the first helium ionization zone, a groove which, when it widens and gets shallower at higher densities, eventually merges with the hy— drogen ionization zone, as is the case for the so- lar track. Widely separated from the hydrogen and first helium ionization zones, we find the sec- ond helium ionization zone. It seems to disappear at the low density edge of the table, but that is only so because the. ridge gets very sharp and is unresolved in temperature, at low densities. Hot— ter stars, that is, stars shifted towards lower R, will clearly exhibit three, more distinct ionization zones when compared with the Sun. Apart from the two helium ionization zones, the differences in the pressure derivatives are very sim— ilar to the pure hydrogen case. The high-R dif— ferences are somewhat smaller though, as are the 36 differences in the hydrogen ionization zone. The differences in h. also displays a very small ripple along log}? 2 —4, which might be due to differ— ences in the differentiation technique (see Sect. 5). From the differences in AT (Fig. 7), we see that the absolute differences in the three ioniza— tion zones are just about the same. If we instead compare the differences relative to the size of the respective ionization ridges, we get 0.16% and 3% relative differences for the hydrogen and helium ionization zones respectively. That is, MHD and OPAL have about 20 times better agreement on hydrogen than on helium. In Fig. 8, 7] appears like what We would antic— ipate from the pure hydrogen case in Fig. 4. The first helium ionization zone is only visible at low densities. as it merges with the hydrogen ioniza- Fig. 5.— The reduced pressure, P/(gT), for the H— He mixture in the upper panel, and its differences (OPAL minus MHD) in the lower panel. tion zone shortly before the solar track is reached. The differences, however, exhibit a much more complicated structure. Along each of the ioniza— tion zones, there is a deep valley in the differences. and along the bottom of these valleys runs a very sharp ridge, bringing the differences up to posi— tive values. This is a clear sign of a broad neg» ative peak minus a sharp negative peak, mean— ing that MHD temperature ionize faster than does OPAL. in the beginning of this section, we found that MHD was also pressure ionizing faster than OPAL, so all in all OPAL is the softer EOS ofthe two. The ridge—in—the—middle—of—the—valley picture is also found in the pure hydrogen case (Fig. 4). but as the hydrogen ionization zone is not frilly covered at low densities, the low—T side of the val— ley is missing. Fig. 6.—— x0, the logarithmic pressure derivative at constant temperature, for the H-He mixture in the upper panel, and its differences (OPAL minus MHD) in the lower panel. 37 3.3. H,He,C,N,O and Ne mixture In this section we add the last four elements considered, namely carbon, nitrogen, oxygen and neon. Comparing Figs. 9—ll of this section with the corresponding Figs. 5—7 of the previous sec— tion, hardly any differences appear, neither in the absolute values nor in the differences between the two EOS. For a few points on the high—R. boundary of the tables, differences in M, and ‘71 have increased dramatically. At least some of these odd points are. the same for X? and '7]. This might. indicate that these points are spurious, possibly associated with convergence problems in either EOS in this difficult region. The heavy elements are just barely discernible in the differences of)“, (Fig. ll). However, for 71, 3414" Fig. 7.7 X7“: the logarithmic pressure deriva- tive at. constant density, for the H—He mixture in the upper panel, and its differences (OPAL minus MHD) in the lower panel. in Fig. 12, the heavy elements appear clearly both in the absolute '71 and in the '71 differences. espe- cially along the low—g edge of the table. Between the first and second ionization zones ofhelium (cf. Fig. 8), we notice some wiggles, which are likely resulting from the third ionization zones of car— bon and nitrogen, and the second ionization zones of oxygen and neon. Above the second ionization zone of helium, we can see all the ionization zones from the fourth ionization zone. ofcarbon right up to the tenth ionization of neon, although they are not all resolved in this g—T—grid. A rough estimate reveals that the relative difference between MHD and OPAL for the heavy elements is of about the same magnitude as the one for the helium ioniza‘ tion zones, i.e. 3%, or about 20 times worse than the 0.16% agreement for hydrogen. This unexpectedly large discrepancy for the heavy elements might be a hint that these differ~ >414” Fig. 8.—~ 71 for the H—He mixture in the upper panel, and its differences (OPAL minus MHD) in the lower panel. ences are primarily caused by differences in the lower excited states. For hydrogenic ions, there are analytical solutions for all states. This might. explain the small discrepancy for hydrogen. For ions with more than one electron there are no an- alytical treatments, except for their higher states, which become nearly hydrogenic. So it might well be that the lower lying states of the non— hydrogenic ions are responsible for the differences noticed here. The Yukawa potentials (Rogers 1981a), which are used to describe bound electron states in OPAL, are fitted to give the correct (ex— perimental) ionization energies. MHD uses exper— imental results for the energy levels. It is no sur— prise therefore to get. quasi—perfect agreement on the location ofthe ionization zones (confirmed by the ridge—in—the—middle—of—t.he—valley picture in the '71 differences), whereas the energies of lower ly— ing excitation levels might differ These differences Fig. 9.— Reduced pressure for mixture 3 (cf. Tab. 1) in the upper panel, and its differences (OPAL minus MHD) in the lower panel. propagate into the partition functions and affect the course of ionization. In addition, the differ— ences in the adopted micro—field distribution, and the mechanism by which they ionize highly excited states, might play a role in this region (Nayfonov & Dappen 1998; Nayfonov et al. 1999). Since the differences occur at the low—g edge of the table, we expect however, that they mainly reveal dif— ferences in the thermal ionization, not in pressure ionization. Let us return to pressure and have a closer look at the non-ideal effects in the high-T—high-g cor- ner. From the dotted iso—R lines in Fig. 9. it is clear that the non-ideal effects are not functions of R alone. Instead it turns out that they are largely functions of 92/713. Comparing the perfect. gas pressure and the fully degenerate, non—relativistic electron pressure Fig. 10,— X9: the logarithmic pressure derivative at constant temperature, for the full mixture in the upper panel, and its differences (OPAL minus MHD) in the lower panel. 39 Pperf = and pas-g 21(1)”3 fi( 9 P 5 87r me pem” we see that the two pressures compete along 92 (x ”Fa—lines. This means that relative to high-R (Coulomb) effects, there are more degeneracy ef— fects in the high—Tag corner of the table, which reveals the nature of the sharp rise of both P and Xe in this region. The correlation with larger OPAL—MHD differences (See lower panel of Fig. 9) prompted us to perform a direct comparison between the Fermi—Dirac integrals from the two codes. We found non—systematic differences a re— assuring eight orders of magnitude smaller than the EOS differences we observe in this region. _B 11111,, An alternative explanation could be the lack of electron exchange effects in the MHD EOS. This is a combined effect of Heisenberg’s uncertainty 3414” Fig. 11. XTY the logarithmic pressure derivative with respect to temperature, for the full mixture in the upper panel, and its differences (OPAL minus MI-ID) in the lower panel. relation (Heisenberg 1927) and Pauli‘s exclusion principle (Pauli 1925): Due to the former, electron wavefunctions are extended, but due to the latter, the wavefunctions oftwo close electrons with same spin cannot overlap. This results in the combined wavefunction either having a bulge or a node at the mid-point between the two electrons, giving rise to two different kinds of contributions to the Coulomb interactions. In the fully ionized, weak degeneracy limit, the first—order e—e exchange pres— sure (DeWitt 1961, 1969) is negative and propor— tional to gz/T. Analyzing the differences in solar solar case, we actually find in Sect. 4, that those powers of g and T are the ones best describing the differences above T 2 2 x 106 I\' 12.— The adiabatic logarithmic pressure Fig. derivative, 71, for the six element mixture in the upper panel, and its differences (OPAL minus MHD) in the lower panel. 40 4. Comparisons in the Sun To study the EOS under solar conditions, have evaluated the MHD and OPAL EOS c g —— T track that corresponds to the Sun u: the respective interpolation routines. Obvim this is a simplified comparison, not ofevolutior models of the Sun, but merely of the equation state for fixed solar—like circumstances. As den strated elsewhere, such a simplified procedur welljustified (See 6.9. Christensen—Dalsgaard e 1988). We use the three chemical mixtures from T 1, bearing in mind that Mix 3 has about twice lar metallicity. In contrast to the comparison the previous section, we now include radiative< tributions. This will ofcourse not change the ferences of thermodynamic quantities, since I, formalism use the well—known additive radia contributions ((‘ox & Guili 1968). In all the figures of this section, we notice I the MHD and the OPAL EOS differs very littl temperatures below 20 000 I\' and above 106 K, they differ significantly in between. And tho the differences are small, above 20 000 K they 1 intriguingly systematic. The wiggles in the differences, most notice: in the region between logT=4—4.5, are almost tainly due to the interpolation schemes. T become quite dominant. in the 71 difference. mentioned in Sect. 3, the tabular resolution of OPAL tables for Mix 2 and 3 is about three ti better than that ofthe corresponding MHD tal This means that most of the interpolation wig comes from MHD. The exception is the pure drogen case (Mix 1), where the tables have same (low) resolution and the respective inte‘ lation errors are of the same order. 4.1 . Pure hydrogen Ifwe take a look at the absolute pressure in 13 a), we notice a bend at log?" = 6.4. This m the bottom of the convection zone. Inside the t vection zone, that. is below logT = 6.4, there is abatic stratification of pressure and temperat Le. Vdz<5|ogT> :71—X0 a (9108]) ad 71XT I When the gas is nearly fully ionized, essenti Vad = 2/5, evidenced as the straight—line pal the curve in Fig. 13 a). ln the ionization zones (the outer 0.02RG), Vad is lowered to about 0.11 at logT 2‘: 4.], again clearly evidenced as a de— pression in the pressure curve. Vad comes back to '2/5 at. logT 2 3.76, but this happens at the top of the convection zone where there is a down- ward bend to a. radiative stratification. The depth ofthe convection zone is about 0.28539. and just slightly higher, at a depth ofOQSRC, hydrogen fi— nally gets fully ionized (fewer than I in 105 are still The Solar case r/R 1.0 0.9 0.3 0.5 0.0 20- L A A A . AAAAAL AA.-1 A -A.‘L . - . - _ . a) 154 L it. ‘ ’ 2 ‘0‘. 7 1 _ MHD I 5.‘ ---- OPAL L 1' D . ........... MHD(-r=l) 0d ' 'V Y """" *1 """"" TV * rffi I l- 4 5 6 7 logT r/R 1.0 0.9 0.8 0 5 o 0 0.0041 ‘ A My”. ‘1 A r . b) I 0.002- I.-._ *- : , x I 0.000 ‘ ,."“’ -.\ s : i _9 —0.0024 r <1 1 C .1 i- -0 004 1 _ OPAL-MHD 7 4.006% ---- OPAL-MHD(T=1) L ........... OPAL-MHD(T(hx)) . -o.ooa— . ,,,,,,,,, T ,Hflfl, ,,,,,,,,, V’- 4 5 6 7 logT Fig. 13* The logarithmic pressure along a solar Q,T-tzrack for pure hydrogen. The upper panel shows the absolute values of the MHD (solid line) and the OPAL (dashed line) pressure. \Ne also plot the MHD pressure, using 1' : 1 to show the effect. of omitting this correction (cf. Sect. 2.2.2). These three pressures are indistinguishable unless we look at the lower plot, showing the difference OPAL minus MHD. Here we show, apart. from the normal MHD, also the version with 1' : l, which seems closer to OPAL, and a version where we have halved the argument of 7'. neutral). at a temperature oflogT : 6.3. For co parison. hydrogen is 99.88% ionized in the mid of the convection zone at logT : 6. So, althor it. is reasonable to say that the hydrogen ioni tion zone is confined to the outermost 2% oft Sun. one should also bear in mind the long tail unionized hydrogen that is extending almost to I This tail has especially large effect on the opacity. since in 1 bottom of the convection zone. visual and UV" only bound states can add opac to the constant “background opacity" from el tron scattering. hi the upper plot of Fig. 14 we can actually. the differences between the absolute values of, It is evident that OPAL has a much smoother a broader ionization zone than the somewhat bun" MHD. Turning off the T—correction (dashed line almost centers MHD on OPAL, but the bumps The Solar case 7/)? 1.0 0.9 0.8 0.5 A A - - -.-..I;A--l .A-l L_‘__A a) _ MHD (190- ---- OPAL «I ........... MHD(T:1) 0.354 . ,,,,,,,, ,, ....... ,,-c-.+--f, 4 5 6 7 logT r/R 1.0 0.9 0.8 0 5 0.005— ‘ ‘ i. A I : I \ ll “ b) 0.000 , i—x’i“ ‘ ‘. “TN 1 -0.005— " >2 Q 1 -0.010: 1 _ OPAL-MHD -0.015 3 - - - - OPAL-Mflbhzl) ; ........... OPAL-MHD(-r(sx)) 4.020: ' ,,,,,, , - , are? - , ......... , 4 5 6 7 logT Fig. 14.— The logarithmic pressure derivat with respect to density, X0, along the solar tr; for pure hydrogen. a) the absolute value, b) t difference (OPAL minus MHD). main the same. These bumps were also noticed by Nayfonov 8x. Dappen (1998) and their analysis showed that in the region where log’l’ : 42*52 the bumps are caused by excited states in hydro- gen. In this part of the Sun. hydrogen is 30% in— creasing to 97.8% ionized, so even a small amount of neutral hydrogen can have a significant effect on the EOS. At. logT 2 6.5 we see how degeneracy sets in. increasing X0 towards it"s fully degenerate value of 5/3. In the lower plot, we notice that degen- eracy is accompanied by an increase in the differ- ences. This could be attributed to the MHD EOS not including electron-electron exchange effects, as pointed out in Sect. 3.3. The behavior of XT (Fig. 15) confirms the pic- ture obtained from Fig. 14, that is. Mill) ionizing The Solar case r/R 1.0 0.9 0.8 0.5 0.0 2,2: l - A A AAAAaal A.Asl AJAAAJ A A4 A 2.0; __ MHD 8) i” ; ---- OPAL ; 1'3“. ........... MHD(-r=1) f 1.83 i- s ‘ p x j : 1.4“j :- 1 2{ _~ 1.0 {» fiw? a n 0.8 -l ' V ' V """" I vvvvvvvv 7 """""" T v '- 4 5 8 7 logT r/R 1.0 0.9 0.8 o 5 0 0 0.04011 ‘ ‘ “ - _ OPAL-MHD b) 0.030- ---- OPAL-MHD(T=1) 0.020 - ........... OPAL-MHIXTUJX” g; 0.010- 0.000- -0.010- -0.020-”,, ,,,,,,,,, 1 ,,,,,,, 1 ,,,,,,, — 4 5 6 7 logT Fig. 15,—— The logarithmic pressure derivative with respect to temperature, X7”: along the solar track for pure hydrogen. a) the absolute value, b) the difference (OPAL minus MHD). faster and being more bumpy than OPAL. Hm ever. since the dynamic. range of \T is much larg' than that of \0. the bumps, which have still abm the same size as those of \ 0. are now being dwarf? by the much larger ionization peak in \T. (‘or paring the differences shown in the lower plot we notice that the overall differences are abor twice as large as for \0- but the ionization per in the respective upper plots is about 10 tim« larger for x7 than for NO We also notice that f XT- as a likely result of the higher dynamic rang the interpolation-wiggles at logT S 4.6, are mm more prominent than in No- We can also distinguish MHD from OPAL i the absolute values of 7] (Fig. 16 a)). althoug they are much closer than in the y’s of Fig. l and 15. This is confirmed in the differences show The Solar case r/R 1.0 0.9 0.8 0.5 ( 1.70-, i ‘ ‘ LTTMLMI ‘ ' i i i a) 180—; 1.50-j r1 1.404 1.30-§ — MHD _ _ _- OPAL 1.20— ........... MHD(T=1) 1.10-3, , ,,,,,,,,, , ,,,,,,,,, T ......... ,fi 4 5 6 7 logT r/R 1.0 0.9 0.8 05 ( 0.00414 .l..--1 AA; A A 0002-: 0800-: 4 3 -0.0021 -0.004 1 _ OPAL-MHD _0 008.: - --- OPAL-MHD(T=1) , ........... OPAL-MHD(r(bx)) .( —o.008 , ,,,,,,,,, , ..... fific 4 5 6 7 logT Fig. 16.— The adiabatic logarithmic pressu‘ derivative, 71, along the solar track for pure h; drogen. a) the absolute value, b) the differem (OPAL minus M HD). 42 in panel b), which are overall smaller by an order of magnitude compared to the P—, X0‘ and X7“— differences. In contrast to our experience with P, X9 and XT. here diminishing the r-correction in MHD (dashed and dotted lines) does not lead to any better agreement with OPAL. This is again convincing evidence that 71 is a very efficient filter for high-R effects. The differences that we see are therefore due to the physics of ionizat ion, except at low temperatures, where interpolation errors seem to dominate. 4.2. Hydrogen and helium mixture, The effect of helium is very hard to see in the reduced pressure shown in Fig. 17 a), and in the shape of the differences in Fig. 17 b). However. a The Solar case r/R 1-0 0.9 0.8 0.5 0.0 8.201 1 i ‘4 “‘““““ ”5“? *4 L, S ..+—~, ; 8.10-j :_ A 1 : ‘2. 2 E 8 a 00 , f” s 2 — MHD E 7.90-5 ---- OPAL E. L ........... MHD(-r=1) E 7.80;”v'n vvvvv 1 vvvvvvvvv . vvvvvvvvv 1-; 4 5 6 7 logT r/R 1.0 0.9 0.8 0.5 0.0 a N O ‘1 a— 4 Q q .1 , _ OPAL—MHD 4 - 4.006% - --- OPAL-MHD(T=1) -_ , ........... OPAL-MHD(-r(bx)) . -o.008- . ,,,,,,,, , ,,,,,,,,, , ,,,,,,,,, , L 4 5 8 7 logT Fig. 17.— The reduced pressure in the H—He— mixture along the solar track. a) the absolute value, b) the difference (OPAL minus MHD). The thin horizontal line in panel a), indicates the fully ionized, perfect. gas pressure. comparison with the pure hydrogen case (Fig. 13 allows us none the less to see a few changes to tl differences in the lower panels: The peak arour logT : 4.7 gets considerable smaller by addir helium. except for the r : 1 case. where the di ference actually increases in this region. Also. tl differences outside the high-R region decrease l adding helium. independently of the choice of r in general, adding helium does not alter tl shape of the differences in P, x0 or X7“: and tl changes due to composition are only manifest l a change of the amplitude of the peak arour log7‘ : 4.7. This is surely due to the fact the most of the ionization in the Sun takes place 5 the high-H. region. so that the first—order higl H. differences due to the ionizations themselw simply dwarf the second—order effects due to d. tailed partition functions. among other. The sol: track does follow the ionization zones to some (1' The Solar case 'r/R 1.0 0.9 0.8 0.5 ( __ MHD - - - - OPAL ........... MHD(T=1) 0.005 1 0.0004: -0.0054: >2 1 q ‘l —0.0lO -_ : — OPAL—MHD -0.015-‘ ---- OPAL-MHD(T=1) : ........... 0PAL-MHD(r(bx)) —0.020-' ,,,,,,,,, ,, ------,c 4 5 6 7 logT Fig. l8.— x0 for the H-He-mixture along the st lar track. a) the absolute value, b) the differem (OPAL minus MHD). gree, and only enters the hydrogen ionization zone “head on”. With the solar track curving along the hydrogen ionization zone in this way the ionization features will be smoothed out over a much larger temperature interval than if we had examined an iso-chore. This smoothing leads to more blending of ionization zones from various elements, hamper- ing analysis. The shape, merging and smoothing of the ionization features is best seen in Figs. 9—1 1. This behavior is clearly illustrated in. 6.9.. Xe (Fig. 18 a), where we observe a rather sharp onset of ionization followed by a much slower transition to full ionization. The second ionization of helium appears as part of the bump around IogT : 5. The bump is somewhat more pronounced than in the pure hydrogen case. A more careful comparison with the pure hydrogen case (Fig. 14) reveals the first ionization zone of helium as a slight extension of the hydrogen peak. on the side towards higher The Solar case fi/R 1.0 0.9 0.3 0.5 0.0 2.0{ _ MHD a) i” . _ - - _ OPAL : 1-3‘. ........... MHD(-r=l) ,- 1.6-5 L Q g C 1.4- :' 1.23 L 1 ''''' ' '1 i- 1.0jr “‘T 0.8; ' 1 """"" I """"" I """""" T _- 4 5 6 7 logT ‘r/R 1.0 0.9 0.8 0.5 0.0 0.040: L ‘ A A TAMPA“ A A A A ' . _ OPAL-MHD b) a 0.030- '" , ---- 0PAL—MHD(-r=1) 0.020% ........... OPAL—MHD(T(bx)) f- “ i a... >< 0.010- 3;, I- 9.1,“... 0.000-:LWV‘ ‘\ ," M -0.01o- -' "0.020' 1 I‘Tfivfi-,1 vvvvvvvvv 1vfrv '**T'_ 4 5 6 7 logT Fig. 19.— X7 for the solar track and the H—He— mixture. a) the absolute value, b) the difference (OPAL minus MHD). temperatures. l-lelium gets almost fully ionize at logT : 6.0. where 1.77% is singly ionized an 98.23% doubly ionized. it only ionizes slowly fu ther at, higher temperature, until at. logT : 6.7 it. suddenly becomes fully ionized. This happer at a depth of0.ti3R(.3, at the edge ofthe hydroge burning core. At logT : 6.5, just slightly abm the temperature where hydrogen gets fully ionize< there is finally no more neutral helium left. For X7“ (Fig. 19 panel a). the bump at logT 2 is a clear sign of the helium added. as opposed t the similar but more entangled bump in Xe (c Fig. 14 and 18). The second He ionization zone very distinct in 71 (Fig. ‘20 a). and the first ioniz; tion zone is manifested by a widening ofthe hydrt gen ionization zone towards the high-T side. Tl differences (panel b) are just as entangled as fc pure hydrogen (Fig. 16) but with lower amplitudl On the descending part . just above logT : 5, ther The Solar case r/R 1.0 0.9 0.8 0.5 0 1.70—;‘ ‘ A a) 1.60—; 1.50-j a: 1.40-3 1.30% __ MHD __-_ OPAL 1.20-5 ........... MHD(T=1) ) .3 5| 33' y 51 o: I I u l ’0-004 '2 .' _. OPAL-MHD l —0.006-i .' --_- 0PAL-MHD(-r=1) . ‘,' ........... OPAL-MHD(T(bx)) .4 -o.003~ ”,1 ,,,,, f”, ,,,,,,,,, V ,,,,,,,,, W 4 5 6 7 logT Fig. 20* 71 for the solar track and the H-He mixture. a) the absolute value, b) the differenc (OPAL minus MHD). are some large interpolation errors, caused by the change to the coarser grid. We also notice a pecu— liar bump at. log?" : 6.6 Looking at. the various difference plots in this section, we see a. correlation between a high am- plitude in the differences and a high R-value, a property we already inferred from the solar track (Fig. 1). The minimum in R is found at the base of the convection zone, where we also find a lo- cal minimum in the magnitude of the differences between the EOS. The location of this local min- imum coincide for all four thermodynamic quan- tities. This confirms our suspicion that: at least. some of the discrepancy stems from T. The reason for this conjecture is that. the differences between MHD EOS with different 7' almost vanishes in this region, whereas they increase in the same way as The Solar case r/R 1.0 0.9 0.8 0.5 0.0 A AAAAIAALLI ‘ AA- 1 9 iv Ci 1 W “I l°g(P/pT) 9° 8 “.1“ . _ MHD ---- OPAL I- ........... MHD(-r=l) D. M .2 < _ OPAL-MHD -o.003- ---- OPAL-MHD(-r=l) »_ 1 ........... OPAL—MHD(-r(ax)) t —o.0033 fl, , q-- 4 5 3 7 logT Fig. 2].—— Reduced pressure for the solar track and the 6—element mixture. a) the absolute value, b) the difference (OPAL minus MHD). The thin horizontal line in panel a), indicates the fully ion- ized, perfect. gas pressure. 45 the MHD—OPA L differences grow for intermediate temperatures. At high temperatures, above a minimum oc— curring at logT 2 6.4, the MHD—OPAL differ- ences grow, but the differences between the three 7' versions themselves remain small. On the solar track, logR : —1.8 at the minimum of the MHD- OPAL difference, and it only rises slightly to -1.4 at logT : 6.8 where the solar track bends to follow more or less an iso-R line. The constant R value is attained around logT : 6.15. The differences be tween the T-versions are indeed the same in both of these regions (this is best seen in the pressure differences eg. Fig. 17), which explains why the three curves with different 7' follow each other so closely at high temperatures. The MHD—OPAL difference in this region can therefore not be ex- plained by the r—correction. it also turns out that The Solar case r/R 1.0 0.9 0.3 0.5 0.0 1.04.: I a A A A . I A l xl A a 1.025 a) g 1 I 1.00 1* .. 7 0.935 :_ >2 0.93% 3 0.94—: _ MHD :- 0.92-j ---- OPAL :_ 090% ........... MHD(T=1) -_ 0'88; ' I I V L 4 5 6 7 logT r/R 1 0 0.9 0.3 0 5 o 0 0.005- ‘ I‘, ‘ “ ‘ - ‘ p . ,' “ b) : 0.000 —-—-A ' ‘ ATVN I """" C —o.0051 _e 9 r . <1 3 : —0.0101 _- j _ OPAL—MHD _ —0.015—‘ ---- 0PAL-MH0(-r=1) '— 3 ........... OPAL-MHD(T(»X)) t -0-020:---,.. - ,W hi- 4 5 3 logT Fig. ‘22.— The logarithmic pressure derivative with respect to density X9 for the 6—element mix— ture along the solar track. a) the absolute value, b) the difference (OPAL minus MHD). in this region, the differences of the g-T plane are mainly functions of 92/7“3 instead of R. In Sect. 3.3 we suggested that. this dependence might arise from electron exchange effects or maybe from pos- sibly different. evaluations of the Fermi~Dirac in— tegrals. However, a third explanation might be based on the quantum diffraction mentioned in Sect. 2.2.3. 4.3. H,He,C,N,O and Ne. mixture Adding C, N, O and Ne to the H-He mixture has two main effects: first, it displaces 4%1—1e(cf. Tab. 1), thereby diminishing the helium features, and second it leads to a slight decrease in the high—R OPAL-MHD differences due to the increased mean charge [see Eq. (9)]. Only in '7'] (Fig. 24), can the heavy elements be observed directly. Comparing The Solar case r/R 1.0 0.9 0.3 0.5 0.0 2.25' ‘ 9 “H“— 2.05 _ MHD a) _ . - --_ OPAL 1-8‘: ........... MHD(T=1) ‘ 1.65 >‘< : 1.4 5 1.25 1.05 0.35 ' - vw 4 5 6 7 logT 'r/R 1.0 0.9 0.3 0.5 0.0 0.040- 1 . A .A..-Aa1...l.....l.. A -- 0 o __ OPAL—MHD b) . 30- ---- 0PAL-MHD(-r=1) 0.020 _ l." fl.- ........... OPAL_MHD(T(HX)) - 5 PI I "1'3 3 0.010 ‘ 4."; l‘ \ ' it?!" y. : '.-.. ................... .J‘ 0300-:WV ‘ I, M -o.010- ‘ —0.020- f,“ ,, , -, ......... 5 6 7 logT Fig. 23.-— The logarithmic pressure derivative with respect to temperature XT for the 6-element mixture along the solar track. a) the absolute value, b) the difference (OPAL minus MHD). 46 with the H—He case (Fig. 20), and going from low to high temperatures, we first notice a slight di- minishment ofthe feature associated with the first ionization zone of helium due to the 4% decrease ofthe helium content. This weakening of the l-le+ feature is counteracted by the second ionization zone of carbon (24.38eV), as well as that of the less abundant Ne (21.56eV) and N+ (29.60eV). The feature of the second ionization zone of he— lium is also diminished, but counteracted by the third ionization zone of oxygen O++ (54.93eV). the most: abundant heavy element. C++, 03+. N++ and Ne++ adds further ionization in this tem— perature region. Continuing towards higher tem— peratures we notice a slight, straightening of the “knee” around logT 2 5.3, due to the intermedi- ate ionization stages of C, N and Ne with ioniza— tion potentials between 47eV and 240eV. Finally, The Solar case r/R 1.0 0.9 0.3 0.5 0.0 1.70-‘ ‘ ‘ CW 160- ' 1.50- ' s 1.40- ' 1.30- — MHD ' _ - - - OPAL 1.20- ........... MHD(T=1) ' 1.10- ,-. -.-- ,, -- -- 4 5 6 7 logT r/R 1 o 0.9 0.3 o 5 0 0 0.0045‘ ‘ A A“ : I” ~ 0.0025 7 0.000 - _ i ‘i I 3 4.0025 ‘5. r —o.0045 _ OPAL-MHD 7 -0.003-: - - - - OPAL-MHD(T=1) -_ , ........... OPAL-MHD(1-(ax)) I —0.006-‘,,_ ......... , ......... H; 4 5 6 7 logT Fig. 24— The adiabatic logarithmic pressure derivative, 7], for the six-element mixture along the solar track. a) the absolute value, b) the dif— ference (OPAL minus MHD). at logT 2 6.2, we find a broad dip supplied by the two uppermost ionization stages ofC,N,O and Ne, having ionization energies in the range between 400 eV and l 360 eV. The only quantity in which the introduction of heavy elements is manifested directly is '71, which is an important key variable in helioseismology (since it is closely related to adiabatic sound speed cg = 71p/g). The promise of these features is that the presence of heavy elements is well marked in '71. Actually, this marking is so distinct (Gong et al. 2001a), that in future solar and stellar appli— cations of the MHD and OPAL equations of state it might be worth to include more heavy elements The influence due to our small quantity of heavy elements is about three times larger than the dif— ference between OPAL and MHD, though we has— ten to add that our heavy element abundance of Z : 4% is chosen too high in order to exhibit the effects more clearly; they would ofcourse decrease with a more solar metallicity around Z = 2%. We have not discussed radiation pressure yet, merely because of the lack of controversy about it. However, it is worth a few notes. The ra— tio between radiation pressure and gas pressure is constant along iso—R lines the two being equal around logR 2 —4.5. The largest radiation ef— fects therefore occur at logT = 6.4 where there is also the smallest discrepancy between OPAL and MHD. The relative size of the effect of radiation is: 0.0007, —0.001,0.003, —0.002 for logP. X9: X7“ and 71, respectively. 5. Discrepancies due to differentiation A closer inspection ofthe derivatives in the per- fect gas region reveals some discrepancies which are likely due to the numerical differentiation per— formed in the OPAL EOS. This is most noticeable in 71, where the OPAL—MHD differences in the perfect gas region are as large as 0.03%, which admittedly is small indeed. Helioseismology will, however, soon be dealing with such precision. This difference most probably comes from problems in the numerical calculation of an adiabatic change as performed in OPAL (note that MHD uses es— sentially analytical expressions for 7], X9 and XT Since an adiabatic change is not rectangular in the T — 9 plane, such an interpretation is consistent with the fact that the differences in the derivatives 47 with respect to g and T (X0 and XT» respectively) are about one order of magnitude smaller. This also means. that in the ionization zones where pressure and entropy are non—linear functions of g and T, this differentiation noise must be much larger. On the. other hand, the differences be— tween OPAL and Mle are still at least an order ofmagnitude larger than this differentiation noise. We hope, however, that future improvements will make OPAL and MHD converge to within that level of the actual EOS, requiring higher numeri- cal standards. The differences in X0 and x7 have a tendency to follow iso—R tracks. while the differences in 7] follow isotherms. These two behaviors are still un— explained. ln Fig. 25, the differences following isotherms are pretty clear, but the iso-R differ- ences are also visible, well below the rising moun- tain at high R—values. Fig. 25* This is a zoom—in on the fully ionized, perfect gas region of a pure hydrogen plasma (cf. Fig. 4), where 71 = 5/3. The upper panel shows the results for the MHD EOS which uses analytical expressions for all first— and second-order deriva— tives. The lower panel shows the same for the OPAL EOS, where derivatives are calculated nu— merically on a grid that are much denser in g and T though, than in the tables published. 6. Conclusion The present comparison of the two MHD and OPAL EOS has revealed the reasons of several dif— ferences between these equations of state. They can be summarized as follows (in order of impor— tance): a) We find the largest. differences at high den- sities and low temperatures, or more pre- cisely, at high R—values. From Sect. 2.2.1 and Eq. (9) we know that. this property is indicative of differences in the treatment of plasma interactions. Comparing the peaks of the differences in (2.9. pressure (See Figs. 13, 17 and 21), we obtain Mix-l-to-Mix—‘Z ratios of 0.797, and Mix-1—to-Mix-3 ratios of 0.788, which agrees very well with Eq. (9). and thus further substantiates our interpre tation. These differences are lowered dra— matically when we put. 7' : l in MHD, indi- cating that it is worthwhile to abolish 7' and reconsider how to get. rid of the short-range divergence in the plasma-potential (See Sect. 2.2.3). b) 1n the higl'i-temperature—liigh—density corner ofthe tables we observe how degeneracy sets in. Along with degeneracy, we also notice how some specific differences are growing. This effect. could be due to quantum diffrac— tion or exchange effects, both included in OPAL but not in MHD. Quantum diffrac- tion is the effect. ofthe quantum mechanical smearing out of, primarily, the electron due to it’s wave nature. The exchange effect is a modification of the quantum diffraction aris- ing from the anti-symmetric nature of two— particle wavefunctions of fermions. c) Differences also appear in the ionization zones, and a great deal of them can be at- tributed to the 1' correction, but not all of it. The causes for the rest. of these differ— ences are not easily identified. They might be due to the basic differences between the physical- and the chemical approach to the plasma. The treatment of bound state en— ergies and wave functions might have an effect in this region. These are highly ac- curate in MHD but. calculated in the iso— lated particle approximation, whereas they 48 are approximate (fitted to ground—state en- ergies), but varying with the plasma. envi- ronment in OPAL. We have also tried ex— perimenting with the assumed critical field strength used in MHD for the disrupt ion ofa bound state [Eq. (4.24) of Hummer & Miha— las (1988)]. However, this intervention had only a very small effect. Earlier investiga— tions by lglesias 8'. Rogers (1995) indicated that a change in the micro-field distribu— tion might have a greater effect, and that highly excited states are more populated in the OPAL EOS, although the OPAL EOS ionizes the plasma more readily than MHD (Nayfonov 8.". Dapper] 1998). d) The evaluation of thermodynamic differen— tials is done numerically in OPA L but ana- lytically in MHD. For the quantities we have examined here, the difference becomes most apparent: in 7]. 1n the trivial perfect gas re— gion of the g—T plane, OPAL is rugged on a 0.03% scale (see Sect. 5), as opposed to the smooth MHD. These 0.03% may sound neg— ligible, but helioseismology is approaching that level. In ionization zones, the discrep- ancies due to differentiation are most likely larger. On the other hand, physical differ— ences between the two EOS are still at. least. an order of magnitude larger. For helioseismic studies of the equation of state it. is a very nice property of the Sun that high— R conditions are found exclusively in the con- vection zone, where the stratification is essen- tially adiabatic, and therefore virtually decoupled from radiation and the uncertainty in the opac— ity (Cltristensen-Dalsgaard & Dappen 1992). As opacity calculations are still subject to errors of 53—10%, we stress the importance of the fact that opacity effects do not contaminate the structure of the convection zone. This means that the solar convection zone is a perfect. laboratory for investi- gations ofthe most. controversial parts ofthe EOS. The difference between 71 from and EOS and that of the Sun can be inferred from helioseis- mology, and that, with an accuracy that by far exceeds the discrepancy between the two of the best present. EOS for stellar structure calcula— tions (Christensen—f)alsgaard et al. 1988). The pursuit for a better EOS is therefore not at. all academic, and we can improve both solar mod— els and atomic physics in the. process (Basu (KI. Christensen-Dalsgaard 1997). We thank Jorgen Christensen-—Da|sgaard For supplying us with a copy of his solar model S Christensen—Dalsgaard et. al. (1996). RT. ac— knowledges support. from NASA grants NAG 5 0563 and 5—l‘2450, and NSF grant 0205500. W. D. acknowledges support from the grants AST- 9618549 and AST-998739l ofthe National Science Foundation. ln addition. VV".D. and RT were sup— ported in part by the Danish National Research Foundation through its establishment of the The— oretical Astrophysics Center. REFERENCES Alastuey, A., Cornu, F., Perez, A. 1994. Phys. Rev. E 49, 1077 Alastuey, A., Cornu, F.. Perez. A. 1095. Phys. Rev. E 5], 17°25 Alastuey,A., Perez, A. 1992, Europhys. Lett. '20, 19 Basu,S., Christensen-Dalsgaard..l. 1997. A8'A 322, LS Basu,S., Dappen,V/V., Nayfonov,A. 1900, ApJ 518, 985 Berrington, K. A. (ed) 1997, project”, Vol. 2. Institute of Physics Publishing “'l‘he opacity Berthomieu, (3., Cooper, A., Cough. D., Osaki, Y., Provost,.l., Rocca, A. 1980. In: H. A. Hill,H. A. and W. Dziembowski.VV. (eds), “Nonradial and nonlinear stellar pulsation”, Vol. 1‘25 of Lecture Notes in Physics, lAU Coll. 38, Springer Verlag, Berlin, 307 Cauble,R., Perry,T.S., Rach,D.R.. Rudil,l<.S., Hammel,B.A., Collins,G.W., Gold,D.M., Dunn,.l., Celliers,P., Silva,l..R.D., Fo- or(J,M.E., Wallace,R..l., Stewart,R.F.., , W'Oolsey,N.C. 1998, Phys. Rev. Letters 80, 1248 Christensen-Dalsgaard,.l. 1991. In: D. O. Gough,D.O. and .l. Toomre,.l. (eds), “Chal- lenges to theories of the structure of moderate- mass stars”, Vol. 388 of Lecture Notes in 49 Physics. lAU Coll. 38. Springer Verlag. l l l Christensen-Dalsgaard, .l.. Dappen, W. A&AR 4(3), 267 Christ ensen-Dalsgaard. .l ., Dappe AjukovS. V. Anderson. E. R., Antia, Basu. S. Raturin, V". A.. Berthomh Chaboyer, B, Chitre. S. M.. Cox, Donatowicz. .l ., l bowski. VV”. A.. Gabriel. M.. Cough. Cuuenther, D. D. Cuuzik. .l. \. Harvey. Hill. F., Houdek. C... lglesias.C. A.. vichev. A. (L, Leibacher. .1. VV’., litt, P. M.C. R.. Provost..l., Reiter,.l., R .lr.. F.. .l., Rogers, F. .l.. Roxburgh. Thompson, M. .l., Ulrich, R. K. 1906. S 272(5266),1286 Demarque. P., Christensen-Dalsgaard. .l., Dappen, VV., ton.Y. 1988. Nature 336, 634 Cox..l. P., Cuili.R.'l‘. 1968.. "Physical princi Vol. I ofPrinciples ofStellar Structure. G and Breach, Science Publishers Dappen. W. 1992. Rev. Mex. Astron. Astrol l4l Dappen, VV’. 1906, Bull. Astr. Soc. India 24. Dappen. VV.. Anderson. L., Mihalas, D. 1987 319, 195 Dappen, VV’.. Lebreton, V.. Rogers, F. 1990. Physics l28, 35 Dappen,W., Mihalas.D., Hummer,D.G., las, R. W. 1988, ApJ 332,261 Debye. P., Hiickel, F,. l9‘23, Physic. Zeit. 24(S DeVVitt , H. E. l96l, Plasma Phys. ‘2, '27 J. Nucl. Energy, P: DeVVitt,H.E. 1969. In: S. Kamar,S. (ed), luminosity stars”, Gordon and Breach, York, 2] l Di Mauro, M. P.. Rabello-Soares, M. C., 384, 066 Christensen-Dalsgaz Basu, S. 200?. Ebeling, W., Forster, A., Fortov. V. 19.. Gryaznov, V. l'\'.. Polishchuk, Y. A. 1.991. “Thermodynamic properties of hot. dense plasmas”. Teubner, Stuttgart. Germany Ebeling, VV., l\'rae1't,W., Kremp. D. 1976. “Theory of bound states and ionization equilibrium in plasmas and solids”. Akademie V’erlag, Berlin. DDR. Eggleton, P. P., Faulkner..l.. Flannery. R. P. 1973. A&A 23, 325 Elliot,J.R.. Kosovichev,A.('1. 1998. In: S. Ko- rzennik,S. (ed), “Structure and dynamics of the interior of the sun and sun-like stars". SOHO (i/CONC 98 Workshop, 453 (in press) Fowler, R... (flnggenheimdi. A. 1956. “Statistical thermodynamics”. Cambridge University Press Gong,Z., Dappen.VV., Nayfonov.A. 2001a. ApJ 563, 419 Gong,Z., Dappen.VV".. Zeidal. 2001b. ApJ 546. 1178 Graboske,H.C., Harwood. D..l., Rogers. F..1. 1969, Physical Review 186(1). 210 Heisenberg. W. 1927, Forsch. und Fortschr. 3(11). 83 Huebner. W'. F. 1986. “Physics ofthe sun”. Vol. 1. 33. D. Reidel Publishing Co., Dordrecht Hummer,D.(l., Mihalas,D. 1988. ApJ 331. 794 lglesias,C. A., Rogers,F..1. 1991, Ap.) 371,1.73 lglesias,C. A., Rogers. F..l. 1995. ApJ 443, 460 lglesias,C. A.. Rogers, F..l. 1995, ApJ 443. 460 1g1esias,C. A., Rogers, F..1. 1996, ApJ 464, 943 lglesias.C. A., Rogers, F..1., Wilson,R. (,1. 1987, Ale. 322, 1,45 lglesias,C. A., Rogers, F..1.. VVilson,R.C.. 1992. ApJ 397, 717 Kippenhahn, R.., V/Veigert, A. 1992, “Stellar struc— ture and evolution”. Springer Verlag. Chp. 18.4 Kraeft,VV., Kremp,D., F31mling,VV., R6pke.G. 1986, “Quantum statistics 01‘ charged particle systems”. Plenum, New York Mihalas. D.. 1)appen.VV'.. l-lummer.D.(T. 1988. ApJ 331.815 Nayfonov. A., Dappen, W. 1998. ApJ 499. 489 Nayfonov. A.. Dappen. VV'.. Hummer. l). (51.. Miha- 1as.D. 1999. ApJ 526. 451 Pauli,VV". 1925. Zeits. 1. Physik 31. 765 Rogers. F. 1981a. Phys. Rev. A 23(3). 1008 Rogers, F. 1994. In: ("1. ('habrier.(il. and F.. Schatzn'lan. 1'3. (eds). “The equation ofstate in astrophysics". 1A1." ("011. 147. (‘ambridgc Uni- versity Press. 16 Rogers. F..1. 1977. Phys. l.ett. 61A. 358. Rogers.F..l. 1981b. Phys. Rev. A 24. 1531 Rogers, F..1. 1986. Ap.1 310. 723 Rogers. F..1.. lglesias,C A. 1992. Ap.1S 79.507 Rogers. F..1.. Nayfonov. A. 2002, ApJ (submitted) Rogers. F..1.. Swenson. F..1.. lglesias.C. A. 1996. ApJ 456. 902 Saumon. D.. Chabrier. ('1. 1989. 1n: (1. VV’egner,C-. (ed). “VV'hite dwarfs". 1AU (_‘011. 114, Springer V’erlag. Berlin. 300 Saumon,D., (.Ihabrier,(‘}.. Horn,H.1V~1.V". 1995, ApJS 99. 713 Seaton..\=1. 1987. .1. Phys. 13 20. 6363 Seaton.lVl..l. (cd.) 1995. “The opacity project”, Vol. 1. Institute of Physics Publishing Seaton.l\’1.J., Zeippen.("..l., Tully,.1.A., Prad- han. A. 1\'.. .V’1endoza.C.. Hibbert, A., Berring— ton. l\'. A. 1992. Rev. Mex. Astron. Astrofis. 23. 19 Shiballashi,l1.. Noels, A.. Ctabriel,M. 1983, A&A 123. 283 Shibahashi, H., Noels, A., Gabrie1.M. 1984, Mem. Soc. Astron. ltal. 55 163 Ulrich,R. 1982, ApJ 258, 404 Ulrich.R.., Rhodes, E. 1983, ApJ 265, 551 This '2-column preprint was prepared with the. AAS ”TEX macros v5.0. 3.3 Pressure Ionization in the MHD Equation of State In this section I elucidate the behavior and importance of an often quoted and ques- tioned term in the Mihalas-H1.111‘1mer-Dappen (MHD) equation of state (EOS), the so—called \Il—term, which provides pressure ionization of neutral plasmas in the high density/low temperature region (Trampedach & Dappen 2004b). 3.3.1 Introduction In stellar evolution computations. and in particular in the case of stars more massive than the Sun, it is generally sufficient to use a simple equation of state. The plasma of the stellar interior is treated as a mixture of perfect gases of all species (atoms, ions, nuclei and electrons). and the Saha equation is solved to yield the degrees of ionization or molecular formation. In the case of low mass stars, however, non—ideal effects, such as Coulomb interactions, become important. For such stars, the most useful equations of state. as far as their smooth real— ization and versatility are concerned, are (i) the so—called Mihalas—Hummer-Dappen (MHD) equation of state (Hummer & Mihalas 1988; Mihalas et al. 1988; Dapper] et al. 1988), and (ii) the OPAL equation of state, the major alternative approach developed at Livermore (Rogers 1986; Rogers ct al. 1996, and references therein). A brief description of these two equations of state is given in Sect. 3.3.3. Although the MHD equation of state was originally designed to provide the level populations for opacity calculations of stellar envelopes, the associated ther- .51 modynamic quantities of MHD can none the less be reliably used even for cores 0 low—mass stars (Chabrier & Baraffe 1997). Low—mass stars harbor the most extreml plasma-conditions, the gas being far from ideal, due to the high densities and lov temperatures. There the Coulomb interactions between particles are more impor tant, as is the destruction of the more fragile excited states. This latter phenomenor will eventually lead to pressure ionization, which in the MHD equation of state i: achieved through an occupation probability formalism for bound species, as detailec in Sect. 3.3.2. The neutral—neutral interactions employed are, however, not strong enough to pressure—ionize a plasma consisting of mostly neutral particles (low tem perature), and approximate higher-order terms are, therefore, included through the \IJ—term, addressed in the present paper. Since the MHD equation of state otherwise includes processes important for low- mass stars and envelopes of white dwarfs (W. Stolzmann, private cormnunication) 6.9., Coulomb pressure and electron degeneracy, the question of the impact anc validity of the \Il—term is more than academic. 3.3.2 Pressure Ionization in MHD In the MHD equation of state pressure ionization is facilitated by occupation proba- bilities, wi, truncating the otherwise divergent partition functions 2 : 2.0.9.6454“ . (38‘ 52 where i labels the levels and E,- and g, are the corresponding energies and statisti weights. The occupation probabilities are then determined by the physical mode the pressure ionizing interactions. As described in Hummer & Mihalas (1983), the formalism incorporates pert bations by fluctuating fields from charged particles, and a first—order approximat to the hard-sphere interactions with neutral particles. It turned out that the apprc mate neutral-neutral interaction alone is too weak to overcome the steeply increas energy of the free, degenerate electrons, and thereby cause pressure ionization. Thi hardly surprising. If, however, there are just a few charged particles present, the fl tuating micro-fields take over, and pressure ionization occurs smoothly. The press1 dissociation of hydrogen molecules avoids this problem altogether, as no electrons .- released, and pressure dissociation does occur in the MHD equation of state, unaic by extra terms. The hard-sphere model, or excluded volume model, has the free energy :2. 1:, Z 487": N,ln 1— (3%) Z .N_,(R.: + 12,)3 (: i i .‘i where both sums extend over all states of all species of particles. Expanding 1 logarithm in 52' = (3%) ZNj1Ri—1— H.713» (. 53 B :7: : [I 9 Z I 8 z. I if z. o o u . The first-order term is already included in the expression for the occupation prob bilities, 10,, apart from a factor of two. As shown by Fermi (1924), occupation probabilities. 10,-, as used in the MH formulation for partition functions, Eq. (38), are accompanied by a term in t1 Helmholtz free energy Fm : f — Zn,- (8f/3n,) , (4‘ where lnw, E — (EU/872,) //cBT. In the case that the non-ideal term, f, is linear in t1 occupation numbers {72,}, the Fw-term vanishes identically, making Fw independe: of the state of excitation. By confining effects that depend on the population excited levels to such linear terms in the occupation probabilities, the equation state is reduced to a [VI >< M~matrix problem, where 11/] 2 170 is the number species (atoms, ions and molecules), instead of a N x N—matrix problem, whe N 2 16 000 is the total number of states considered. The first term in the expansio Eq. (41), of the hard—sphere free energy therefore result in no extra terms in the frc energy. Retaining the N X N-order of the problem, the second—order term can therefo only be included in an approximate way, 6.9., assuming that all particles of a give 54 species are in the same state. Using the ground—state, Eq. (40) simplifies to 4- . €11 : (3%) Z Nlu1R1u ‘1' Rlulsv 1 ' u where 1/ and ILL label particle species. The second-order term is then of the form F, : —,lJBTZfV,/(—Q€3) : —1cBTZN.1n\II.. 111111.: wei. t This artificial term was introduced to allow computation of numerically con: tent results for very cold and dense plasmas. in order to compute tables of th modynamical quantities in rectangular temperature-density grids with less regard physical accuracy in this difficult region of limited astrophysical interest. It is a legitimate concern that the presence of the artificial pressure ionizat? mechanism, the so—called \II-term, might contaminate the results for less extre conditions relevant for stellar interiors. The present analysis shows that this is I the case, and that the pressure ionization in the MHD equation of state for st lar conditions, is caused by the decreasing occupation prol')ali)ilities with increas? density. The \Il-term only affects the high density low temperature corner of the T)-plane, and was merely introduced to ensure numerical stability and convergei in this region. Before presenting the results of the study, which is based on a systematic SWltt ing on and off of the \Il—term, for the convenience of the reader, the specifications 55 the MHD equation state and its relation to alternative formalisms is summarized ir Sect. 3.3.3. 3.3.3 A brief review of the MHD equation of state Historically, the MHD equation of state was developed as part of the internationa. “Opacity Project” (OP, see Seaton 1987; Seaton et (ll. 1992). It was realized in the so- called chemical picture, where plasma interactions are treated through modifications. of atomic states, 216. the quantum mechanical problem is solved before statistica mechanics is applied. It is based on the so-called free—energy minimization method This method uses approximate statistical mechanical models (for example the non- relativistic electron gas, Debye-Hiickel theory for ionic species. hard-sphere atoms tc simulate pressure ionization via configurational terms, quantum mechanical model: of atoms in perturbed fields. etc). From these models a macroscopic free energy is constructed as a function of temperature T, volume V. and the concentrations N1,...,NM of the [VI components of the plasma. The free energy is minimizee subject to the stoichiometric constraint. The solution of this minimum problem ther gives both the equilibrium concentrations and, if inserted into the free energy and its derivatives, the equation of state and the thermodynamic quantities. More specifically, in the chemical picture, perturbed atoms must be introducec on a more—or-less ad-hoc basis to avoid the familiar divergence of internal partitior functions (see 6.9., Ebeling et al. 1976). In other words, the approximation of un- perturbed atoms precludes the application of standard statistical mechanics, i.e. the attribution of a Boltzmann—factor to each atomic state. The conventional remee of the chemical picture against this is a modification of the atomic states, 6.9., 1 cutting off the highly excited states as a function of density and temperature of t1 plasma. Such cut-offs, however, have in general dire consequences due. to the discre nature of the atomic spectrum, t.€., jumps in the number of excited states (and thi in the partition functions and in the free energy) despite smoothly varying extern parameters (temperature and density). The occupation probability formalism er ployed by the MHD equation of state, however. avoids these jumps and delivers ve: smooth thermodynamic quantities. Specifically, the essence of the MHD equatic of state is the Hummer & Mihalas (1988) occupation probability formalism, whie describes the reduced availability of bound systems when immersed in a plasm Perturbations by charged and neutral particles are taken into account. The neutr contribution is evaluated in a first—order approxil‘nation to hard—sphere interaction which is an adequate description for stars in which most of the ionization in the i: terior is achieved by temperature. For colder objects (brown dwarfs, giant planets higher-order excluded-volume effects become important (Saumon & Chabrier 199 Saumon et al. 1995). In the common domain of application of the Saumon et e (1995) and MHD equations of state, Chabrier 8; Baraffe (1997) showed that the tv developments yield very similar results. Despite undeniable advantages of the physical picture, the chemical picture a proach leads to smoother thermodynamic quantities, because they can be writte as analytical (albeit complicated) expressions of temperature, density and partic abundances. In contrast, the physical picture is normally realized with the unwiele chemical potential as the independent variable, from which density and number abun- dance follow as dependent quantities. The physical-picture approach involves, there— fore, a numerical inversion before the thermodynamic quantities can be expressed in their “natural” variables: temperature, density and particle numbers. This increases computing time greatly, and that is the reason why, so far, only a limited number of OPAL tables have been produced, and then only tables that are suitable for stars more massive than N 0.8 [146). 3.3.4 Examination of the influence of the \11 term To quantify the effects of the \II—term on the MHD EOS. two tables were calculated, including and neglecting the 1Il-term, respectively. The pressure difference between the two is plotted in Fig. 4. The occupation probabilities in the MHD EOS are very dependent on the pres- ence of charged particles in order to ionize. At sufficiently large densities and low temperatures there are no such seed charges from temperature ionization, and pres— sure ionization due to the occupation probabilities cannot occur. Even the slightest amounts of metals, i.e., some elements with very low ionization potentials, will supply seed-charges in a larger region than will a pure hydrogen plasma. Introducing helium will, of course, work in the opposite direction due to its high ionization potential. To investigate the most extreme case (fewest seed charges) of astrophysical relevance, the calculations were performed for a hydrogen/helium mixture with Y = 20% (by mass), including no metals. 58 10g(,0/1g (rm—3]) Fig. 4. lower left to upper right are tracks of stellar structure for stars in the mass range 0.6— 1.5]VIQ, as indicated above each track. The dotted contour lines mark the hydrogen and second helium ionization zones respectively, with contour lines in steps of 10%. Fig. 4 shows that the pressure in the 0.6 [Ir/[9 model is altered by less than 0.003% by introducing ‘11. In the temperature range between logT = 4.3 and 5.2, it seems as if the stellar tracks follow the iso—difference curves in a significant fraction of the star. This, however, is the region just below the photosphere. and the large temperature This figure shows contours of the difference in pressure between including and ignoring the \lJ-term, in the sense logmlAP/PI. The dashed lines going from llLJllJJIlllllllllllllll I l l 1 llLllllllllIlll I’t. loglAlD/Pl ........... x(H*), 23(I-Ie”) — — — Stellar structure nnnnnnnnn nnnnnnnnn ......... ......... I I I I I I IWIWI I I I TWT I I I I I I I I I I I I I I I I I I I I I I I l I I I 4 5 e 7 L0g(T/1Kl) 59 and density gradients there makes this region very thin. For the 0.6 11/16, model, this peak in the pressure-differences is centered around a. relative radius of r/R = 0.987, and already at 0.975 and 0.995, difference has fallen by an order of magnitude from its maximum. The extent of the affected region in a star is indeed thin. and the location is below the super-adiabatic uppermost part of the convection zone and far enough from the photospheric transition zone, not to affect it. These transition regions are even more shallow, but nevertheless very important for the stellar structure, as they act as upper boundary conditions for the whole star (Trampedach 61‘ (ll. 2004(1). The thin region affected by \I/ does not have a corresponding importance and the implications for the star as a whole are insignificant. The inner 94% of the star (by radius—99.8% by mass), on the other hand, pressure ionizes perfectly by means of the occupation probabilities without any in— tervention by the lII—term. 3.3.5 Conclusion As can be seen from Fig. 4, the \Il-term barely affects any stellar models. By ex- trapolating the stellar model-tracks in the figure a little, it can be estimated that the pressure in a 0.5 [149 star is at most changed by 0.01% and for a 0.41119 star, about 0.03%. As explained in Sect. 3.3.4 the largest change introduced by the \IJ—term occurs in a very thin layer, in a region of the star which responds linearly to such changes, i.e., it hardly affects the remainder of the star. 60 The lIJ-term is not a cause for serious concern regarding the \l’alidity of the MHD EOS for stars down to masses of 0.4/149. 3.4 MHD 2000 Prompted by helioseismic analysis and detailed comparisons with the principal com- peting equation of state (EOS) project, I here present an updated MHD EOS in— cluding a number of new phenomena, as well as the previously published updates of improved electric micro-field distributions and relativistically degenerate electrons (Trampedach & Dappen 2004a). The main new features include higher order den— sity terms in the Coulomb interactions; quantum effects. including relativistic effects; and a new formulation of interactions that. involve particles with zero net—charge, abandoning the use of hard—sphere interactions. These changes expand the region of validity towards higher temperatures and densities, to easily encompass normal stellar cores. I also include di- and poly-atomic molecules (in addition to the already present H2 and If;I molecules), in order to include stellar atmospheres in the domain of the Mle EOS. 3.4.1 Introduction The equation of state (EOS) is an integral part of most astrophysical analyses, having the two r6les of supplying the thermodynamic state of the plasma, and of serving as a foundation for opacity calculations by providing ionization- and dissociation—balances and detailed populations of all electronic states. In the late 19808 two very successful EOS emerged, both being parts of efforts to improve on the opacities available to the astrophysical community: The Mihalas— Hummer-Dappen EOS (Hummer & Mihalas 1988; Mihalas el al. 1988; Dappen el al. 62 1988) as part ofthe international Opacity Project (OP) which is nicely summarized in the two volumes by Seaton (1995) and Berrington (1997). and the OPAL EOS (Rogers 1986; Rogers et al. 1996), pursued at the Lawrence Livermore National Laboratory, being the foundation for the OPAL opacities (Iglesias el al. 1987; Iglesias el al. 1992; Iglesias & Rogers 1991; Iglesias & Rogers 1996; Rogers & Iglesias 1992a). From detailed comparisons between these two, fundamentally di'lferent, approaches to the EOS (Dappen el al. 1990; Trampedach el al. 2004c) and comparisons with he— lioseismic inversions (Basu el' al. 1999), a number of problems with the MHD EOS have been identified. In particular, experiments with the partition functions (Gong et al. 2001a) suggested that the fluctuating micro-fields in the Mle EOS were too efficient at destroying highly excited states. This result was supported by direct com- parison between a number of formulations "for the micro—field distributions, including those used in the MHD and the OPAL EOS (Iglesias &: Rogers 1995). Recent work by Badnell & Seaton (2003) indicate, however, that the equivalent level-populations in the OPAL EOS are improbably high at high densities. Nayfonov et al. (1999) in- troduced an improved treatment of micro—field distributions, dubbed Q-MHD, which also includes the effects of electrons screening the charges of the ions (and vice verse). I have implemented my own version of Q—MHD, as detailed in Sect. 3.4.5. Helioseismic investigations by Elliot & Kosovichev (1998) showed a clear sign of the Solar radiative interior not only being slightly degenerate, but also slightly relativistic. This motivated the inclusion of the work by Gong el al. (2001b) on relativistically degenerate electrons, as well as a number of other relativistic effects, as detailed in Sect. 3.4.3. In the Sun, relativistic effects become significant in the 63 radiative zone, increasing towards the center. Another effect of about the same magnitude at the solar center is the exchange interaction between identical particles, which I include in Sect. 3.4.4.2. Trampedach el al. (2004c) recommended introducing quantum diffraction (see Sect. 3.4.4.1) and higher order terms in the Coulomb interactions (see Sect. 3.4.4.3) to improve the agreement between the MHD and the OPAL EOS. Both these effects were included in OPAL, but not in the MHD EOS. In Sect. 3.4.4.2 I present the contributions from the exchange interactions be- tween identical particles, and derive an expansion for the relativistic exchange integral in terms of the generalized Fermi—Dirac functions. In the original MHD EOS, perturbations by neutral particles (atoms) were i11— cluded by means of the hard—sphere model. In order for the problem to remain tractable this term was only treated up to first order in density, and second—order terms were included in an approximate way, as discussed by Trampedach & Dappen (2004b). The hard—core interaction model is inherently flawed, as it introduces a di- vergence for high densities and is undefined for densities higher than the close-packing density of the spheres. The interactions involved when neutral particles interact are obviously still caused by the charges of the particles constituting the atoms; in close encounters between atoms, the electronic wave-functions will overlap and the nuclei will therefore no longer be completely screened from each other, and net—forces arise. This phenomenon is described through the introduction of the concept of effective charge in Sect. 3.4.4.4. In Sect. 3.4.6 molecules (other than H; and H?) are introduced into the MHD 64 EOS, through parameterized partition functions. Inspired by the original MHD EOS. analytical approximations to the occupation probabilities of molecules are derived, based on the two molecules, Hg and 11;. that are treated in detail. 3.4.2 Chemical Compositions Comparisons between the original MHD EOS and the present updates will obviously depend on the chemical mixtures used in the calculations. and I therefore introduce the reader to my choices of chemical compositions, before proceeding to the equation of state issues. To facilitate comparisons with earlier work (6.9., Trampedach el (Ll. 2004c) and with the OPAL EOS, a H/He-mixture and a 6-element mixture are used, as listed in Tab. 1. The logarithmic the table lists abundances relative to hydrogen, [N,-/NH], which is normalized to 12. Mix 1 and 2 have helium mass—fractions of Y = 20% and 16%, respectively, and a. hydrogen mass-fraction of X : 80%. The right—most column in Tab. 1 is the ionization potentials for the first ionization stage, X11 to indicate the temperatures for which they start to add electrons to the plasma. The 6—element Mix 2 has a metallicity more than twice that of the Sun, Z9 : 1.8% and a distinctly sub—solar helium abundance, Y9 = 0.245 (Basu & Antia 1995), exaggerating metallicity effects with respect to most stellar applications. This is, of course, useful for analyzing the effects of metals on the EOS. Next, I need to address the issue of the best solar composition, i.e., metal— mixture, to use for the final table calculations. The computational cost increases 65 Table 1. Chemical mixtures in [N,/N1.1] (see text) element Mix 1 Mix 2 GN96 Optimized xl/[eV] H 12.00000 12.00000 —— 13.59844 He 10.79902 10.70211 —-— — 24.58741 C — 8.90307 8.55000 8.55884 11.26030 N _ 8.30307 7 97000 7.97000 14.53410 O # 9.23307 8.87000 8.87000 13.61806 Ne * 8.72077 8.08000 8 11100 21.56454 Na — — 6.32000 6.32895 5.13908 Mg W -— 7.58000 7.58000 7.64623 Si # A 7.56000 7.56000 8.15169 K h 7* 5.13000 5 13000 4.34066 Fe fl — 7.50000 7.65924 7.87052 0.4. . A I n L A l A . . . l 1 n . . I . A . . l n n A41 - truncated GN96 r - - - _ optimized mix 0.2 " | _ :i ' ‘. E '. <1 , i x4000 0.0- ------ , —0.2 Fig. 5. The relative difference in mean—electron-mass, he, as a function of tempera- ture in a solar stratification. The solid line is for the GN96 mixture truncated from 30 to 11 elements, and the dashed line is the optimized 11 element composition, as listed in Tab. 1. dramatically with the number of elements and the complexity (atomic number) of the elements considered, but the characteristics of the full mixture needs to be preserved. 4.0 4.5 5.0 long 5.5 66 I I I I I 1' I I I I 6.0 For the full mixture, I use the compilation by Grevesse et al. (1996) (GN96). As measure of the characteristics of a particular composition, I use the average charge Z = N. Za/ 1v... (4:: < where the sum is taken over all ionization stages of all metals. The population ( all ionization-stages of metals along a solar stratification, is estimated, using simpl Saha-equations (Saha 1921), with ground—level statistical—weights for the partitior functions. The approximations here should be immaterial, since only the location and relative strength of ionization features are of importance, here. The GN96-mixture is truncated, and the remaining elements are padded, to bee reproduce the (Z) of the full mixture, as well as the average mass of particles in th plasma. The optimization is performed on ln(Z) in order to capture the details c ionization in the parts of stellar atmospheres where hydrogen is still only a mine electron donor. A good compromise between speed and precision can be struck with a 11—elemer. mixture, comprised of the elements as listed in Tab. 1. The column labeled GN9 is the CN96—mixture truncated to 11 elements, and the last column is the optimize version. The optimized mixture has a RMS—deviation from the full GN96—mixture c 3.8%. As is evident from Tab. 1, most of the padding has been added to iron, som is added to neon, and the rest is shared between carbon and sodium. Potassiurr with its ionization potential of 4.341 eV, the lowest of all elements up to 37Rb, whic. is more than 300 times less abundant, is crucial for reproducing the details at 101 67 temperatures. Leaving it out results in over—all discrepancies of 22% and differen approachng 100% at low temperatures. Optimizing just the four metals of Mix 2 results in an over—all RMS—deviatf of 45%. Below 6300 K the deviation increases rapidly to 100%, i.e., this mixt‘ supplies none of the electrons donated by the low—ionization—potential elements 1 sodium and potassium of the 11 element mixture. As can be seen from Tab. 1, 1V 2 has no elements with ionization potential below 11.26 eV. From here on, the optimized nine metal—mixture will be referred to as the so (metal) composition. 3.4.3 Relativistic Electrons This section is concerned with the consequences of an appreciable fraction of t electrons having relativistic speeds, i.e., pe ,2 mec, where the masses of particles z always understood to be rest—masses. I let the treatment of relativistic electrc be dealt with prior to the other EOS improvements of the present analysis, as influences most of these. Due to their much larger masses, ions only become relativistic at temperatu: exceeding 1012 K. These temperatures are four magnitudes higher than considered the MHD EOS at present, and they can therefore be safely ignore for now. For an ensemble of particles of species oz, I use kBT mac2 68 as a measure of the importance of relativistic effects —— the ratio of kinetic- to r1 mass-energy. As mentioned above, only electrons are significantly relativistic, : the sub—script can be dropped, to simplify notation; 13 2 [3... The present sect is, however, kept general and applies to all fermions. Note that. there is a lacl- consensus in the literature on the definition of 6; sometimes the inverse is used. There are a number of effects arising from particles having relativistic spee First of all, special relativity complicates the relation between momentum, p, 2 kinetic energy, E, changing the general distribution function _. 871' p2dp 1(1). Eldl) = 5:3,- 1 from (p2 = 2mE) , 1/21 ftp. E1du = 1‘; t to (p2 = m2c2[(,3u -1— 1)2 — 1]) [3n +1 2 — 1 f(p, E)du : n" \/( 1+ 62:77 (flu +1)du . ( 71.. : 4—(2'm/cBT) : 69 The thermal de Broglie wavelength for particles of mass. m, is /\0 = h/e/27rkaT , (52) and the subscript indicates that the thermal average is evaluated at 7] = 0, 27.6., using Boltzmann statistics. The two distribution functions, Eqs. (48) and (49) have inspired the definition of the (non—relativistic) Fermi-Dirac functions 00 qu'u. F110) 2/0 W a (531 and the generalized Fermi-Dirac functions 00 u e/l + ifind-u. )2/0 2 for 1/ > —1 , (54) 1+ 611—7} 1/ (773/8 neither of which have analytical integrals. I solve the generalized Fermi—Dirac func- tions numerically, as described by Gong el (Ll. (2001), who also list all derivatives up to order three. For convenience, I will list the two first—order derivatives, as they will be needed in Sect. 3.4.4. The differentiation is straight forward, and gives _7__,,a(7,fl )_/000 U (N ‘1‘ 73515 (In. (55) 1+ 621—77 1+ 672—11 and _77,_(8F,,(fl)1/000 U V1 + $311 udu (56) :4 1 + eu—n 1+ gel... ' 70 From Eqs. (49) and (54), I see that the number density of relativistically degen- erate particles, is n : [000f(p,E)(lp (57) = n.[mememateriel . (58) which is needed for normalization of ensemble averaged kinetic quantities, as e.g. the 6- and O—factors discussed below. The largest effect of the change of distribution function, is a change of the trans- lational free-energy of electrons, F3. I follow the formulation of Gong el. al. (2001b) for F3 and its derivatives. The change ofdistribution function, also means that ensemble averaged momenta and energies will change. This has an effect on the degeneracy factors, 06 and O... modifying the electron screening contribution to the Debye—length, Eq. (61), and the ensemble averaged De Broglie wavelength, Eq. (68), respectively. These factors will be treated in detail below. Yet another effect due to relativistic speeds, is that of retarded potentials, ac— counting for the finite speed of propagation of electric fields. It can be shown (see, e.g., Landau & Lifshitz 1989, Chapter 8) that the potential produced by a moving charge, is (59) where R is the vector from the charge, to the point of observation. 71 3.4.4 The Bee Energy of the Coulomb Interactions In the original version of the MHD EOS (Hummer & I\="lihalas 1988; Mihalas el al 1988; Dappen et al. 1988), the Coulomb free energy, F11, is approximated by the Debye-Hiickel free energy ATBTI/ — -— A73 0 FDH 1271’ D (6 where AD is the Debye length 4aeg ~. . A72 : N96,, No.22 61 D [CB 71"; _, + 0;; 0' 1 . Since most expressions involving AD, do so in negative powers. I introduce the in- verse screening length. RD : A51. In the summation in Eq. (61). 0 runs over al particles with a net charge, except for the electrons, for which I also include effect: of degeneracy through the 0.3-term. The non-relativistic (NR) 0,, used in the original version of the MHD-EOS, C 6N7 = F—1/21711/2F1/‘2177l 7 (62: is based on the interaction propagators, G2, for the electrons (Cooper & DeWitt 1973) G; _ 87r V '00 p2dp 1 Ne _ 53 NC 0 6“"? + 1 1 + e"’“ i 0, = (63‘ I which can be generalized to the relativistic case by using Eq. (49). From Eq. (55) 1 72 now recognize the integral as (917,,(77, [fl/an, so that 6 : 8171/2191 51/871 71" (78173/20], mm?) C F1/21773l3)+l3F3/2(77,,13) ” (64) which through the recursion relation, (1/ —1— 1)F,,(77) : 8F,,+1/c)7]. yields the previous result in the non—relativistic limit. This new expression for 0... however, does not account for other relativistic effects in the plasma interactions due to, e.g.. the retar— dation of the potential. For lack of a better theory. I choose to limit ourselves to the relativistic effect on 66. In Fig. 6 I show the 06—factor in the non-relativistic limit (solid curve) and for 5 up to 1.0 in steps of 0.1 (dashed curves). Degeneracy of electrons strongly inhibits them from contributing to the screening of the ions, i.e., since they can hardly change state, because none are available to them, the electrons will act more like a uniform, charged background. Relativistic effects only make the electrons slightly more efficient at screening. In the following sections, I will present improvements to the Coulomb free en— ergy, based on quantum mechanical effects, higher-order density effects, and efleclive charges to account for extended particles. The final expression for F4 is summarized in Sect. 3.4.4.5. 3.4.4.1 Quantum Diffraction In the original MHD EOS, FD” was multiplied by a correction-factor to account for limits on the distance of closest approach between two particles. This correction 73 l l l l 1 l l 1 1 1L 1 J L l l [A I l I l l l A I I l I l l l l l l 1 l l l l I I I I I I I I I l I I I I I I I I I I I I I I I I I I I l I I I I I I I I I l -10 0 10 20 30 77 Fig. 6. The two degeneracy factors, (9,, Eq. (64), and 96 Eq. (78). The solid line shows the non—relativistic case, for which 0.217710) : Gem, 0). The dashed lines display the relativistic 6,, for [3 up to 1.0 in steps of 0.1, and the dash—dotted lines show the same for O... The behavior in )3 is monotonic in both cases. factor has since been disputed (rI‘rampedach et al. 2004c). and it was suggested to abandon it for the more physical phenomenon of quantum diffraction. Quantum diffraction removes the short range divergence of the Coulomb poten— tial by means of Heisenberg’s uncertainty principle. The wave nature of the particles results in finite charge densities at the position of particles, thereby leading to zero potential at r = 0. This effect results in a suppression of the Coulomb interactions at distances comparable to a deBroglie wavelength, which I accomplish with 74 1.2 - 1' 1" 7 I I l I 1-0: ,' Fit from this paper :' . i — — -- Riemann et al. (1995) . -1 2 .’ ...... _ ' ' 0.8 a 1, Emall 7. expansion _ . -. ........... arge 7, expansmn . P: 7 L x“; 0.6 . . l i 0.4‘ — 0.2— _ 0.0 ' L ....., ..., .... . . ....., . . ..-..., 0.01 0.10 1.00 10.00 100.00 1000.00 10000.00 7. Fig. 7. The quantum diffraction term, as evaluated from Eq. (65) and (66), compared with the Padé—formula given by Riemann el al. (1995) and with analytical expansions for small and large 7. where _ b .I' 2 (1’7 [1 -1'- w] (00) with coefficients a : 1.5393236 and l) : 02204-7583. (67) As is evident from Fig. 7, this form of T and 1', makes a good fit to Eq. (10) of Riemann el al. (1995), but without the singularity at 7 2 1570. From 7 ’2: 4 and up to the singularity, the two fits differ by up till 0.6 %. but my fit is much closer to the large y—expansion than to Riemann et al.’s fit in this region. I suspect that their fitting coefficients are merely published with too little precision. The asymptotic «.1 CI! behaviour of 7‘ for large, as well as for small 7', is the same as given by Riemann et al. (1995). The argument, 7, is the center of mass de Broglie wavelength in units of a Debye length. Defining the single particle de Broglie wavelength as A” = 73/1)... I can define the thermally averaged, relative or center-of—mass de Broglie wavelength between two particles, as .13, 2 13+ 13. = t2? [ + 1 52 O... 95 j — 213377 1m.O + 551 (67) where I use 1 for a 75 e (90. = (69) / for = e The general expression for (79—2) is 00 “2 ' E 1' (170—2) : A P ftp. )cp , (70) 16” f (7). E141) Where E is the energy of the particle, 11 its momentum and f(p, E) is its distribution function. For the general case of relativistic Fermi—Dirac-statistics _ 87r 792 _ 53 1 —1— e5—17 f0). E) . (71) the usual dimension—less energy 11 2 E/kT is introduced, resulting in the relativistic 76 expression for the momentum p : mcfifiu —1— 1)‘2 — 1 (72) dp :2 mc , . (73) Using \/(flu + 1)2 — 1 = «2611.01 + 1311/2, I can rewrite (1}), which is also the main part of the integrand of the numerator of Eq. (70), as r 1 1 szm.C(/g 73“ C” (74) 2(/1+,3u/2 fl. Realizing that this is the /3—clerivative of the generalized Fermi-Dirac function, Eq. (56, I obtain dp 8F_3/2 8111/2 ———-— : 2’ 2 3— . 7F 1 + eu-n m 7 1 8.6 +’ 8,6 1 3) For the denominator of Eq. (70), 7207.13), 1 similarly obtain PQdP 3 3/2 1 + 6%,, Z @0770) E [Fl/2 + flFB/Qj - (76) Combining Eqs. (75) and (76), finally results in 2 8F. 0 8E. (713 (7)—2 : , 3/2/ 5 ‘1‘ [3 1/2/ . . (77) mekBT F1/2 + 3173/2 In the non-relativistic limit the terms in [i can be left out. The classical limit results \I \J in (p52) 2 (mekBT)_1, and consequently 8113/2/06 + (3017—1/2/813 Oe : 2 171/2 + (3173/2 which accounts for the degeneracy of electrons. in the thermally averaged De Broglie wavelength. It is trivial to also account for Fermi-Dirac statistics of ions, but because of the mass difference I have chosen to ignore this, as also stated in Eq. (69). The behaviour of 6),, is shown in Fig. 6. I notice that in the non—relativistic limit, 68(7), 0) : Oe(n,0), which can easily be verified by setting 1’3 : 0 in Eq. (56) and compare with the non—relativistic 6., in Eq. (62). From Fig. 6 I see that the average cle Broglie wavelength of electrons decrease both with increasing 7) and increasing ,6. In both cases the momentuni-distribution will be peaked towards higher values and hence decrease the average de Broglie wave- length. For 6 = 1, quantrim-diffraction of protons will be equally important to that of electrons at 7),, 9: 51. At the ,8 ,3 0.017 explored in the present paper, this does not become a concern. even for an 77., which is orders of magnitudes larger than the value above. Dividing )10/3 with AD, I now get (79) _ r. e. g ”2 7”” _ AD./—_2chT ' m 0. mg In order to incorporate this one—component plasma (OCP) quantum diffraction term into my multi—component plasma, I need to write out KID in the expression for 78 the Debye-Hiickel free energy, Eq. (60), and relate T(70,3) to the pairwise interactions between particles of species a and 6. I rewrite F4 as 2\/—e3 F NZQIYO 03+ \/. .. 80 firs—Wye A): Wye-(A, {New (1 and note that the T’s have to be squared to agree with the one-component plasma. case for which 7' was developed. In order to save computing time, I define an average 7' that can be taken outside the summation over elements ,_ 2 :0 1V0 Zz\//Ve 0e T211706?) + Z/BN 32/37' 211/06) (81) Ea NO’ZEY \/N€6e + 213 A7213? where the two sums in the denominator are decoupled, as in the classical expression for the Debye—Huckel free energy, 3.6., Eq. (60). A very good fit to the above expression can be made from only two contributions rB/Vflt 171+(1—N117'12)» (82) where the first contribution is the combined electron—electron and electron—nucleon contributions and the second is the nucleon-nucleon term. The RMS deviations of the two—component fits are less than 0.25% throughout the compositional X/ Y—plane. The three new parameters, N1, ’yl/cyee and 72/7,... can furthermore be fitted to 79 binomials of order 5 in the X/ l”'—plane, 6.9.. l = Z(1,1Xil“"l (83) Tee ij '7—2 : Zb,,-\”’i'l’ (84) ”Yer: ij N1 = Zc,,-.\”l'1. (85) 20' I carried out this fittin ‘ )rocedure. usin ‘ relative Inetal—abumlaiIces from Grevesse 8 l 8 & Noels 1992), and obtained the fitting coefficients listed in tables ‘2-7—11. . RMSKT'm‘T')i l frl I I '0 C a o 0 G ¢ 0 '. ’. '. ’2‘ '. '. ‘. I I I l I I 0.0 0.2 1.0 Fig. 8. The RMS—deviations of the 7" evaluated from the fifth-order binomial fits, from the exact expression, Eq. (81). The dotted contours are spaced 0.0125% apart and solid contours are 0.05% apart. Using the 7’s and Nl’s calculated from Eq. (83—c) results in overall RMS devi- Table 2. Coefficients, (£0, for Eq. (83). ’Yl/‘Yee 0 1 2 3 4 5 0 9.1782e-01 —1.1232e-02 1.7127e-03 -4.2358e—04 6.0588e—05 —2.4894e—0! 1 —4.7200€-02 4.3322e—03 -2.6263e—04 8.9617e—05 —1.12360-05 4.3540e-0‘ 2 1.0418e—02 —8.7279e—04 1.6545e—04 —7.2533e—05 9.7177e-06 -4.1137e—0‘ 3 -1.9114e-03 1.5821e—04 -4.6280e—05 2.2475(3—05 -2.9216e—06 1.1581e-0' 4 2.1899e-04 -1.2414e-05 4.4749e-06 -2.4821e-06 3.0554e—07 —1.1137e-0: 5 —8.6693e—06 2.1068e-07 -1.4091e-07 9.1226e-08 —1.0532e—08 3.4806e—11 Table 3. Coefficients, ()0. for Eq. (84). 72/7... 0 1 2 3 4. 5 0 5.0211e—03 1.7181e—04 -7.1055e-05 3.0452e-05 —4.5304e-06 2.6472e-0' 1 3.5494e-04 —1.1401e-04 2.4850e—05 -6.1806e-06 9. 74-1 7e-07 —4.4853e-01 2 —2.3679e-04 1.2572e—04 —1.4812e-05 3.5951e—06 —6.3073e-07 3.2206e—0: 3 1.0054e—04 —4.41.03e—05 5.8336e—06 —1.5860e-06 2.7030e—07 —1.3612e—02 4 -1.5879e—05 6.2142e—06 —7.6992e-07 2.2073e-07 —3.6915e—08 1.7905e-0! 5 9.5203e—07 296300-07 3.1661e—08 —9.4868e-09 1.5661e—09 -7.3255e—1i Table 4. Coefficients, CU", for Eq. (85). N] 0 1 2 3 4 5 0 1.42136—01 1.2357e-02 —9.3669e—04 6.5955e—04 —8.8170e-05 5.2504e-06 1 2.7666e—02 7.9054e-04 5.8247e—04 —1.9683e-04 2.8651e—05 —1.5253e—06 2 —1.8022e-05 8.3095e—05 —9.2704e—05 3.3238e-05 -4.1539e—06 1.5133e—07 3 5.6812e-04 —1.6656e—05 6.4215e-06 —3.4200e-06 1.6857e—07 7.5673e-09 4 —6.1064e—05 —2.9772e—07 —2.8268e-07 1.2400e—07 2.4253e—08 —2.7578e——09 5 3.00186-06 —8.9259e-08 7.8630e-09 1.6651e-09 —1.9424e-09 1.4887e-10 81 ations of less than 0.25% and single T(’)’ee) deviations of less than 2%. The largest deviations, as seen in Fig. 8, occur at 766 2 102 around X = 0.7, Y : 0 — a rather rarely encountered region of the X-/Y—plane. For the moment I have neglected electIon—degeneracy effects in the fitting of 1;, but this can be remedied by using fl (")9 "/1 ~ —> ———,/—’ — 86 ’1 Am/‘ZkT me (wee) ( l ,1. l “/2 > .. —> ————,/ , 8f ,2 /\D V 257T 777%: (flee ( ) where the 7—1‘atios are the no-degeneracy fits, that can be fitted with the coefficients from Tab. 2 and 3 respectively. This whole procedure might seem a bit complicated, but bear in mind that the three parameters, NI, "fl/“fee and 72/766, need only be calculated once for each composition, and during the bulk of the calculations, the number of T—calculations will be reduced by an order of magnitude. 3.4.4.2 The Free Energy of the Exchange Interaction The free-energy of the first—order exchange interaction is Fex : 2: Na_—_:_ [1070) + j(7707fl0)l (88) (Kraeft et al. 1986), where the sum extends over all particle species. The thermal de Broglie wavelength, )‘cyO [See Eq. (52)], is based on Boltzmann statistics, and all 82 degeneracy effects are included in the exchange integrals, 1(770) and .7070, 60). I reformulate Eq. (88) by collecting all rye-dependent terms in one factor, with the relativistic effects included as a Taylor expansion [1262 . 22 ——— [V2 01 Jex ' as o - 89 16rr1«"'kBT g: “mo .(77 6 I I ) Fex: In this form, it is easier to recognize the dependence on the independent vari- ables, particle—mass and particle—densities, the latter of which indicates the binary- interaction nature of the exchange-term. The exchange integral, Jex, can be written 1(0) + .7071/3) Jex , 3 E , , , 90 (77 / ) [fir/207,3) + (3F3/2(7]a/3)l2 ( I where _ , _ 1(77) 'jeX(7790) — [MUD _ [Tl/2(7))? (91) is the non-relativistic part. DeWitt (1961) derived an expression for 19x, amenable to numerical integration 2 fl... 1’31/2(77’)d77’ F12/2l77) ’ [eX(77) which is shown as the upper solid line in Fig. 9. 83 I find that [eX can be approximated by the analytical expression — B 1 — 3"“? I... E 2x [CB + X8] 1/ X : —°_ (93) .47] using _ 4 A = 1.0129034 , B = 1.5384386 and C = 9—4 . (94) This fit agrees with a numerical integration, within 0.75 %. The difference between -10 O 10 20 3O 4O 50 77 Fig. 9. The exchange integral, ch(7}, 6), as function of degeneracy parameter 7]. The upper solid line shows the numerical integration of [6,.(77), indistinguishable from the fit of Eq. (93). The dot—dashed line (right hand axis) is the relative deviation between the numerical integration and the fit for [ex. The other lines show the behavior of Jex with increasing 6 (in steps of 0.1) going from top to bottom. The fit of Eq. 103 is shown with dashed lines, but only discernible above 77 2 30. my fit and the numerical integration is shown with the dashed line in Fig. 9, on 84 the right—hand-side axis. The fit is constructed so as to give the correct asymptotic behavior in the non—degenerate limit: lim [ex : 2. and in the completely degenerate 7}—>—c\3 limit: lip] Iex : 9/(277). This expression is an alternative to Eqs. (39) and (40) of 77 oo Kovetz et al. (1972). Following Kovetz (at al. (1972), 1 write the integrals as .1"?qu “3 .rgdx'g eul—n +1 0 6112—!) +1 [ 2 1 (9192 —1+J‘1J'2)] X — ——111 9192 1‘1.l/14l‘2'.1/2 yllJ'Z _ 1 _ 1'11‘2 where .r, : pi/mc, y,- = ,/1 + .17? and u,- : ei/chT are the dimension-less momenta, 1+.7 : W] (95) 0 energies (including rest-mass) and kinetic energies for the two interacting (identical) particles. 1 re—write in terms of u, and get in the general case m (1111 m (la-2 I z / __ __ . + j 0 6111—7) + 1 .0 6112—77 + 1 (96) >< [46,/u1u,-2\/l + %,t’3ztl\/1 + $621.2 — lnA where the argument of the logarithm is A _ 611.111,? + u] + U2 + 2,/u1u2\/1 +%/3u1\/1+%fiu2 _ ’7 ,. . (97) 6am? +111 +212 — 2,/'u,1u.2\/1+ %,13'll.1\/I + $6112 From Eq. (54), I recognize the first term in the square-bracket of Eq. (96) as 4F12/2(77, 6). 85 In the non—relativistic case I get I Z _ [“3 i; m ”(11me (9 0 61‘1“" +1 0 GHQ—77 +1 with u1+ 21,-; + 2 11.1212 ul + 212 -— 2,/u1'u.2 ANR : (9'. Since the non-relativistic part has been solved through Eqs. (91) and (92). I on seek the relativistic correction dul 0° (Ill-2 1 < A > (10' eul‘?7 + 1 0 GHQ—77 + 1 ANR .7 = 4flF12/2(n,fi) ~— / This integrand lends itself to expansions in the non—relativistic Fermi-Dirac i: tegrals, as carried out by Kovetz 61 (L1. (1972). Since my code already evaluates re ativistic Fermi-Dirac integrals, 1 instead perform. the expansion in terms of FV(77,6 Taylor—expanding A/ANR \/1 —|— 15621.1\/1 + flhl’z i in ul and U2. After collecting terms with the same combinations of powers of ul and U2, finally obtain mm =—— 3.6Ff/2 (112) _r/r 1 7. 2 7. : ZQ’Q. ‘l" AZQE 0‘ 72‘ <_) + —' +1 where AZQ. is the atomic number of particle a minus its net charge, 20,0. I then average this over the volume V and obtain the average effective charge. felt by ho— mogeneously distributed particles, Z 3 fig" 22( )1 (113) o. :2 —, 7‘ o, r (7* RE, 0 r- . 3 36—1/po 1 2 . 3 '4 2 20,0 + Add >< 60p; — (5 + 3,00, +10p0 + 20pa + 20100) (114) where pa, : ra/Ra and R3, : 3V/(47rNa). The term inside the bracket, makes the transition from 0 to 1 in the interval p 2 003—3. In the previous version of the MHD EOS, the interactions with neutral particles, were accounted for by means of a first—order approximation to hard sphere model. As this could be done using the occupation probabilities, it was possible to include the effect of highly excited atoms, having large radii. The second—order term was also included, although assuming all particles to be in their ground—state, through the \IJ—term, analyzed by Trampedach Sc Dappen (2004b). This rather ad-hoc term can now be abandoned. The hard—sphere model is an unphysical model, as the model is undefined for 92 high densities. Apart from wrecking havoc with numerical schemes for solving the equations, this is also unphysical, as the interactions at small distances should merely approach the Coulomb interactions between the bare charges of the nuclei. This is exactly what is described by using the effective charge, Z, from Eq. (112). 3.4.4.5 The New F4 1 summarize my changes to the Coulomb free energy, 1~1= 4.))me + F... (115) where g(/\) accounts for the higher—order terms in the Coulomb interactions, arising from the long-range nature of the Coulomb potential. The quantum diffraction-term, ?, accounts for the overlapping of particle wave—functions in close encounters, effec— tively removing the short-range divergence of the Coulomb potential. The exchange free energy, Fex, accounts for the interactions between wave—functions of identical particles. In addition to these changes, I have also Changed all occurrences of particle charges from their net—charge, Z0, to an effective charge, Z0, of particles of species a. This effective charge is the volume—averaged charge, assuming .5- wave-functions for all electrons bound in species oz, going from Za at low density, to the charge of the bare nucleus for inter-particle distances much smaller than the extent of the electron orbitals. This effective charge supplies small seed-charges, which are needed for the pres— 93 sure ionization to set in, via the occupation probabilities. Using just the fixed net— charge, Z, the plasma cannot pressure ionize at low temperatures and the original MHD EOS therefore included a hard—sphere—term in the occupation probabilities, plus an ad—hoc approximation to higher-order hard-sphere interactions, \I’. The lat- ter term was analyzed by Trampedach & Dappen (2004b) and found to be of very little consequence for stars more massive than about 0.4 MG»- 3.4.5 Improved Micro-field Distribution The occupation probabilities of the MHD EOS are based on the ionization of excited atoms by fluctuating electric fields, as described by Pillet et (Ll. (1984). The proccess is iterative and goes as follows: the Stark splitting of levels moves the electron up while the next higher level is moved down, crossing the occupied level and allowing the electron to cross over. The field then changes sign, moving the level with the electron up to meet the next, down—shifted level, and the process continues until the atom is ionized. The effeciency and range of this process clearly depends on the distribution of amplitudes of this fluctuating electric field. A higher probability for large fields will make it less likely for even low lying states to survive, and will make pressure- ionization very efficient and vice versa. The original MHD EOS uses a linear approximation to the Holtsmark micro- field distribution, which itself (Holtsmark 1924) is the result for an interaction free plasma, z'.e., the low density limit. Nayfonov et al. (1999) implemented analytical 94 fits to a higher order theory, and I, in turn, have implemented their work as part of the new MHD EOS. My implementation differs a little from that of Nayfonov et; al. (1999), in that it is not optimized for speed, but for flexibility. My modular approach allows for any formulations of the micro—field distribution to be employed in the occupation- probabilities. Following Hummer & Mihalas (1988), the probability of an electronic state i in particles of species 5 to survive despite fluctuating electric fields, is Flcsr 16,, =/ P(F)dF E Q(Ff,“) , (116) 0 where Ff; is the critical field—strength for ionization of this state, and P(F) is the micro-field distribution-function, and Q(F) is the accumulated micro—field distribution— function. The field—strength is usually expressed in units of the field-strength due to pre— turbers at the average ionic distance from the perturbed particle, F0, 27.6., ,8 : F/FO. In the MHD EOS the approximation 16 2 ”3 7T) , (117) F0 : (13M, (9N- is used, with Nion : 2095,, Na, and (10 being the Bohr radius. The critical field-strength for ionization of level i of particle—species s with 95 ionization—energy X“, is , 2 4 ’2/3 fl. : 16-2/3__1123x.-./e zs Zi3+1 Z NoCo‘ agée The quantum—correction factor, [(3, as given in Hummer S; Mihalas (1988), i 16 < n )2 n + 7/6 , ‘— —(—————‘— 11 > .3 If), 2 3 n+1 72.2+n+l/2 l n g 3 where n is the principal quantum-number of state '2'. (um (nm The original MHD formulation used a rough fit to the Holtsmark distribution —3/2 Qorig( r825) : 6_(315 whereas Nayfonov et al. (1999) introduced : f(fi137 ZS) (L) 1 +f((313.25.a) ’ QI/Bis) with f(6..,Z..a) = _ Cl and Q=aX,X=n+aW% 96 (mm (an (122) (we The coefficients were fitted to calculations based 011 the micro-field distribution func- tions by Hooper (1966; 1968). 3.4.6 Including More Molecules I wish to include more molecules in this EOS, in order to be able to use it as a foundation for calculations of atmospheric opacities. The current calculations of atmospheric opacities are often based on crude and old fashioned EOS, sometimes lacking thermodynamic consistency. Also, the bound—free metal opacities are often outdated, which I wish to remedy by merging the metal opacities from the Opacity Project (OP) and the improved calculations carried out as part of the Iron Project (IP) (Hummer et al. 1993) with those of modern molecular line databases (Parkinson 1992; .Iergensen 2003, and references therein). Molecules other than I12 and its ion will necessarily exist in much lower con- centrations, and will therefore have a smaller net effect on the thermodynamics. I therefore relax my precision requirements a little, and settle for parameterized parti— tion functions, instead of including a detailed treatment of all bound roto-vibrational and electronic states. The entropy associated with a partition function Q is rs = E + NchTan (125) 97 and therefore the Helmholtz free energy is F : E — TS Z N [E0 — ATBTIIIQ] .. (126) where E0 is the dissociation energy of the molecule with respect to the zero-point of the implicated elements. For CH2, for example, I have EOICH2) : DOICH2) + E0(Cl + 2EO(H) :- DQ(CIIQ) + DO—‘ where scattering, 0),, is included. For T ——> 0, on the other hand, the intensity weighted mean (15).] = [m 5:,\J,\(l)\ (141) 0 without scattering, reproduces the fluxes of the monochromatic solution (Mihalas 1978, Chapter 3.2). An ad-hoc bridging function is used to interpolate between the two cases, resulting in the standard opacity _3 T. —\ T. I”: b azc OT (H)*+(l—e 307 )1»: . (142) The >1< furtermore indicates that Eq. (14-0) is evaluated from the continuous opacity only, excluding lines, and —T ‘2 2: NA] .13) 6 91/ IDA] .7 —T), /2 E . 1 J * . [\Je ivy) J is Weighted towards optically thin wavelengths. The (pseudo) source functions and the interpolated opacity are then used for the radiative transfer in the 3D simulations. The angular integration to obtain the Inean intensity is evaluated by keeping the simulation box fixed and interpolating it 114 to a tilted grid, exploiting the periodic horizontal boundaries. Only the rectangular part of the box having standard optical depth 7' < .300 is used for the radiative transfer calculations, and this part is furthermore rescaled to optimally resolve the temperature structure. For each angle, the radiative transfer is thus solved Nb," times; typically Nbin : 4 is used, and N9 = 2 and NC, = 4 for the latitudinal and azimuthal angular resolution, respectively. There are a number of weak points in the method, as implemented. The bin- mernbership is determined from the 1D average structure of the simulation, which will most likely not correspond to the stratification experienced by any of the rays of light going through the simulation. Since radiative transfer is a highly non-linear problem, the heavy reliance on this reference stratification is troublesome. The choice of bridging function between the optical deep and the optically thin opacities is rather arbitrary. The right limits are of course ensured with this expres- sion, but the details of the transition are not. Instead of the actual opacities of each bin 22(2) .17,\Ju.‘,\1, the log-equidistant opac- ities, 55,, are used. This approximation was used to minimize the size of the table, in the early days of these sin'lulations, and can be abandoned now (Skartlien 2000). Another weakness of the implementation of the opacity binning, is the calculation 0f Rosseland mean—opacities. The bf— (bound-free) opacities are calculated for the Whole table, whereas the bb (bound—bound, or line-) contribution is only included in the 1D calibration stratification. The effect of lines is then extrapolated to the rest Of the table, assuming the same ratio between opacities with and without lines along 115 iSO-T contours. It turns out that points in the simulations that have the same optical depth, are highly correlated in g and T, and fall along a narrow line. There is no reason that the factor correcting for lines, should depend in any simple way on the optical depth, but the bridging function between optically thick and thin, should, and I simply employ the same extrapolation scheme for both. I abandon most of these. approxin‘iations in my SOS method, except for the reliance on a 1D calibration stratification. I show that. the SOS method depends Inore weakly on the calibration model, than does the opacity binning, thus resulting in a more accurate local radiation field in the simulations. This will be addressed in Sect.416. 4.1.3 Primer on Opacity Sampling Opacity sampling (OS) consists of performing monocl‘n'omatic radiative transfer for a large number of random wavelengths. It is a statistical method because of the random selection of wavelengths (in practice, equally spaced in. e.g., logA). Due to the extremely complicated behavior of stellar opacities, more than 104 wavelengths are needed before the procedure converges. For early type stars and metal-poor stars a larger number of wavelengths is needed, in order not to miss the rather few but iITlportant lines in their spectra. Full OS is therefore. prohibitively expensive for 3D hydrodynamical simulations. Conventional 1D stellar atmospheres have been modeled with opacity sampling by, e.g., Plez et al. (1992), Kurucz (1995) and Asplund el; al. (1997). 116 4.1.4 SOS My aim is to reproduce the N198 ~ 105 OS solution, with orders of magnitude smaller number of wavelengths, NEOS ~ 50. The straight-forward method. of going through all combinations of NASOS wavelengths, and finding the set with the smallest RMS deviation from the full solution, is rather prohibitive; It would require the computa- tion of (NSSVVASOS RMS-differences. A pro—conditioning of the problem, by dividing it into smaller sub-problems, can make it more tractable. This is done by dividing the spectrum into Nreg regions and each region into Nb,” bins, and then choose one wavelength that best represents the total of the bin. The Rl\.’IS fitting is then carried out only within a bin or between two bins, reducing the dimension of the problem to the order of lV/{DS/NEOS or (1 [IPS/N/(SOSV, with N808 : Nreg >< Nbin. My scheme for selecting the NEOS wavelengths is performed on the ID reference stratification, and proceeds as follows: 1) I divide the spectrum into 1)”?ng regions, by requiring each region to radiate the same amount of energy at 7' : l. 2) In each region, the member—wavelengths are grouped together, or binned, according to the position of the minimum of the monochromatic radiative heating, (1,\ = MUA — BA) 7 (144) where J,\ is the usual zeroth moment of the monochromatic intensity and B) and K.) are the monochromatic Planck function and opacity, respectively. Note that this ex- Pl‘essmn includes scattering despite the absence of the actual source function (Mihalas 117 1978). The location ofthe maximum cooling (minimum of qy) is calculated as the cen‘ of-mass of all cooling above Ty : 1 (i.e., Ty g l), and wavelengths with no cool above that, are counted as continuum bins, i.e., having a cooling peak around Ty e The limits of the bins are determined for each region, depending 011 the dis bution of loci of heating minima. With equidistant bins. some bins will invaria turn up empty, but I also want to make sure that even sparsely populated region: the T—scale are represented well. My solution is to fit Nb,” Gauss—functions, to the logarithm of the distribut function, substituting -2 for log(0), to make any non-zero regions of the distribut function stick out prominently against. the background. I also add in a small repuls potential between the Gaussians, in order to avoid degeneracy. The limits of the l: are evaluated as the mid-point between the centers of the Gaussians. The wh procedure has proven very robust in real applications. 3) I have now chosen the Nreg regions and the Nb," bins, and need to decide on criteria for choosing the NEOS wavelengths representing these bins. The key quant for the simulations, is the radiative heating/cooling, (1, which directly determines effect of the radiation field on the hydrodynamics. I also need to ensure a physicz Ineaningful connection to the equilibrium state and observables, in this case the ‘11 In principle the radiative flux can be determined from integration of the heating, 0 since the spatial deviation of (1,3315 from (18,8, will in general fluctuate wildly, a separ Constraint on the flux is important. Secondarily, I also want to be able to calcul COHSIStent radiative contributions to the pressure and internal energy of the ste 118 plasma. For this, the zeroth and second moments of the radiation field; J and Is respectively, are needed. I also find it useful to demand that the total Planck functio is close to the nominal B = 0T4/7r. A sixth quantity that could prove useful, is th Rosseland averaged opacity. This turns out to result automatically for 7' 2, I whe the five previous quantities have been fitted for, as is evident from Fig. 12. I consequently fit for f (145 { (1(7) MT) 3(7) «1(7) 1(0)} |min(q)la [[efl' lJ(102)IJ(102)711(102) l Where Hefl' : aTe4ff/(47r) and a is Stefan—Boltzmann’s constant. The normalizatior. in Eq. (145) are necessary in order for the different quantities to carry similar weigh in the fitting. The use of J, Ii" and B at 7' = 100 is somewhat arbitrary, but it seem to result in a fair weighting compared to H and q. 4) I define the weight of bin j of region i, as a)”: Z 20y (I46 /\ Eban and f‘z'fl-Uz‘j = Z fxwx, (147 AEbinU and I similarly define w,- and f; of region i. 5) For each region, I find the two bins with the most members and fit the Nbin — C“Sher bins individually, minimizing RMS(f,-j — fy), and using the weight, wij, on th * resulting representative of that bin, ,j. 119 6) The two remaining bins, jl and jg, of region i, are optimized jointly to rep1 duce the remainder of the region Nbin fi,12wz',12 = fl? — Z fijk’lt’z'jk 7 ”—‘1’ l€.f;1 ‘l‘ (1 _ €)f,-’:2lwi,12 (14 A123 where an,” : wil—twig and c = {0, 0.2, 0.4, 0.6, 0.8,1}. For each combination, {j1,j; I compute the RMS deviation for the six values of 5, find the interpolated minimt and the associated emin(j1,j2). The smallest of these interpolated minima. then (I termines jl and jg and is used with 6min to finally find .0712- 7) This ends the task of finding the Nb," wavelengths /\(i,j),.), lc : {1, . . . , Nb,“ reproducing region i, and the whole procedure is simply repeated for the other regior 8) In the end, I refine the fit by optimizing the weights, wij —> 103, by means a Xz-minimization performed on the logarithmic weights (to ensure positive defini weights) and enforcing conservation of the total weight, 2010:]. : XXV" wyl. This la . q q . . step does improve the agreement between fog and f”O", and I prefer to leave 1t 1 despite the weakening of the definition of the weights. 4.1.5 1D radiative transfer I have carried out the wavelength selection procedure, detailed in Sect. 4.1.4, on number of MARCS models (Gustafsson et al. 1975; Asplund el' al. 1997), as shown Fig. 12. The atmosphere parameters are listed at the top of the plot. I have plotti both the full 105 wavelength MARCS results and my corresponding results with ju 50 wavelengths, but the differences are only immediately visible in KROSS for whiclf I20 Teff 10g109 [Fe/H] 0.2 ”E 0.0 o "a? E -O.2‘ —O.2‘<‘ E E —o.4‘ —0.4 g «5 '6 g —O.6* 4.63: z —O.8‘ —O.8 —1.0- m ”U Q) .2 Tc E L: 0 Z :c “U :: EU 5 ed '0 Q) .L“. '5 E L4 0 Z 52’ *8 :23 .2 —4 -2 0 2 -4 -2 O 2 -4 —2 O 2 -4 -2 O 2 logloT loglo'r long logloT Fig. 12. Heating q, flux H, J, B, 311' and Rosseland opacity for four of the test Cases (horizontally). In the three upper panels, the dashed lines (right axis) are the differences in the sense (SOS — monochron'latic). In the lower panel the dashed line is the Rosseland opacity from my method. The horizontal dotted lines are the zero— points for the differences, and the vertical ones show the location of T = I. From left to right, the columns represent a red giant, a cool solar-type star, a metal—poor “Sun” and Procyon—a sub—giant. I21 did not fit. For the other quantities I also plotted the differences (dashed lines, right hand axis). The photospheric dip in the heating, (1mm. arises from the transition from convective to radiative transport of the flux. In the second panel from the bottom. 311" is larger than B, which is larger than J, only diverging above the photosphere. In this panel I only plot the J—differences. which are representative of the B- and K—differences too. Fig. 12 shows that my method is successful in reproducing the full radiation field to within a percent. The metal—poor atmosphere (middle—right) is easiest reproduced with few wavelengths, whereas the red giant. (left) causes more problems due to the many molecular lines in the relatively cool atmosphere. Notice that this situation is opposite that for pure opacity sampling (see Sect. 4.1.3). The reason being that with SOS, I have the freedom to choose my \lvavelengths to cover the few but prominent lines, whereas with OS the wavelengths have to be chosen randomly. The SOS method is obviously not applicable to conventional ID stellar atmosphere models, as it relies on a full OS evaluation of the radiative transfer in the reference stratification, as detailed in Sect. 4.1.4. 4.1 . 6 Spanning convection In Fig. 13 I test the old and the new method on the heating in a vertical slice of a solar simulation snapshot (Stein 85 Nordlund I989; Stein & Nordlund I998; Asplund et al. 1998). I plot deviations from the “monochromatic” calculation, averaged over the horizontal dimension of the slice. Solid and dashed lines denote plain- and RMS— averages respectively, the old opacity binning is shown with thin lines, and the new 122 0.10 ' ‘ ' [.1 III II‘ L [I ' -—-— (SOS - monochromatic) f, /: . I - — - - - RMS(SOS - monochromatic) ,\l 1‘ I], . l / 1 . — (bin - monochromatic) "H l/ l I — 0'05 - — - - — RMS(bin - monochromatic) '1“: , -6 -4 -2 0 2 long Fig. 13. Comparison of the heating from the opacity binning (bin, thin lines) and my new method (SOS, thick lines), with the “monochromatic” case (see text). A <. . .) denotes horizontal averages over a 1.5><6 Mm, 82X200 point vertical slice from a Solar simulation (solid). RMS denotes the corresponding horizontal root-mean-square averages (dashed lines). SOS method is shown with thick lines. The reference (“monochromatic”) model in this case is based on the ODFs used in the ATLAS atn‘1osphere-models of Kurucz (1992c), with 230 ODFs of 12 points, corresponding to 2760 “wavelength”-points. This is a far smaller set of wavelengths to choose from, and the fit for the ID average of the slice is about 10 times worse than what is shown in Fig. 12 for the MARCS models (Note that this does not say anything about the accuracy of MARCS or ATLAS models). In the future, the plan is to base my wavelength selection on an OS calculation for the 1D average stratification of the simulation with 105 wavelengths, 123 but for now, I can test my new method on the ODF case. The accuracy should increase with more choices of wavelengths. Fig. 13 shows that my wavelength selection is more stable against convective fluctuations than the opacity binning method—both the average— and the RMS— deviations over the slice are smaller with my new method. Some of that is due to my new method explicitly fitting for q(T) of the average structure, whereas the opacity binning is a forward calculation, with no feed—back mechanisms (i.e., no iterative fitting scheme). The transition from convective to radiative transfer of the flux, is accompanied by a. radiative cooling—dip. Fig. 12 displays the difference between the monochromatic and the binned case, exhibiting both a positive and a negative “bump” around logT ~ 1. This is the tell-tale sign that the cooling-dip in the binned solution is located (on average) a little higher in the atmosphere, than is the case in the full OS solution. Such a feature is not present in the SOS solution. 4. 1.7 Conclusion This first feasibility analysis has shown that I. It is indeed possible to find 50 wavelengths that can represent the radiative transfer of a full OS calculation with 105 wavelengths. 2. Both the heating and the first three moments of the intensity J, H, K are well reproduced. 3. Although the Rosseland opacity does not enter into my fitting criteria, it is 124 reproduced to within 1% in the optical deep layers. 4. The proposed selective opacity sampling method is more stable against the convective fluctuations than the opacity binning. I have yet to implement the new method in the 3D convection simulations, but these first tests are very promising. I25 Chapter 5 Improvements to stellar structure models, based on 3D convection simulations As a part of the grander scheme, to better understand stars, improved models of stellar structure and evolution are obviously key ingredients. In chapter 3, I presented an analysis of one of the leading equation of state projects, as well as some improvements aimed at broadening the range of applicability and increasing the accuracy of the MHD equation of state. This work is equally aimed at the convection simulations and the general stellar structure and evolution problem. These equation of state improvements are to be followed up by new calculations of opacities, as outlined in Sect. 3.5, which likewise has implications for both the convection simulations and stellar structure models. 126 With the progress in the evaluation of atomic physics, by far the most un— certain aspects of stellar structure models, are those involving dynamical processes which typically lead to mixing. Convection is the most important of the dynamical processes, as it not only is the most efficient mixing mechanism, but also directly affects the stratification of the star. This has been realized for almost a century, but the preferred formulation of the problem is still the rather simplified picture of the mixing—length formulation (MLT) (Bohm—Vitense I953; Bohm—Vitense 1958). In Sect. 5.2 I explore some of the reasons why MLT, after all, does a pretty good job at describing convection, and I point out some of the differences between MLT and what we have learned from 3D simulations of convection. Since the MLT formula- tion has a (principal) free parameter, a, a first. step towards a better description of convection in stellar models, is a calibration of 0 against the convection simulations, as undertaken in Sect. 5.2. In order to be able to perform this calibration of a, m5erything other than the ac- tual convection, has to be the same in the stellar structure models and the convection simulations, including the equation of state, opacities and the atmospheric boundary condition. The problem of stellar atmospheres is computationally very expensive, and cannot be incorporated directly into stellar structure calculations. Instead, the results from atmosphere calculations can be used in the form of T—T relations; tem- perature as function of optical depth, as confirmed in Sect. 5.1. T-T relations are derived from the simulations and applied to the stellar structure models, to facilitate the calibration of a in Sect. 5.2. The T-T relations from the simulations have also been fitted in the atmosphere parameters so that they can readily be used in stellar I27 modeling. One of the short-comings of conventional stellar atmosphere models is their re— striction to ID space, with. some dynamical effects included after the fact, e.g., micro— and macro—turbulent velocities. These atmosphere models also use the MLT for de- scribing convection, but as the largest deviations between the MLT picture and the 3D simulations occur in the atmosphere, the structure of an atn‘iosphere that incor— porates MLT-type convection, is unlikely to resemble a real stellar atmosphere in any detail. Using the combination of the oz—fitting and the T-T relations in stellar structure and evolution calculations, carries the promise of more reliable stellar models, at least in the solar neighborhood of the HR—diagram. The implications for stellar evolution will be studied in a future paper (Trampedach el al. 2004b). 128 5.1 T-T relations from convection simulations T-T relations are normally used for describing the photospheric transition from op— tically thick to optically thin in stellar models. This is well justified, but the impor— tance of the T-T relation as an upper boundary condition for stellar envelops seems not to have been fully appreciated. I assess the effects of employing often used as- sumptions about T—T relations on stellar models and illustrate the interplay between atmospheric stratification and the depth of an outer convection zone. Convection in the stellar structure models is described by the mixing-length theory (MLT). I present calculations of T—T relations based on 3D radiation—conpled hydrodynamical (RI-ID) simulations and give simple fits to my results for easy use in stellar structure calculations. 5.1.1 Introduction The calculation of stellar atmospheres is so complex that it has formed its own sub-discipline. Most complications arise because radiative transfer in the transition from optically thick to optically thin is hard to treat in a simplified manner without losing essential features. To treat this region properly, the radiative transfer has to be solved for hundreds of thousands of wavelength points. This obviously renders atmosphere calculations time consuming and impractical to incorporate directly in stellar evolution codes. A solution to this problem is to use the results of stellar atmosphere modeling (semi—empirical or fully theoretical) as upper boundary conditions for stellar struc— 129 ture models. Since T-T relations can be derived from limb-darkening observations, the use of semi-empirical models based on such observations often has been consid— ered the safest choice. Knowing pressure and opacity as functions of Q and T, and assuming hydrostatic equilibrium, the system of equations can be closed by the T-T relation without having to solve the frequency dependent radiative transfer — i.e., the stratification of the detailed atmosphere calculation can be recovered with a grey opacity, as described in Sect. 5.1.2. I proceed by giving a short overview of the 1D structure calculations in Sect. 5.1.3, where I also elaborate on the implementation of T—T relations. Theoretical T—T relations from 1D stellar atmosphere models have been published in connection with, e.g., the ATLAS (Kurucz 1992c; Kurucz 1995), the MARCS (Gustafsson et al. 1975; Asplund et al. 1997) and the NEXTGEN (Hauschildt et al. I999a; Hauschildt et al. I999b) grids of stellar atmospheres. These are grids in effective temperature, surface gravity and metallicity and dense enough that simple interpolation is safe. The level of sophistication is very impressive, the only weak point left, being the treatment of convection. In late type stars the modeling of photospheres is further complicated by convec— tion. Not only are the atmospheres no longer one—dimensional, but the fluctuations are also well outside the regime of linear perturbations. So far, the only way to deal with the combined problem of radiation—coupled convection in stellar photospheres, is to perform realistic radiation l‘1ydrodynamic (RHD) simulations. I describe such 3D convection simulations in Sect. 5.1.4 and go through the various ways of extracting the average structure from the simulations in Sect. 5.1.5. In Sect. 5.1.6, I compare the 130 results of my “high-precision” solar simulation with T-T relations from the literature. In Sect. 5.1.7, I present a single fitting formula for T(T), which works well for all the simulations. I furthermore give a list of fits, linear in logTeH and loggsurf, of all the coefficients, for use in stellar structure calculations. I apply the T-T relations, derived from the simulations, to stellar structure calculations in Sect. 5.1.9. I compare the effects of using some of the most common combinations of T—T relations and opacities, and point out some of the often encountered inconsistencies. I also explore how changes to the physics in the atmosphere changes the depth of the outer convection zone. This paper is the first in a series dedicated to the improvement of stellar structure and evolution calculations. These improvements are based on lessons learned from 3D radiation—coupled hydrodynamical simulations of convection in the atmospheres of a handful of solar-like stars. Here I present results on the radiative part of the mean— stratification of the simulations in the form of T-T relations. Paper II (Trampedach el' al. 2004a) deals with the convective part of the mean—stratification by calibrating the mixing-length and presenting results easy to implement in stellar structure codes. The radiative and the convective parts of the problem are strongly interdependent, as discussed in the two papers, but I also show that separating the two parts is possible and has a significant effect on stellar structure models. Paper III (Trampedach el‘ al. 2004b) will address the consequences of applying the above improvements to stellar evolution calculations. The effect of T—T relations on stellar evolution, has also been studied by Chabrier 85 Baraffe (1997) for low-mass stars with solar composition. They argue that the use 131 of T-T relations imply grey radiative transfer in the atsmophere. From my analysis in Sect. 5.1.2.2 I show that this is not the case; a T—T relation can fully describe a non- grey atmosphere, but the zeroth, Eq. 153, and first, Eq. I54, moments of the transfer equation are modified. These modifications imply the use of Rosseland opacities, also being the natural choice for the stellar modeler. Chabrier 85 Baraffe are rightfully concerned about the proper implementation and interpretation of the T-T relation in the transition between radiative and convective zones—an issue which is often overlooked. I hope that my discussion in Sect. 5.1.2 below will clarify matters, and justify the use of T—T relations, even in this region. Ludwig et al. (1999) have performed a similar calibration of the mixing—length based on 2D simulations of convection, , but employing a completely independent method. Their T—T relations are computed and implemented in a way similar to what I present here. Paper II provides further comparisons between the two methods and the results. 5.1.2 The Basis for T-T relations I describe the 1D, plane-parallel, radiative transfer in terms of the usual moments of the radiation field I It”) : %/_lu”1,\(u)du , (149) Where the intensity, Iy, only depends on the angle with the surface normal, ,a = cos 0. Dependence on optical depth, T, is implied throughout this section. Extension to the 3D case is dealt with in Sect. 5.1.5. 132 The first three moments are also called Fra r ‘ JA 2 [(0)7 47rd = Hy : [1(1) and [\y = [1(2) . (150) The transfer equation is d], , 11 C122”) : Iy(,u.) —— .Sy , (151) where the source—function, Sy, is isotropic. The corresponding radiative heating (cooling when negative) is Qrad,.\ = 47FQH,\(JA — SA) a (152) where the 47r comes from the angular integration of an isotropic quantity and Qrad is the extensive version (per volume) and qmd : de/Q is the intensive quantity. A solution in radiative equilibrium obviously obeys de : f6” Qrad,,ycl)\ 2: 0. 5.1.2.1 Grey Radiative Transfer In a grey atmosphere, the opacity is independent of wavelength, and all the A— subscripts can be dropped. Integrating the transfer equation over angle then gives Elf—{— __ idFrad dT _ 47r dT =J—s, (an 133 whereas the first moment of the transfer equation gives d1" Frad 0T4fl ( Ifconv) —— : ['1 Z : —€_ 1 _ 1'r4 dT 47r 47r Ftot ( J ) 01' T4» T Fconv ’ ; . .. (._ / 4.2.1.) (.55) "7f 0 F101 where the integral contains the convective effect. The convective flux is the sum of the enthalpy and the kinetic energy fluxes. It turns out that only AI\'(T) : I\'(T) — [{(0) is needed, and the value of the constant, l\'(0), from the integration of Eq. (154) is immaterial. Assuming local thermodynamic equilibrium (LTE), as I do in both the simula— tions and the stellar structure calculations, S = B : 0T4/7T. l\-"Iultij‘)lying both sides of Eq. (155) by 47rS'/(3AI\'UT:H) therefore results in the T—T relation 4 T 4 5' T E‘onv(7-I) — —- = — +1 ’ . 175 3 (TN) BAK IT ./0 Fm. ”l ( 3 l The T/AK—factor is convergent, since Alf increases approximately linearly in T and I can therefore describe finite temperatures. Using If instead of AH, would have introduced an extra term for the temperature at T = 0. All subsequent references to the T-T relation will only deal with the radiative part, 4 §(T) =,S_T, an) 134 unless otherwise noted. 5.1.2.2 Average Radiative Transfer If, on the other hand, the opacity does depend 011 wavelength, integrating Eq. (154) over wavelength gives oo 1 11' F”, f f ”(n = H = 4 d, (158) 0 iny dz 7r where I substituted dTy : QRyClZ. Fwd and H are just the results of direct integration over wavelength. As Jy ——> Sy and Ky —> %Jy for Ty —> 00 and .S'y : By in LTE, Eq. (158) leads to the usual definition of the Rosseland opacity ~00 1 dB 1 — ——————° 57TH” (159) . _ 00 dB ' ' h‘Ross 0 TTTTdA The differentiation with respect to z and T can be freely interchanged, as they are monotonic functions of each other.) From this, I can now define dli' 00 I dK Fm E mm] — A (n = H = J, (160) CITROSS 0 Icy CITROSS 47r where CITROSS 2 QKROSSCIZ defines the Rosseland oplical depth. The tilde refers to quantities that are averaged over wavelength in a way that makes them obey the grey transfer equations, Eqs. (153) and (154). In a similar way, the zeroth moment 135 of the transfer equation gives M E moss/mi dH" (0 = J— s. (161) O dTROSS I; ,\ C17-13055 where I again substitute 3 = B. In the non—grey case, the tight link between the zeroth and first moments of the averaged radiative transfer is broken, as H 7E H in general. The expression for the T-T relation, Eq. (156), is, however, unchanged when sul‘)stituting If and T3035. With the radiative temperature being given by Eq. (I57), and the temperature from the simulations by Eq. (156), I can reduce the temperatures from the simulations to the purely radiative ones by . F... .1 ‘1“ ITrad : TTl/4 [T —/0 #dTl] . (162) This is the stratification in the case of no convective flux. Deriving TradfT) from the simulations does, however, retain secondary convective effects, ie, the radiative equi- librium, as influenced by convection. As mentioned earlier, Trad is the temperature to be used in T—T relations. 5.1.3 The 1D-envelopes The T—T experiments have been performed on models of stellar envelopes (Christensen- Dalsgaard 85 Frandsen 1983), that each cover the range from a relative radius of T/R = 0.05 and out to an optical depth of T : 10‘4. I36 All time dependent and composition altering processes, e.g., nuclear reactions, diffusion and settling of helium. and metals, have been left out. This renders the envelopes functions of the atmospheric parameters. Teff and gsurf (and composition) only, but it also rules out any abundance gradients. Helium and metals slowly drop— ping out of the convection zone build up an abundance gradient, just. below the convection zone, which is counteracted by diffusion. As diffusion and settling only generate abundance gradients below the convection zone, my results on the T-T re- lation and the a—calibration should be independent of these processes, whereas the depth of the convection zone might be affected slightly. Radiative levitation of high- opacity elements, on the other hand, would have the largest effect in the photosphere and the change in composition would alter the opacity and equation of state so as to change the efficiency of convection. The 3D simulations, however, show that the convective overshoot sustains velocity-fields. at least out to T 2 10‘4, that would immediately wipe—out any chemical gradients in the atmospheres of the stars I have explored here. For the envelopes, I have used the same equation of state (EOS) chemical compo- sition, and, for T S, 104 K, the same opacities as for the simulations (cf. Sect. 5.1.4). For higher temperatures, I used the OPAL opacities (Rogers 85 Iglesias I992a). The difference between the two opacities is generally small at this temperature, and the transition is smoothed and always takes place in the adiabatic part of the convection zone, minimizing the impact on the structure of the model. Convection is treated using the standard MLT as described in Bohm—Vitense I37 (1958), using the standard mixing length I Osz (163) and form factors 1 = :1, and 1/ = 8. The T-T relations are supplied through Eq. (169) (see Sect. 5.1.7), and different choices for the coefficients then constitute the various cases listed in Tab. 7 below. 5.1.3.1 Implementation of the T-T relation The Hopf function, q(T), introduced by Henyey el’ al. (1965), is part of the T-T relation 4 T 4 § (100 for T ——> 00, and recovers the diffusion approximation. In deriving the T—T relation in Sect. 5.1.2 I make the transformation from actual to radiative T-T relation by changing T, as shown by Eq. (164). In the envelope calculations I need to re—introduce convection in the T—T relation, this time described by MLT. This is done with the inverse transformation, which is accomplished through a modification of the optical depth dT : f,,dT , (165) 138 so that (1(Tl‘l‘ T 2 (1(7) + T + (Iconval - (166) Differentiating both sides 0qu. (166) with respect to T, I get (1’0) + 1 + 61$...(T) q’tt) +1 ’ (167) fH: where primes indicate differentiation with respect to the argument and T is found from solving Eq. (166). Employing hydrostatic equilibrium and, as usual, defining the radiative gradient as the gradient that would be caused by radiative transport of energy alone, i.e., assuming that ‘T is given by q(T) + T, I find dlnP Hp 3Ftot V.., E = — 1 ’ . 1 ' d (alnT)md g 160T4( T (1) ( 68) The actual gradient, V, can similarly be found by using Eq. (166) and the transfor— mation to T, resulting in f,, = V/de. With these relations the T-T relation can be used throughout the stellar envelope model, without the (common) artificial transition between atmosphere and interior. 5.1.4 The 3D-simulations The fully compressible, transmitting boundary, RHD simulations are described by Nordlund 85 Stein (1990; 1989). Since the matching to ID envelopes demands a high degree of consistency between the simulations and the envelopes, I found it necessary 139 to bring the micro—physics up to the same level as that used for the envelopes. Direct comparison to observations of the Sun also necessitated an update. I therefore revised most of the opacity sources and added a few more sources, as described in detail in Trampedach (1997). The line opacity is supplied by the opacity distribution functions (ODFs) of Kurucz (1992a; 1992b), and the EOS is changed to the Mihalas-Hummer- Dappen (MHD) EOS (Hummer 85 Mihalas 1988; Dappen el al. 1988). Table 7. Parameters and derived convection zone depths for the seven simulations. name Star A 02 Cen B Sun a Cen A Star B 77 B00 Procyon Spectral class M 5 IV K 1 V G 2 V G 2 V F 8 V G 0 IV F 5 IV-V Teff 4851 K 5362 K 5801 K 5768 K 6167 K 6023 K 6470 K loglo gsurf 4.095 4.557 4.438 4.295 4.035 3.753 4.035 [VI/[149 0.600 0.900 1.000 1.085 1.240 1.630 1.750 a 1.8705 1.8313 1.8171 1.8032 1.7360 1.7383 1.7193 (lcz 0.5600 0.3063 0.2861 0.3070 0.1966 0.2087 0.1035 dcz solar 0.5616 0.3085 0.2861 0.3057 0.1870 0.1927 0.0925 (lcz 5000 A 0.5505 0.2971 0.2647 0.2795 0.1552 0.1618 0.0639 ([02 HSRA 0.5498 0.3016 0.2705 0.2866 0.1579 0.1623 0.0692 I have performed simulations for seven sets of atmospheric parameters, five of which correspond roughly to actual. stars, as listed in Tab. 7. The Star A-simulation was added to get a better coverage in the Tea/gsurf—plane and Star B was a. simulation of Procyon that turned out too cool (they therefore have the same gsurf). Mostly for historical reasons, I used a hydrogen fraction by mass of X : 70.2960 % and a metal fraction by mass, Z 2 1.78785 %. Each of the simulations was performed on a 50 X 50 X 82—point grid and covers about 4—6 granules horizontally and 13 pressure—scale—heights vertically, with 20% being above (T) : Tea. After relaxation to a quasi-stationary state, I calculated mean models as described in Sect. 5.1.5. 140 5.1.5 Averaging procedures I evaluated T-T relations for both the Rosseland optical depth, TROSS, and the 5000 A monochromatic optical depth, T5000, each with three different averaging procedures. The convective fluctuations in density and especially temperature are so large that the opacity and the EOS (e.g., gas pressure, Pg) are non—linear in the fluctu— ations, which has the consequence that (Pg) 74 Pg((g), (T)). This means that the average gas pressure is, in general, not related to the average density and temper— ature through the EOS. The relative difference amounts to about 5% in the solar photosphere. This effect is much larger for the opacity where the relative difference reaches more than 90%. This has led me to evaluate the three averaging procedures compared in Fig. 14, and detailed below. Case a): Calculate optical depth—averaged models by interpolating temperature, T, onto iso-T surfaces, equidistant in log10 T, and averaging over these surfaces. These T-averages were then subjected to temporal averaging and I denote this procedure (. . .)T. The radial p—modes, which are excited in the simulations, are to a large extent filtered out by this method because of the dependence of opacity on density and temperature. With increased temperature and density, the opacity and optical depth increase, moving the T—scale outwards with respect to the z-scale, in much the same way as the column mass scale. The temporal RMS-fluctuations of gas-pressure on a T-scale are about 10 times smaller than on a fixed height—scale. This effect is also clearly seen in the results of Georgobiani el al. (2003a), where the p—modes excited 141 Sun: Te": 5770 and g,u,,=2.74e+04 LIJIIJIIIIIIIIIIIIIILILIIIIIIIII11111111111111lllIlIIIllllII' 9000: , : 8000 _ 3) (T(Tn...» /’ _— : , -' ' :1 """"" b) () / .- 7000{ F E \ _ R - 6000{ E 3 5000: r .1 3,,. 4000? IIIIIIWITTIITIIIIIIIIIIIIIII[IIIIIIIIIIIIIITIITITTITT—IT—IIT —4 —3 —2 -1 0 1 2 Logtr) Fig. 14. The effect of various averaging techniques applied to a solar simulation (see text for details). Notice how cases a) and b) follow closely in the atmosphere down to (T) : Terr. at long 2 —0.2. in a solar simulation are very prominent in temperature, sampled at a fixed height, but almost vanish when sampled on the undulated iso—T surface. This averaging procedure is motivated by the form of the radiative transfer equation, Eq. (151). Since T is the only quantity entering the radiative transfer equation in a non—linear way (a division), this is also the only quantity that cannot merely be replaced by its horizontal average. I therefore recommend this averaging method for use in conventional 1D stellar models, and it is the method used for the rest of my analysis. This averaging procedure corresponds to observations of the 142 “radial” or disk-center T(T) in the sense that it is the average of T(T) along radial rays. Also, looking at disk center, it is not possible to “observe” the height—scale— only the T-scale. Case b): Calculate geometrically averaged models by mapping the horizontal averages onto a column mass scale and performing the temporal averaging on this scale instead of on a direct spatial scale. This approach filters out the main effects of the pmodes excited in the simulations. I refer to this procedure as Lagrangian averaging, (. . .)L. This Lagrangian averaged T as a function of the Lagrangian averaged T corre- sponds in a sense to limb observations, observing the Sun “horizontally”, instead of “radially” on an optical depth-scale. Real limb observations would also contain a case a) component due to the Sun’s sphericity. As noted above, this method is not compatible with the radiative transfer equation. Case c) Calculate (T)L as function of a T, based on integration of n((T)L, (Q)L) # (a). The non—linearity of the opacity is taken into account by this method. If Pturb and the non—linear effects on Pg were known, the correct g,T, P-structure could be recovered with such a T-T relation. This method mixes the convective and. optical parts of the problem, though, making it more difficult to improve their treatment in stellar structure calculations in a consistent way. Fig. 14 shows that the T—T relation can be derived fairly unambiguously from limb-darkening observations, for (T) g Teff. There the inhomogeneities are small enough for the opacity to be linear on the scale of the fluctuations, rendering the three averaging procedures equivalent. At greater depth this is no longer the case as 143 the inefficient convection at the top of the convection zone causes large and non—linear fluctuations, splitting case a) and b) apart. The T—T relation still has a significant effect on the envelope model at these depths, though, requiring a specific choice to be made. Based on Sect. 5.1.2.2 I choose case a), as also advocated by Ludwig (:1; al. (1999). 5.1.6 The solar T-T relation Apart from the simulations of the seven stars presented in Tab. 7, I also made a sim— ulation for direct comparison with solar observations. This simulation has a slightly different composition: YES. : 0.245 in accordance with helioseismology (Basu 8 Antia 1995), and ZQ/XQ = 0.0245 in agreement with meteoritic and solar photospheric metal to hydrogen ratios (Crevesse 85 Noels 1992). This ratio results in the hydrogen mass-fraction, X = 73.6945% and the helium-hydrogen number ratio, He/H:0.0837, instead of the often assumed lie/H: 0.1. The latter ratio is what I, for historical reasons, used for the seven other simulations, listed in Tab. 7. The Z / X ratio is reduced by about 4% from the other simulations, primarily decreasing the line—blocking by metals, resulting in a decreased (1. Helium, being an inert element, contributes little to the opacity. Its own opacity is very small in the solar atmosphere, and the high ionization potential means that no electrons are donated for the formation of H‘—ions, which is the most important source of continuum opacity in the Sun. The larger X (lower Y) therefore leads to increased opacity and the lower Z leads to less line—blocking. These changes to the composition I44 also alter the mean molecular weight, ,a, which enters both in Eq. (178) for the convective flux and in the hydrostatic equilibrium through its effect on the pressure. Sun: Te”: 5770 and 9“,,22.74e+04 600. * : a) l 400 ‘1 — Tan—sim_ THoI-eger-uuner :- : - — — T30-.1m— TAUOIQ - 200 " - — - - 3D-sim_ Tan-aim(Tsooo) L Io ~ ~ —. \ l — 0“ ~ 0 I IIIOOOODIDI “In-ICOOCIUI‘DIOID o O o— “wi':'£.¢=. f' u- \ '0'. INS - ‘.‘. I AT/IKI P C p H x l' \ '- B p <1 ' 200 ' _ 3D-alm—Tflolweger-Mlfller """ l- " - — — Tan-im-Tmuo ' 400 '1 ...._ TSD—sim_TVAL - 1 ' - - — T3D-alm— THSRA 1- """" TSD lIm—TKI'IIIIDI Swo '- _ -1 — - my 600 """ I """"" I """"" T fififififififi I """"" - 3 —2 - 1 0 l Loews...) Fig. 15. Comparison of the T—T relation from the simulation, with some often used solar T-T relations. a) The difference in Rosseland T-T relations between the simulation and an Atlas9 atn'iosphere model (Kurucz 1993; Castelli et al. 1997’). b) Differences, 011 a T5000—scale, between the simulation and an Atlas9 model, and four semi—empirical atmosphere models: The model by Holweger 8: Miiller (1974), the classical fit by Krishna Swamy (I966a; 1966b), the HSRA model (Gingerich el al. 1971) and the VAL model (Vernezza el al. 1981). The effective temperature of the solar simulation is Teff = 5777 j: 9K, which 145 is in excellent agreement with irradiance observations (Willson 85 Hudson 1988): 5 777 :1: 2.5 K. The horizontal spatial resolution of this simulation is twice that of the seven other simulations ie, 100 X 100 X 82 points. I compare the resulting T-T relations with various 1D models, in Fig. 15. Panel a) shows temperature differences on the moss—scale and panel 1)) is for the T5000- scale. The vertically hatched area around the zero-line (in both panels) shows the temporal RMS-scatter of the T—T relation of the simulation, confirming that all the differences are statistically significant. Panel a) also shows the difference between the temperature measured on the two T—scales for the simulation (dot—dashed curve), the sign of which means that we can see deeper into the Sun at 5000A than on (a Rosseland) average, and therefore that the Rosseland opacity is larger than the 5000 A opacity. The past decade or so of work on compiling and computing line—data for atoms and molecules has added a lot of line—opacity in the UV, which has increased the Rosseland opacity with respect to the 5000 A opacity. This in turn has increased the difference between TROSS and T5000. This is clearly expressed in the 200 K difference at T 2 g. This difference is two times larger than the corresponding differences among modern atmosphere models (cf. Fig. 15b), so it is no longer justified to assume the two T-scales to be equal. For stellar structure calculations, it has been common practice to use a T-T5000 relation combined with a Rosseland opacity. With the current opacities and today’s demand for accuracy, this no longer seems a valid approximation and I recommend not using it. 146 For the theoretical Atlas9 atmosphere models (Kurucz 1993), both monochro- matic and Rosseland T—T relations can be obtained, illustrated as the dashed line in both panels of Fig. 15. The peculiar wiggles in these curves are features from the Atlas9 model. Using the “overshoot”-option, later rejected by Castelli el al. (1997), these wiggles combine to a larger but smoother dip, compared to the simulation. The rather close agreement in the radiative part of the atmosphere is expected, since the same line-opacities were used, and since the convective fluctuations have only a small effect on the averaged T-T relation above the convection zone (i.e., all averaging methods give the same results, cf. Fig. 14). The model by Holweger 85 Miiller (1974) presented in Fig. 15 (solid line) is about 100 K warmer than the simulation above the photosphere, on both T-scales. For this model the T—difference between the two T-scales is, however, less than 60 K, so below T 2 0.5 the behavior in the two panels differ by the (TROSS — T5000)-difference for the simulation. The Holweger—Muller model is hotter by up to 300 K on the moss—scale, but is very similar in the photosphere and getting increasingly cooler with depth on the TsOOO-scale. The two other semi—empirical atmosphere models presented in Fig. 15, VALC for the quiet Sun (Vernezza at al. 1981) and the Harvard-Smithsonian solar reference atmosphere (HSRA) (Gingerich et al. 1971), differ significantly and essentially in the same way from the simulation results. High in the atmosphere, the differences can most likely be attributed to non—LTE effects, heating from a hot corona and possibly also magnetic fields in the atmosphere, none of which is included in the simulations. It could, however, also be due to misinterpretation of temperature proxies. Since 147 the UV Planck-function, ionization balances and level populations in atoms and ions involve exponential terms in T, the temperature derived from spatial and temporal averaged quantities will be higher than the correspondingly averaged temperature, as shown by Carlsson 85 Stein (1995) in the (:lynamic ID case. This is irrespective of my finding from Fig. 14, that the convective fluctuations have little effect 011 the average T-T relation above the photosphere. Individual lines and the UV brightness can still behave non-linearly. The immediate outcome of this would be larger T, but as the opacity increases strongly with T, the T—scale could easily change enough so as to result in a lower T-T relation than the actual, as found from Fig. 14. The differences between the HSRA, the Holy-veger-h’liiller and the VAL models may be due to temporal variations in the solar atmosphere between the observa- tions, ten years worth of improvement in the handling of non-LTE effects, as well as differences in the opacities used for the T-scale. At intermediate optical depths, from T 2 0.1 down to (T) : Terr. the agreement between both theoretical and semi-empirical atmospheres and the simulation is very good. Differences are also smaller than the difference between the 5000A and the Rosseland T—T relation from the simulation. This agreement is reassuring, as all four approaches to the solar atn'iosphere should be about equally valid in this region: 3D~ effects are small, the velocity—field only contributes a small turbulent pressure to the hydrostatic equilibrium, the convective flux is less than a percent of the total flux, non—LTE effects are small and the hot corona is too far up to have a significant effect here. 148 With the onset of convection, the T-T relations in Fig. 15 diverge, with the various ID—models, sharing the MLT—formulation of convection, all differing from the simulation in more or less the same manner. Below (T) : Terr. 3D—effects become important and the turbulent pressure contributes up to 14% of the total pressure, rendering the simulation the best choice for an atmosphere model. The validity of the simulations in this region has to be assessed from comparisons with more direct observations of the Sun, e.g., measurements of the flux—spectra, limb-darkening and line—profiles. This serves to stress that T—T relations from semi-empirical models are not observations and are only unambiguous with respect to the underlying limb— darkening observations above the point where (T) : Teff' This, of course, rests on the assumption of complete knowledge of the opacities in the atmosphere, which, despite the last decades progress (Kurucz 1992b), still seems a rather unrealistic assumption (Kurucz 1992e; Lester I996). The last T—T relation presented in Fig. 15 is the one by Krishna Swamy (1966a; 1966b), which is still used as an upper boundary in some stellar model codes. It is interesting to note that the behavior below the photosphere is opposite that of the ITIOI'C ITIOCIGI‘II T-T relations. 5. 1 . 7 Fitting formulas T-T relations for seven individual stars are of rather limited value when concerned with stellar structure calculations in general, unless there is a way to interpolate between these individual stars. 149 I have therefore fitted the individual T-T relations, through the corresponding Hopf functions, Eq. (164), to expressions of the form q(T) = (11+ (126—837)“ (108 T - £15) + 77 [(16 - (1777] . (169) where _(Iconv(T) 170 (MT) ‘1’ T ‘1‘ qconv(7-) ( i ) ,7 = is the negative ratio of convective to total Hopf function (cf. Eqs. (156) and (164)]. This definition makes 77 increase with increasing convective ratio of the total flux. The first term serves to give the limit of optical deep layers with radiative energy transport. At large optical depths, the Rosseland opacity suffices to describe the transport of radiation and the stratification is therefore that of a grey atmosphere, provided all the energy is transported by radiation. In terms of the Hopf function, q(T), this means that q(T) —> (100, mimicked by the q1 + T term in Eq. (169). In optically thin layers, the transport of radiation can no longer be described by a single opacity, and the T—T relation deviates from the grey case. The c(2(log10 T — q5) part of the second term provides the asymptotic behavior for T —> 0, and the e‘lq3flq4- factor interpolates smoothly between the two cases. None of the simulations show any sign of leveling off to an isothermal atmosphere, which must be due to the cooling by very strong lines, extending further out than the simulation domains. Therefore, Ie have not included an isothermal term, as is normally done. I still, however, encourage the use of an isothermal boundary condition at T = 0 in stellar structure models. 150 Proper non—LTE calculations would most likely produce temperature minima. for all seven atmospheres, which would be closer to the photosphere than the transition to an isothermal atmosphere, but the global influence on stellar structure models will be small. If some of the energy is transported by convection, the T—T relation will be cooler than the grey case at large optical depth, which is described by the last term in Eq. (169). This term therefore describes the transition from radiative to convective energy transport. It should not be used with normal stellar models, as these often assume a purely radiative T-T relation (set (16 : q7 = 0 in that case), and incorporate convection subsequently. It might seem natural to use the same method as in the envelope models, and use the effective optical depth, T, as given by Eq. (165). This, however, turns out to underestimate the effect of the transition to convection. The reason is that, as pointed out by Nordlund 85 Stein (2000), the high temperature sensitivity of the opacity hides the warm up-flows (granules) from our view, thereby cooling the T—T relation in the convection zone. This is a pure 3D effect due to the large in-homogeneities at the top of the convection zone, which are maintained by the sharp increase in opacity with temperature. The horizontal average of the temperature in the simulation is therefore higher than in the ID—models, although it is lower on a. T—scale. That this is actually the case for the Sun has been supported by comparison with helioseismology (Rosenthal et al. 1999). In 1D models, the 77—terms should be considered part of qconv. Whether they could be included with 7] calculated from the ID fluxes, without compromising consistency, is not resolved. 151 Table 8. Coefficients for Eq. (I71). (1719 am (ing dang, /dY 0.87029 —0.78450 0.17643 0.00256 -1.14200 —0.24250 0.06275 -2.28313 0.44612 —0.94208 —0.70312 4.71555 0.75131 -3.25925 0.21293 —3.29934 2.38410 —11.25957 1.44532 33.49655 1.30317 —82.75698 5.98615 2.66277 1.65473 3.80075 -0.68706 1.86574 Klamsb-OOIQI—‘fi 5.1.8 Results and Discussion The T- and time-averaged temperatures for each of the simulations were fitted to Eq. (169) with standard deviations of 3-7 K (except for Procyon which could not be fitted to better than 14 K—I comment on this below), yielding the coefficients qnm. Each of these coefficients were then fitted to expressions gsurf da 7161 (r — 1:.) (171) ff "1' “iii? 10g fn = an... + an 10g , CITJj) gsurf‘fi: Cl Y where the fn’s are related to the coefficients in Eq. (169) through (mm for n. = I, 4, 5, 6 log (1mm for n = 2, 3, 7 . The numerical values of the coefficients are given in Tab. 8, and the standard deviation between the actual T—T relations and the global fit is 10 K. The relative deviations are everywhere less than 40 K (cf. Fig. 16). With 10 parameters, it is not trivial to make unambiguous fits to Eq. (169). 152 A lot of effort has therefore gone into splitting degeneracies in parameter-space. I furthermore performed the fitting in two different ways: first by fitting the individual T-T relations and then fitting the resulting 70 parameters to Eq. (171), and second by fitting the T—T relations directly from Eq. (171), i.e., making global fits that assume the fins to be linear in Teff and gsurf. I iterated between the two methods, using the results from the first method in the second method and vice versa, until the two sets of results converged. I used four iterations. The accuracy of the fits is illustrated in Fig. 16, from which it is seen that the deviation from the actual T-T relations lies mostly within the shaded area of the temporal RMS—fluctuations. A Comparison with Fig. 15 shows that the deviations are smaller than the differences between the current range of solar atmosphere models and smaller than the differences between the T—T relations from the seven simulations, as shown in Fig. 20. From Fig. 16 it is noticed that Procyon fits the worst. Allowing for a change in position and slope of the zero—line, the n'1ain-component of the disagreement is recogniced as a sharp dip at T slightly smaller than one. Looking at the panels for 7'; B00 and Star B, similar features are noticed, although less pronounced. In an earlier paper, Trampedach el al. (1998) presented a fitting formulae that also contained a negative Gaussian in T. This expression was very good at fitting the above mentioned features, and brought down the RMS deviations to below 8 K, for the individual fits. It proved difficult, however, to parameterize the coefficients for the Gaussian in T93 and gsurf in a way that both fit the simulations and resulted in physically plausible T— T relations outside the immediate Tea/gs,,rf-r‘ar1ge of the simulations. I have therefore abandoned the Gaussian—term in the present work, and obtain improved global fits l I I L l I l I I l I l l I 60 ‘. r :a CenA aCenB nBoo Procyon starA starB Sun I; 40 'j .7 20 7M ,1 '_ "" I ' " A I :4 0.‘"‘rr\ ‘T 111“ . 7IIll .71 VV" '1 71.- 7- g, - I . . v fl v M: -20: r 2 1 I -40 1 r _60 _‘i 1'1 ‘77 I l l I " 1 f I " I V I I I l l L —20 —20 —20 —2O —20 —20 -2O LogT Fig. 16. Plot of the differences between the actual T—T relations and the global fit (solid line) and the individual fits (dashed line). The horizontal solid and dashed lines are the corresponding RMS deviations from the fits. The shaded areas show the RMS of the temporal variation of the temperature. and a more widely applicable fit. The coefficients of the global and the individual fits of the last iteration are compared in Fig. 17. The coefficients listed in Tab. 8 are the results of the last global fit. The coefficients in Tab. 8 are based 011 the seven simulations with X = 70.2960 % and Z 2 1.78785 %. Assuming that a different composition offsets all the coefficients independently of Teff and gsurf, the ands have been changed to correspond to the solar simulation of Sect. 5.1.6 with the more modern abundances X = 73.6945% and Z = 1.8055%. The dq,/dY—term in Eq. (171) is derived from this difference in composition and is therefore only valid for helium changes accompanied by the rather unusual 154 metallicity change, Z — Z6 2 —(Y — YQ)/193.6. This probably proves less useful, but, nonetheless, gives an idea of the changes with composition. Changes along more normal composition-change vectors would of course be very illuminating, but is beyond the scope of the present paper. Notice that my fitting expression, Eq. (169) differ from that presented in Trampedach et al. (1998) in a number of ways. First, as mentioned above, I have abandoned the Gaussian term in order to improve the global fits and making applicable in a wider Tefi/gSUFf—range. Second, I have changed the formulation of the transition from. op— tically thick to thin, to make the coefficients more linearly independent, and third, I have improved the separation of the convective effect on the T-T relation (the last term in Eq. (169)). The standard deviations of the old fits were in most cases larger than those for the present fits. In Fig. 18 I present the behavior of the T—T relation with stellar mass on the zero-age main-sequence (ZAMS). The atmospheric parameters for the ZAMS were derived from the stellar structure models by Christensen-Dalsgaard (1982; 1983), which are also shown in Figs. 1 and 3 in Paper II. The higher mass stars have steeper T-T relations in the photosphere, but they also have a larger curvature, making them shallower further out, as compared to the low—mass stars. The low—mass and higher mass T—T relations, cross at T 2 0.05. Fig. 19 illustrates the gravity dependence of the T—T relation. Going to lower gravity, the T-T relation gets steeper in the photosphere and shallower further out—a similar effect as when going towards higher masses on the ZAMS. In both Figs. 18 and 19 I also show the atmosphere with a grey opacity, i.e., LJJ_I_LLIJ A I 1 L LI I l l J_I nmln I I I I LA I I A I I I I I A I I I I I L I I I L I 1 1 I I 1 I A I I I l I 0.95— -- a CenB e—1.20 . cm; a T > O 90 _ ‘_ starA SUD ‘ 4 _r ,_ _ . starA Sun , i 1.22 4 X ._ a CenA _ . a CenA .> \A\ . 0.85- - N - 1. h 0" _ - O Q. "_ 1" _I.24 b3 . . starB o —1 0°80“ starB -7 t . 7 -- Procyon Procyon, . “I. 77 Boo "' “1.26 0375‘ B00 r . Tl )- \4 ‘ K 0'701‘1H4}:iiI::4{::4%;¢:}:4:§‘.::}73‘r14:I1i:§¢::{:::I4+f :{ii lib—1‘28 1.11 7, Boo '1‘— : 0 9 _ F \T\ . starA t 1 0 _ E a Cen B I 0 8 1 starA i: I I ' 5 \£\ ProcyonEE _ 09? starB I 1;” ; .. 2 :E :0.7 t: . Sun ' 2:3 0.81 1? : D“ 3 3 :E : i 0‘ C9“ A 15 a Cen A ,- 0.6 07-: \e\ 2:— t : I 1: : I '1. r - Sun ‘1)- h : 3; starB L0 5 015-: or CenB \§\ :5- : ' 5 i 35 Proc on: 3 :2 .. .3;\ y ; 0‘51 1 I I J_L LJ 1 1 1 I 1 1 1 I 1 L1. In L j I 1 L1 I -4-_— _ 0'4 4 5 _‘ V I 7 V I' If T r f T‘I’ I I l I I I I I V 1 l V V V l VI—I' I I I I I V I l l I Y I l I V V I I r' I Y ‘ V I 7 I 7 I V 1 '3.68 3.70 3.72 3.74 3.76 3.78 3.80 3.82 . a CenB - log T 4 - I0 cf! 4 starA \L _ 4.0 - - , Sun _ 3.5“ - ID 'l a Crank o~ . I 3_0 —4 .— ‘ starB ; . Procyon- 2.5- _ 17 Boo . 2.0A _ 41-..,...,..fi...,...,...,...,. 3.68 3.70 3.72 3.74 3.76 3.78 3.80 3.82 logIOTeff Fig. 17. Plot comparing the global (triangles) and individually (filled circles with error-bars) fit parameters for the T-T relation—fits. The line segments show the log—loT gradient of the global fits. The error bars are the RMS—scatter from T—T relation fits to single time-steps of the simulations. 156 1.40 M = 0.85 Mo“ 0 90 Mo\\\ 1.30 1.00 M.\\\ 1.10Mo\\\ .. e? 1.20 1.20 Me\\\\\ > 133%\ .. l~ 0 \ L“ 1.10 1.50 MOEQS 1.60 Mo.\\ \\_\ 1 00 ................ 1...?Q.M..‘T. ........... / 0.90 ”””” Fig. 18. The change of the T—T relation with stellar mass, on the zero-age main— sequence. The horizontal dotted line indicates the effective temperature, and the vertical dotted line, optical depth unity. The dashed line shows the grey atmosphere. T(T )/Terr ~3.0 —2.5 —2.0 -1.5 -1.0 -0.5 0.0 0.5 log(TRou) Fig. 1.9. The change of the T-T relation with gravity, for fixed, solar TelT = 5777 K. The horizontal dotted line indicates the effective temperature, and the vertical dotted line, optical depth unity. The dashed line shows the grey atmosphere. 157 §(T/Tefi‘)4 2 g + 7' (dashed curves). The T-T relation are seen to approach the grey case, when going to larger masses on the ZAMS. This is consistent with hotter stars having fewer spectral lines. The gravity sequence in Fig. 19 is a little more complicated. The steady decline of temperature in the high atmosphere, irrespective of gravity, is a sign of strong spectral lines that decouple (have monochromatic T/\ ’2: 1) very high in the atmosphere. Explaining the behaviour at intermediate optical depth, in terms of radiative effects alone, is not possible. Comparing the gravity sequence in Fig. 19 with 1D atmosphere models (Asplund 2003, private communications), rather large differences are found. In 1D the change with gravity is more than a magnitude smaller than in Fig. 19. This might be due to the factor of five extrapolation of the gravity, beyond the range of the seven simulations. Of the five parameters, (11 — ([5, for the radiative T-T relation, only (12 and q3 have their logarithms fitted linearly to logTeq and logg. At closer inspection, (12 is increasing with g and is therefore bounded at small 9 by 0 and around the ZAMS the fit is guided by the simulations. The other parameter, ([3, is decreasing with g and therefore unbound towards the giants. The difference between fitting logqg and (13 itself, is less than a factor of two at logg:3.0, so even with this large an extrapolation, the result is not diverging. 1 keep the logarithmic version, since q3 is required to be positive definit. The fit in the gravity-direction is mostly guided by the 77 Boo—simulation, having the lowest gravity and being only 220 K hotter. The fit at the gravity of 77 B00, but solar Terr, is indeed very close to 77 Boo’s T—T relation. The main difference between the two simulations is the much stronger convection in 771300; The turbulent— to 158 total—pressure ratio is about 20% in the photosphere of the nBoo—simulation and only about 12% in the solar simulation. As pointed out by Asplund et al. (1999), in realistic. convective stellar atmo— spheres the radiative heating and cooling is competing with the adiabatic cooling of rising plasma, expanding from the large density gradient in the atmosphere. The adiabatic stratification is typically more than 1000 l\' cooler than the radiative equi— librium solution. A stratification in-between these two extremes. will therefore ex- perience radiative heating and convective cooling. From Eq. (152) for the radiative heating, it is seen that a larger Planck-mean opacity and a larger density results in more efficient radiative heating. In stars of lower gravity the density will be lower, and the (per mass) opacity red—ward of the Balmer jump will be lower, due to fewer free electrons (from H~ionization) and hence, fewer I“I—-iOIIS. On the other hand, the convective velocities, overshooting into the stably stratified layers, are larger and the adiabatic cooling therefore more efficient. All in all. three effects (density, opacity and velocities) work in the same direction, cooling the atmospheres of stars with lower gravity. With each heating mechanism there is an associated flux F. related to the heating rate Q through (IF (1 z , (173) 6]., Eqs. (152) and (153). The two processes, radiation and convection, produce heating which is propor- tional to the deviation from their respective equilibria stratifications, 1.6., the ra— diative equilibrium with J = S and the adiabatic equilibrium with V : Vad. The stratification which is closest to the adiabatic equilibrium will therefore have the least convective cooling and and the smallest convective overshoot flux (which is negative). This is precisely what is observed in the simulations, with the 77 Boo-simulation dis— playing less than half the overshoot flux of the solar simulation. A detailed analysis will be presented in a future paper. On the backdrop of the analysis above, I therefore feel confident in the T—T rela— tions presented in Fig. 19. The reason that the ZAMS sequence in Fig. 18 approaches the grey atmosphere and not the adiabatic stratification is that the continuum opac— ity increases (as higher levels of hydrogen become more populated) as the lines fade away. The radiative heating therefore remains stronger than the convective cooling and the equilibrium state is determined by the radiative transfer. It is worth noting here that the “approximate overshooting” introduced by (Ku— rucz 1993; Kurucz 1992c) entails a positive overshooting flux, which is at odds with the convection simulations, as well as the basic expression for the convective flux; Fcom, o< (V — Vad), as also discussed by Castelli et al. (1997). 5.1.9 The depth of outer convection zones The depth of an outer convection zone depends in a complex way on the surface boundary conditions. With some simplifications, however, a rough idea of the mecha- nisms involved can still be obtained. I convert the equation of hydrostatic equilibrium 160 to an optical depth scale (1P __3_ = Q 3 (174) (17' K. and integrate from 7' : 0 and inwards Pg 2 9: , (175) K, where R and 7" are some appropriate averages. This results in a first~order estimate of the effects of changing various parts of the physics. The precise form of the averaging is immaterial to the present discussion as it is only concerned with the differential response to changes in the physics. Using some average of the inverse T-T relation for 1;, a relation between T and Pg is obtained. An increase in T(T) will decrease T(T), as the T-T relation is mono- tonically increasing, and will therefore have the same effect on Pg as will an increase in the Opacity. I will now assume that a change in the atmospheric opacity will change the pressure by the same factor in the whole convection zone. The change in depth of the convection zone can be derived from the response to such a pressure change at the bottom of the convection zone. The Schwarzschild criterion for convection to OCCUI‘ Vrad > Vad a (176) is mainly governed by de, as the adiabatic gradient is very close to the ideal— and fully ionized—gas value of Vac; : g at the bottom of deep convection zones. 161 The radiative temperature gradient 3% HFmd Pg 160 9T4 de = (177) will decrease by a decrease in pressure and the bottom of the convection zone will therefore move outward a little. Since \7rad depends strongly 011 temperature, and has a very steep gradient at the bottom of the convection zone, the pressure-change hardly affects the location of Vrad : Vac) 011 the temperature-scale. The largest effect is therefore due to the (almost unchanged) temperature at the bottom of the convection zone occurring at a smaller pressure. If on the other hand (11 is changed, the decrease in pressure is accompanied by a decrease in temperature, more or less counteracting the effect of the smaller pressure. This trend is confirmed by the experiments. I calculated envelope models with small changes (0.001) to o, (11 and the atmospheric 11118”: in order to find the differential changes to the relative depth of the convection zone (Table 9). The magnitude of %elz is smaller than ad—CIZE-E and may even have the opposite sign. From Tab. refenvlist it is noticed that ”7192* and gill? have the same sign, except for StarB and 77 B00, which have opposite and rather small responses to changes in ql and opacity. These same two stars also react the strongest to changes in a. The convective flux in the MLT formulation may be written as 5 (Sp/V1,, P2 , _ Fconv : E I/k.‘ 0’2 T3g/2(V — v )3/27 (178) 162 Table 9. Response to changes in ([1, 11m. and o. adcz adcz adcz 8a aql @1111: Star A 0.0965 —0.0878 0.2404 0.5600 0 Cen B 0.0616 —0.0542 —0.1430 0.3063 Sun 0.1171 —0.1064 -0.2264 0.2861 a Cen A 0.1414 —0.1220 —0.2630 0.3070 Star B 0.2479 0.0192 —0.0684 0.1966 77 B00 0.2763 0.0577 —0.0617 0.2087 Procyon 0.1706 —0.2766 -0.2836 0.1035 (fez/R. 11211116 where k is Boltzmann’s constant, M'u is the atomic mass unit, 6 : —(8lnp/BlnT)p, and it is the mean molecular weight, both of which are fairly unaltered with changing conditions in the atmosphere. The average temperature gradient is called V and V’ is the gradient in the convective cells. Their difference is almost equal to the super-adiabatic gradient, V — Vad. Based on Eq. (178), an increase in temperature and / or a decrease in pressure (brought about by changes to the T—T relation or the atmospheric opacity) will therefore be accompanied by an increase in V — Vad in order to maintain the total flux (and the fixed Teff of the model). This increase of V — Vad corresponds to a decrease of the efficiency of the convection, which will lead to a smaller convection zone. The effect can be counteracted by increasing 0, increasing the efficiency by increasing the distance traveled by convective elements. An increase in the efficiency of convection will enlarge the convection zone, as also seen from. Tab. 9. As far as global observables are concerned, uncertainties in the atmospheric opacities, line—blocking and mismatches between the T-scale and the grey opacity 163 may be hidden in 0 together with pure convection effects. This is one reason that so many values for the solar oz can be found in the literature (Another reason being the lack of concensus on the values for the auxilliary MLT parameters). As illustrated by Tab. 9, such flaws in the treatment of the photosphere will not have the same effect on all stars and an incorrect differential behavior would be expected, possibly masking real convection effects. In Tab. 7 I compare the convection zone depths obtained from envelope models with different choices for the T-T relation but constant, individual 0’s (as found in Paper II). The first row of (1C2 is for the individual T-T relations, derived from the simulations as described in Sect. 5.1.5 and using the same Rosseland opacities in the envelope model as in the computation of the T—T relation. This is the proce- dure I recommend and which separates the radiative contribution to the efficiency of the convection from intrinsic convection effects. The individual T-T relations are compared in Fig. 20. The next three rows of Tab. 7 give the convection zone depth for some of the commonly used approximations. The first uses the solar T—T relation for all the stars, where I use the T-T relation from Sect. 5.1.6, which is based on a simulation with higher resolution and lower Teff than the one listed as “Sun” in Tab. 7. The next approximation uses individual T-T relations based on the n‘ionochromatic, 5000 A, T- scale but still using the Rosseland opacity for the grey opacity in the envelope. The opacities, 1.35000 and H3038 are calculated with the same code, physics and abundances. Last I use the T-T relation from the Harvard—Smithsonian reference atmosphere (Gin— gerich et al. 1971). Fig. 15 shows that (1%0200 and (11531“ for the solar simulation differ 164 111.1 111111111 l 111111111 I1... .I1111 111111111 200 / * 0 I—_| k: L__| \ s. _ . <1 _ _ — a CenB \ ‘ —200- -_-_77Boo _ ' ...-_ Procyon \r - __ _ starA - - starB r -400 ..fi-... . . .,... .2 , , fl , ........ -3 —2 -—1 0 1 L0g(TRoss) Fig. 20. Scaled T—T relations for the seven simulations, relative to the solar simu— lation in the sense: T, X (Egg/Tea.) — lb. The more vigorous the convection, the flatter the T-T relation. from (lcz, approximately in proportion to the corresponding differences in the T-T relation around 7' = 1. This also confirms the trends of the linear analysis listed in Tab. 9, as the various convection zone depths listed in Tab. 7 stem from changes to the T-T relation only. 5.1.10 Conclusion I confirm that the use of T-T relations is indeed a. reasonable way of incorporating the effects of full radiative transfer in stellar structure computations—even in the non—grey case. Based on that I proceed to compute T-T relations for a small number of 3D simu— lations of radiation—coupled convective stellar atmospheres. Each of the T-T relations were fitted to analytical expressions, the coefficients of which were fitted to linear 165 expressions in the atmospheric parameters, Teff and gsurf, for easy implementation in stellar structure codes. I have investigated how changes to the radiative part of the outer boundary af— fects the structure of a star, using the depth of the outer convection zone as a global measure. I evaluated the linear response of the change in depth of the convection zone caused by changes in atn'iospheric opacity, T-T relation and mixing length, respec- tively. My analysis shows that the convection zone is about equally sensitive to the three kinds of changes and, consequently, different parameter—triplets can easily result in the same global properties of a stellar model. References to a particular mixing length are therefore less useful unless accompanied by references to the atmospheric opacity and T-T relation. I also compare the effects of various commonly used assumptions about the T~T relation, and conclude that scaled solar T-T relations introduce systematic effects, while the use of an 5000 A T-T relation with a Rosseland opacity has an even larger effect and should be avoided. To separate the effects of convection from those of the radiative transition in the photosphere, and to avoid unnecessary systematic effects, I recommend a consis— tent usage of T-T relations and their corresponding opacities in stellar structure and evolution calculations. Extrapolation of the parameterized T—T relations towards lower gravity (see Fig. 19) prompted a closer investigation of the interplay between convection and radiation, above the convection zone. I have gained new insights into the physics governing the overshoot of convective flows into the stably stratified parts of the atmosphere; 166 when radiative heating is inefficient, the temperature will approach the adiabatic stratification and the (negative) overshoot flux will diminish. A more detailed analysis will be presented in a future paper. The effect on atmospheric stratification is, however, included in my fit to the T-T relations as presented in Eqs. (169), (171), (172) and Tab. 8, and it is ready to be implemented in stellar structure and evolution codes. As precision and scope of modern observations of stars steadily improve, and as we are entering the age of astero—seismology, higher demands are placed on the mod- eling of stars. With improved understanding and treatment of the interplay between radiation and convection, it will be possible to isolate other effects that so far have been shrouded in the uncertainty of the atmospheric part of stellar models. With improved outer boundary conditions, combined with the mixing—length calibration in Paper 11, we can have more confidence in predictions about the depth of convective envelopes. This, in turn, will allow the study other mixing processes, such as convec— tive overshoot at the base of the convection zone, rotational mixing, g-mode mixing, etc., and compare with observations of chemical enrichment from dredge—ups and the destruction of volatile elements such as Li. 167 5.2 Calibrating the mixing-length The mixing—length parameter a is calibrated using realistic 3D radiation hydrody— namical (RHD) simulations of seven stars in the solar neighborhood of the HR- diagram. The calibration indicates a small variation of o with surface gravity, qurr, and a larger variation with effective ten‘iperature, Terr- I give a simple fit to these results. 5.2.1 Introduction Due to the lack of a better theory of convection in stars, the mixing—length the— ory (MLT) has been used for almost 60 years. By far the largest part of the so— lar convection zone is very close to adiabatic, and the stratification in the bulk of the convection zone is therefore determined by the adiabatic temperature gradient, Va.) 2 (alnT/BlnP)ad. Convection is so efficient that only a very small excess gra— dient, or super-adiabatic gradient, V, = V — Vac) is sufficient for transporting the entire energy flux. In most of the convection zone the super—adiabatic gradient is tiny, V, S, 10’5, which hardly adds up to anything significant even integrating over the whole convection zone. We therefore have no need for a theory of convection here. This picture changes dramatically near the boundaries, especially near the up— per boundary of an outer convection zone. Here convection becomes exceedingly inefficient in transporting the energy flux, as radiative energy transport takes over. For the sun, this layer of appreciable superadiabaticity only takes up the outermost 168 1Mm, just. below the photosphere, in the region where the gas becomes optically transparent and the radiation escapes. This layer, however thin, is crucial for the star as a whole, as it is the stars insulation against the cold of space. With the advances in atomic physics as applied to astrophysics, i.e., the EOS (Hummer 85 Mihalas 1988; Nayfonov et al. 1999; Gong et al. 2001b; Rogers & Nay- fonov 2002; Saumon et al. 1995) and opacities (Seaton 1995; Berrington 1997; Iglesias & Rogers 1996; Kurucz 1992b; Alexander & Ferguson 1994), by far the most uncertain aspects of stellar models are associated with dynamical phenomena: semi-convection, rotational mixing, mixing by g—modes, convective over—shooting and the most promi- nent; convection itself. The present paper is part of an effort to improve on stellar structure models, by using results from a number of 3D convection simulations of stars in the solar neighborhood of the IIR—diagram. The first paper in this series, deals with the radiative part of the stellar surface problem, and presents T—T relations derived from the simulations (Trampedach et al. 2004d, paper I). The present paper is Paper 11. Paper III will address the consequences, of the results from the first two papers, for stellar evolution (Trampedach et al. 2004b). This paper is not a justification of the MLT, nor is it aimed at describing the structure of the surface layers of stars. Rather, I provide a way to use MLT and a non- constant a to correctly model the depth of outer convection zones. MLT in general, and my calibration of a in particular, has limited relevance to stellar atmosphere calculations. The present paper is also a continuation of the work presented in Trampedach 169 et al. (1997), which unfortunately turned out to be flawed by inadequate T—T relations. 5.2.2 Mixing-length vs. 3D convection The conventional interpretation of the mixing—length formulation of convection, is that of bubbles, eddies or convective elements, that are warmer than their surround- ings, rising due to their buoyancy. The Schwartzschild criterion for stability against convection: Vlad < Vad, is equivalent to the statement that convection will occur, QC when gas which is warmer than its surroundings, is buoyant. The rad” and “ad” subscripts indicate the V : dlnT/dlnP in the case of radiative and adiabatic strati— fications, respectively. These bubbles of gas are then envisioned to travel for one mixing-length—hence the name of the formulation—before they dissolve more or less abruptly (Béhm— Vitense 1958). This picture has conceptual problems at the edges of convection zones or in small convective cores, where the distance to the edge is only a fraction of a mixing-length. Most often, the convective elements are also ascribed an aspect—ratio around unity, confounding the problem. The mixing-length is typically chosen to be A : GHQ or oHp, where a is the main free parameter (of order unity) of the formulation, and H is the density— or pressure— scale-height for locally exponential stratifications. It has also been suggested to use I : CYZ, where z is the distance to the top of the convection zone (Canuto & Mazzitelli 1991; Canuto & Mazzitelli 1992). This choice would solve the conceptual problems listed above, but it introduces physical problems since there are strong reasons for 170 real convection to have a stratification similar to an MLT model with A 2: GHQ, as mentioned below. There is also a notion of these convective bubbles travelling in a background of the average stratification. A little like the bubbles rising in a glass of beer, for instance; The cOncept of a. background liquid is rather obvious, with isolated and distinguishable bubbles rising in it. The 3D simulations of convection, on the other hand, display a very different phenomena (see also Stein & Nordlund 1989; Nordlund & Stein 1997). The convection consists of continuous flows; the warm gas rising almost adiabatically, in a background of narrower and faster down—drafts, forced by sheer mass—conservation. A fraction of the up—flows is continuously overturning in order to conserve mass on the back—drop of the steep and exponential density gradient. Locally, the density can be approximated by pocez/HC’. (179) When a vertical “slice” of the up-flow has traveled AZ, the slice would therefore be over—dense by a factor of eAZ/Hé’ , if the up-flow was confined horizontally. There is of course no such confinement, and the fraction, (eAZ/HQ — l) ——> Az/HQ for AZ —> 0, will overturn into the down—drafts. The up-flow will therefore be “eroded” with an e—folding scale of H 9. The result of this concept is actually the same as that for the mixing-length picture described above (with A = ozHQ), but without the same conceptual problems, since the flows 171 are continuous and the overflow from the up-flow and into the down-drafts, likewise happens on a continuous basis. Renaming it the erosion- or dilution-length formulation, it could be a first-order approximation to convection, as observed in the 3D simulations. This is probably the reason that the MLT formulation has worked so well despite its short-comings; It is based on simple mass conservation. The above argument neglects vertical velocity gradients. A positive gradient outwards would accommodate more of the up-fiow and result in a smaller fraction of the tip—flow overturning. With some assumptions on the velocity-gradient, this could be included by means of a factor on the erosion-length, and various geometrical conversions could be included here as well. This is the or being calibrated in the present paper. The 3D simulations also display a nearly laminar Lip-flow, due to the density gradient smoothing out most of the generated turbulence. The down-flows are nar- rower and faster, and since they work against the density gradient, they are also more turbulent. The clown-drafts are not compressed adiabatically, since there is continu- ous entrainment of hot plasma from the neighboring Lip-flows. The down~drafts are therefore super—adiabatic much further in than the Lip-flows which mainly become super—adiabatic from radiative loss of energy around the local T 2 1. There is also a lateral exchange of energy, extending the super—adiabatic peak in the up—flow to larger depth than would have been the case with a purely vertical loss of radiative energy. The super—adiabatic peak produced by the combination of these three phenomena is difficult to reproduce within the MLT frame-work (See Sect. 5.2.6.1). 172 The convective motions in the 3D simulations, are prolific above the convection zone, with the velocity decreasing with a scale—height which is larger than the pressure scale—height. This introduces a new contestant in the atmosphere, and radiative transfer will have to compete with adiabatic cooling for the equilibrium stratification, as discussed in Paper I. The asymmetry in the up—flows and the clown—drafts have some profound effects: In the photosphere, for example, the highly non-linear opacity coupled with the large temperature contrast, results in a T = 1 surface which is very undulated; Over the hot granules, the photosphere is located at larger geometrical height than the cooler inter— granular lanes, and the observed (disk—center) temperature contrast is therefore much smaller than in a horizontal cross—section. This introduces a convective back-warming on the geometrical scale, which has no counterpart on the optical depth—scale. In the convective layers, the density and velocity differences between the up-flows and the down—drafts, give rise to a net kinetic energy-flux Fkin I —Q‘U21)Z , (180) which amounts to about a tenth of the total flux. The assumption of symmetry in the MLT formulation, precludes such a kinetic energy-flux, and is probably the biggest cause for disagreement with the simulations in the deeper, almost adiabatically con— vective layers. As convection quickly approaches the adiabatic stratification, an actual theory of convection is hardly necessary in the bulk of a convection zone that encompasses 173 more than a few pressure—scale—heights. There I only need to determine which adiabat the star is following. Since o is not fixed by the MLT formulation, the answer has to come from “outside” calibrations, e.g., through matching of the radius of a solar model evolved to the present age (Gough & Weiss 1976), or as performed in the present paper. 5.2.3 The simulations The fully compressible RHD simulations are described by Nordlund & Stein (1990), and general properties of solar convection, as deduced from the simulations, are discussed by Stein & Nordlund (1998). Among the code features, important for the present analysis, are radiative transfer with line-blanketing (Nordlund 1982; Nordlund & Dravins 1990), and the transmitting top and bottom boundaries. The bottom is kept at a uniform pressure (but not constant in time), to make a node in the radial pmodes and minimize wave generation by the boundary conditions. The entropy of the in-flowing plasma is adjusted to result in the desired effective temperature, and the outflow is left unchanged. The convection in the simulations, consists of a warm, coherent up—flow, with its entropy virtually unaltered from its value near the bottom of the convection zone. Because of the density gradient, mass conservation forces overturning of the up— flows, on distances of the order of the density scale height. The overturning plasma is entrained into narrow, fast and turbulent, entropy deficient down-drafts, generated by the abrupt cooling in the photosphere. Since only a small part of the convection 174 zone is simulated, open boundaries are necessary for obtaining realistic results. Another requirement for comparison with observations is a realistic treatment of the radiative transfer in the atmosphere, and a corresponding quality of the atomic physics behind the opacities and the equation of state (EOS). Compared to the simulations cited above, I have therefore employed the so-called MHD EOS (Hummer & Mihalas 1988; Dappen et at. 1988), updated most of the continuous opacity sources and added a few new ones, as described in detail in Trampedach (1997). The line opacity is now supplied by opacity distribution functions (ODF) by Kurucz (1992a; 1992b). Each of the simulations were performed on a 50 X50 X 82 grid. After relaxation to a quasi—stationary state, I calculated mean models for the envelope fitting (cf. Sect. 5.2.5). The temporal averaging was performed on a horizontally averaged column density scale, instead of a direct spatial scale, to filter out the main effect of the p—modes excited in the simulations. The seven simulations investigated are listed in Tab. 10. Five of them correspond roughly to actual stars. The chemical composition is, mostly for historical reasons, X : 70.2960 % hydrogen by mass and Z : 1.78785 % metals by mass. In Fig. 21 I have plotted the seven simulations in an I—IR-diagram together with evolutionary tracks (Christensen—Dalsgaard 1982; Christensen—Dalsgaard & Frandsen 1983), for stars with masses between 0.85 and 1.7 1146,, as indicated. The two fictitious stars, Star A and B, can of course be shifted up or down in luminosity, depending on what mass actually corresponds to the given atmospheric parameters. 175 Table 10. depths, dcz. Parameters for the simulations, and derived as and convection zone name a Cen B Star A Sun a Cen A Star B 7) B00 Procyon Spectral class K 1 V K 2 G 2 V G 2 V F 8 G 0 IV F 5 IV—V Teff 5362K 4851 K 5801 K 5768K 6167 K 6023K 6470K loglo gsurr 4.557 4.095 4.438 4.295 4.035 3.703 4.035 .M/[l/IG 0.900 0.600 1.000 1.085 1.240 1.630 1.750 max(Pturb/Ptot) 7.8% 8.2% 10.7% 11.2% 15.5% 19.5% 21.0% a 1.8313 1.8705 1.8171 1.8032 1.7360 1.7383 1.7193 dcz 0.3063 0.5600 0.2861 0.3070 0.1966 0.2087 0.1035 A5 F0 F5 G0 G5 K0 1 l l L 1 l f ‘ 77 B00 7 1.0 — ‘ a — ‘ M=1.7MO .b ~ 1.6 Procyon - i 1.5 I i .40 0.5 _ 1.4 star B _ S _ _ g 1.3 no . . .9. a Cen A ‘ 1.2 _ 0_0 — 1.1 _ q Sun _ star A ‘ 1.0 O _ _ or Cen B L 7 0.9 ” ”0-5 ‘ 0.85 7 I I I I I I I I I I I I I I l I I I I I I I I I l T I I I I 3.95 3.90 3.85 3.80 3.75 3.70 3.65 logic Terr Fig. 21. The position of the seven simulations in the I-IR—diagram. The size of the symbols, reflects the diameter of the star. I have also plotted evolutionary tracks, with masses as indicated, to put the simulations in context. 176 5.2.4 The envelope models The simulations were fitted to 1D, spherically symmetric envelope models (Christensen- Dalsgaard & Frandsen 1983), which extend down to a relative radius of r/R : 0.05, and up to an optical depth of T 2 10‘4. I used the same MHD EOS, and in the atmospheric part of the envelopes I used the same opacities, as in the convection simulations. These atmospheric opac— ities are smoothly joined with the OPAL opacities (Rogers & Iglesias 1992a.) in the temperature interval, T = 2.6—4.1 X 104 K. Convection is treated using the standard MLT as described in Bohm-Vitense (1958), using the standard mixing length [ZO’HP (181) and form factors y = % and z/ = 8. I use the notation introduced by Cough (1977), in which the form-factors are (I) : z//4 = 2 and 77 = 4y/(3\/17) = \[2/9. The photospheric transition from optically thick to optically thin is treated by means of T—T relations derived from the simulations. 1 calculated temporal and TROSS (Rosseland optical depth) averaged temperatures, and fitted these to analytical expressions which were used in the envelope models. The fitting formula and the coefficients are given. in Paper I. The point that I use individual T—T relations instead of scaled solar T-T relations is crucial for the present calibration, as discussed in Paper I. The pressure in the simulations is not purely thermodynamic; turbulent pres— 177 sure also contributes to the hydrostatic equilibrium. I therefore include a turbulent pressure in the envelopes, based on the MLT convective velocities 1)t11r1‘),ll) : (3030.11,? 7 (182) where )3 is a constant. adjusted as part of the calil‘n'ation procedure. as described in Sect. 5.2.5. I suppress the turbulent pressure in the envelopes with a smooth cutoff just above the matching point. to avoid spurious effects on the structure of the envelope model. This has two reasons: The practical one, is that most stellar structure calculations do not include such a turbulent pressure, and a calibration of (1, including 13“,“,le in the whole convection zone. would not apply in these cases. The second, and more physical reason, is that 110),“. in M LT models has a very unphysical behavior, and gives rise to an even less physical PLuerD. The turbulent pressure in the part of the envelope above the matching point would change the outer boundary condition for the envelope, which can have a significant global effect on the model. The velocity drops from an unrealistically high, supersonic maximum. down to zero, in a. fraction of a pressure scale height, giving rise to a devastating pressure gradient. To enable integration of hydrostatic equilibrium, it is necessary to introduce some cutoff, which undoubtedly will introduce an unphysical differential behavior of this MLT turbulent pressure. In the simulations on the other hand, turbulent pressure peaks about half a pressure—scale-height below the top of the convection zone, drops off smoothly both above and below and is non-zero everywhere. 178 I recommend not including a lD-turbulent pressure which is confined to the convection zone, i.e., do not incorporate an overshoot region. Well below the super—adiabatic top of the convection zone, Pturb,1D does however match the turbulent pressure of the simulations rather well, giving an almost differ— entiable match. I take this as evidence, that envelope models including PtuerD, with [3 and a fitted as described in Sect. 5.2.5, gives a realistic extension of the simula— tions towards the center of the star. This fact was exploited in an investigation of convective effects on the frequencies of solar oscillations (Rosenthal et at. 1999) by analyzing eigenmodes in a model combining the simulation and a matched envelope model. I can now proceed with the matching, with confidence. 5.2.5 Matching to envelopes In order to derive as from the simulations, 1 matched horizontal and temporal av- erages of the 3D simulations to 1D envelope models at a common pressure point deep in the simulation. The matching is performed by adjusting [3 till the 3D- and 1D-turbulent pressures agree, and or till the temperatures agree. This method demands a high degree of consistency between the simulations and the envelope models at the matching point, which is the reason for using the exact same EOS (and chemical composition) in both cases, and for including a turbulent pressure in the deep part of the envelope models. The matching is furthermore per- formed at a depth in the simulation, where boundary effects are small and fluctuations in the thermodynamic quantities are small. The latter to ensure that the mean Q, T 179 and Pgas are related by the EOS, i.e., that direct 3D-effects are negligible, as is of course always the case in the envelope models. In order to filter out non-convective effects from this calibration of a, I also demand consistency in the treatment of radiative transfer at the top of the convection zone. I accomplish this by using the Rosseland opacities from the simulations, in the atmospheres of the envelopes, and using T—T relations derived from the simulations, based on this opacity (See Paper I for details). The combination of the average stratification of the simulations in the atmo— sphere, and the matched envelope in the interior, have recently been used by Geor- gobiani et al. (2003b), as a basis for computing the excitation of stellar p—modes for the stars listed in my Tab. 10. 5.2.6 Results and discussion The results of this envelope matching is listed in Tab. 10 where I list both a for the matched envelope, and the relative depth of the convection zone, dcz. The simulations are listed in order of increasing vigor of the convection, translating into decreasing a’s (decreasing efficiency of convection) and decreasing relative depth of the convection zone. As a measure of how vigorous the convection is, I use the maximum of the turbulent to total pressure ratio, as listed in the fifth row of Tab. 10. In Fig. 22 I show the as from the envelope matching as function of Teff) as *— symbols. Those values are plotted with error-bars, corresponding to the RMS scatter 180 I 41 I l I l I I I I J l I l l I I l l I I I l l I l l l 2.0 starA . <> _ \3\ oz Cen B . Sun \Q\ “Ci!” f 1.8 - _. star B 77 B00 ‘ Procyon - :5 t - O 1 6 - K _ I r I I I I I I I I I I I I I l l I I I I I l I I I I I I 3.68 3.70 3.72 3.74 3.76 3.78 3.80 3.82 logiorrerr 1.4 Fig. 22. A plot of the 0’s found from the matching procedure (stars), compared with the linear fit from Eqs. (183) and (184) (triangles). The (lower) diamonds show the 2D calibration by Ludwig et al. (1999), which I have also multiplied by 1.14 to agree with my result for the Sun (upper diamonds). The line-segments show the local LogT—derivative of the fitting expressions. in a derived from envelope matching to individual time-steps of the simulations. This scatter is rather small; 1—3><10—3, and is hard to see in the figure. I fit the derived Oz’s to a simple function of the form 7:3 (sur (la , “’fitrl‘eafllbsr. fl +Blog—J—‘+—/(Y_r@) (183) CINE) gsurf,® (I) 181 and find 1 . QC.) 21.81795 . L“, = —1.32012. (11 A = —1.11208 and B = 0.07454 which is shown as triangles in Fig. 22, with the log T—gradient indicated with line— segments. This linear fit has a. standard deviation of 0 = 0.012, whis is more than an order of magnitude less than the variation of a over the seven stars. It is also an order of magnitude larger than the RMS scatter in a from envelope matching to individual time—steps. This indicates that the variation of (1' is not optimally described by a plane. I find it prudent, however, to keep the fit linear in order to avoid divergences outside my rather limited fitting region of the seven simulations. A similar calibration of 0 against 2D RHD simulations has been performed by Ludwig (I, (1.1. (1999, from here on LFS), using a. method completely independent of ours. In Fig. 22 1 have displayed their fit to their results. as applied to the atmospheric parameters of my simulations (lower set of 0’3, with line—segments indicating local logT-derivatives). Based on the requirement of matching the radius of the present age Sun, they allow for a scaling of their results by 1.1—1.2 to translate from 2D to 3D. My results do indeed agree in the solar vicinity. after a scaling by 1.14, as shown by the upper set of <>—sy1nbols in Fig. 22. The disagreement for Star A is most likely due to differences in the opacities. They base their opacities on the Atlasfi line— 182 IIIj’ IIIT IIII IIII IIII IIII I I I I I l I 1.85 Fig. 23. This figure shows ol’s behavior with Toff and gsurf. The surface is the fit given in Eq. (183). The dashed lines shows evolutionary tracks for stellar masses as indicated. The stars show the 0’s as listed in Tab. 10 and the little circles, connected by lines, show the projection onto the fitting—plane. opacities whereas I use the newer Atlas9 line—opacities (in the form of opacity dis— tribution functions). The difference, as outlined by Kurucz (1992d). consists of the addition of molecular opacity (hydrides and (N, Cg and T10) and improved calcula- tions for the iron—group elements —-*— all in all a. factor of 34 more molecular, atomic and ionic lines. These opacity changes should affect the hotter stars the least, but they still have an effect on the solar model—that was, after all, the main motivation behind 183 the opacity updates (Kurucz 1992b). 1 therefore suspect that the factor translating LFS7 results from 2D to 3D should be closer to 1.21. to bring the Procyon results into agreement. The differences for the other simulations would then be due to the opacity update. LFS used grey radiative transfer in the bulk of the 58 2D simulations going into their analysis, adding another systematic difference (also decreasing with Teff) between my results. The fit by LFS have more degrees of freedom. as warranted by the much larger set of simulations. They also cover a larger area of the 'R.ff/gsurf-(liag1'a1n making non-linear behaviour more pronounced. It is, however, interesting to note that my results are rather well described by a bi—linear fit in logTefr and logg. In the future, the present analysis will be extended to larger temperatures and lower gravities. It seems natural to expect that quantities other than Teff and gsurf, would be more relevant for describing the efficiency of convection. The optical depth at the top of the convection zone. for example, seems much more relevant and directly related to the issue. This quantity, and other related ones, appear to give worse fits, than the one presented above. The reason for this is still unclear. It can be argued that the set of atmospheric parameters I have chosen for my simulations lies close to a line, instead of delineating an area, and therefore only the change of 0' along this line can be significant. Again I refer to the relative agreement with the findings of Ludwig el (1,]. (1999), who explored a larger area in the Teff/gsurf' plane. Also the 77 B00- and Star A-simulations are rather far away from that line and yet they both lie very close to the fitting—plane. In Fig. 2.3 I show the linear fit of Eqs. (183) and (184), as function of both 184 effective temperature and gravity. Contours, in grey, are shown for every 0.01 and labelled for every 0.02. I have also plottet values found from the envelope matching procedure (*—symbols), and connected them to the corresponding points on the plane, with small vertical line—segments. The dashed lines show the evolutionary tracks by Christensen-Dalsgaard (1982) also shown in Fig. 21, with masses decreasing upward as indicated. Fig. 23 shows a very interesting trend; The evolutionary tracks show that the stars on the zero age main sequence have 01’s that depend strongly on mass (0 :5 2.11 — 0.27/VI/Il/IQ), but their a’s converge towards a :3 1.82 going up the Hayashi track at the low—T/low—g corner of the plot. Fig. 23 also shows that the oz Cen binary system has kept an only slightly increasing a—difference between the two components (as the B—component evolves slower), of about 0.03. In a recent calibration of the azCen—system, Morel 6t (11. (2000) find GAB = (186,197), which is a bit larger difference than what 1 find, but it has the same sign. Using individual a’s and T-T relations, evolving with time, might help solve the longstanding problems with the modeling of this system (Fernandes & Neu‘forge 1995; Lydon et‘ al. 1993). 5.2.6.1 Depth of the Solar convection zone I have also run a simulation for direct comparison with solar observations. This simulation has a. more modern composition: a helium mass-fraction of Y@ = 0.245 in accordance with helioseismology (Basu & Antia 1995), and ZQ/XQ = 0.0245 in agreement with meteoritic and solar photospheric metal to hydrogen ratios (Grevesse & Noels 1992). This ratio results in the hydrogen mass-fraction, X : 73.69% and the 185 helium—hydrogen number ratio, He/II=0.0837. I have also increased the resolution to 100 X 100 X 82, to better capture the granulation structure, and I have carefully ad— justed the entropy of the inflowing gas (a constant) to obtain an effective temperature of 5777 31: 14K, in agreement with that derived from solar irradiance observations: Ten“ 2 5777 i 2.5K, (Willson & Hudson 1988). Matching this simulation to an envelope—model, gives a : 1.8670. [3 :2 0.75237 and a depth of the solar convection zone, dcz : 0.2869 :l: 0.0009 R5,. This is in good agreement with that inferred from inversion of helioseismic observations: ([02 2: 0.287 i 0.003 RC.) (Christensen-Dalsgaard et al. 1991). The uncertainty I quote for my result is merely the RMS scatter resulting from performing the full fitting of T-T relations and envelope-matching for the individual time-steps of the simulation. As indicated below Eq. 181, there are two more parameters to standard MLT. The two form-factors, the aspect ratio of convective elements, (I), and 7}, which is related to the radiative exchange of energy between the up- and down—flows, are, however, not linearely independant and I therefore limit my discussion to 7]. Fitting 77 with respect to the height of the super—adiabatic peak, I get a = 1.9618, 13 : 0.72792 and 7} = 0.065948, resulting in a. convection zone of depth 0.2872 Hg. This value is still safely within the uncertainty of the helioseismic result. The height of the peak in the super-adiabatic gradient is increased from 0.5 in the model above, to 0.7 with the new value of 7}. If on the other hand I adjust the form—factor, 7], so as to obtain the same log- arithmic temperature gradient, V, at the matching point, then I get of = 3.8545, /3 : 0.46403 and 77 : 6.3943 X 10'4, and a convection zone of depth 0.2894 HQ. This 186 result is not within the uncertainty of the helioseismic result. Furthermore, the peak of the super-adiabatic gradient becomes unphysically large, reaching a value of 2, about 100 km below the photosphere. That T, Q and V cannot be simultaneously matched at a. common pressure-point (with plausible parameters), indicates that the MLT formulation converges rather slowly, if ever, towards the super-adiabatic gradient, V — Wad, of a. real convective envelope. This might be due to the neglect of kinetic—flux in the MLT formulation, as detailed in Sect. 5.2.2. Notice that the depth of the solar convection zone, as found above, results from ab initio calculations, from the E08 and opacity calculations, to the RHD simulations. The adjustable parameters that enter the simulations are the resolution, the viscosity coefficients, and the size of the time step relative to the Courant time. These are tuned to resolve the thermal boundary layer at 7' = 1 and the convective structures, to minimize numerical diffusion, and to minimize the computing time. Nothing is adjusted to fit solar observations. 5.2.7 Conclusion I have calibrated the MLT parameter a, by matching 1D envelope models with 3D RHD simulations, and established a significant variation of a with stellar atmospheric parameters Teff and gsurf. My results point to a common value of oz ’1 1.82 at the beginning of the ascend up the I—Iayashi track and a decreasing with mass on the zero-age main-sequence. 187 There is of course still the possibility that 1, despite my efforts, have overlooked a source of systematic errors in my calibration, but the absolute agreement with the seismologically inferred depth of the solar convection zone, found in Sect. 5.2.6.1, strengthens my confidence. Although various values of a have been considered in the modeling of stellar evolution, an or varying during the evolution of a star has, to my knowledge, not been tried yet. Results of such evolution calculations are presented in Paper 111. I stress that the choice of a depends on the choice of atmospheric physics, 27.6. T-T relation and atmospheric opacity. Employing a scaled solar T-T relation will alter the effect of oz, as shown in Paper I. I recommend that my fitting formula, Eq. (183) be used with individual T-T relations. As ground—based and soon also space-borne, asteroseismology is beginning to provide strong constraints on the structure of stars, other than the Sun, stronger demands are placed on the theoretical models. An absolute calibration of the mixing-length parameter, a, is the first step to- wards improving the treatment of convection in stellar structure models. A funda— mentally improved formulation of convection is of course desireable, but has proven rather difficult to come by. Various attempts have been made to rectify this sit— uation. Canuto (1993) present a formulation based on fully developed turbulence, which however, does not account for the steep density gradient and the inherent as- symetry between up- and down—flows. (Lydon et al. 1992) base their model on 3D hydrodynamical simulations of convection, and this is probably the most promising way forward. A number of approximations render their results less than optimal for 188 the next generation of convection models, however. With the connection between MLT and the 3D convection simulations, found in Sect. 5.2.2, I believe a properly calibrated mixing-length formulation, with the mixing—length being proportional to the density or pressure scale-height, to be the best choise for the time being. 189 Chapter 6 Conclusions 6.1 Summary Numerical instabilities first discovered in a deep solar simulation, have severely ham— pered progress in production runs. The reasons and possible remedies are explored in chapter 2.2, but no conclusions have been reached yet. The research has led to a small number of possible solutions that will be implemented and explored in the near future. The two equation of state projects, most popular with the astrophysical com— munity, were compared in much detail in Sect. 3.2. This analysis revealed a number of differences and considering the success of the OPAL EOS, there seems to be room for improvements on behalf of the MHD EOS. Analyzing the \IJ—term in the MHD EOS, which is an approximation to the second—order (in density) hard-sphere interaction between neutral particles, it was found that it has very little effect under solar circumstances, and was not the main 190 reason for the very successful pressure-ionization in the .\'III D EOS. The term is still suspicious, however, as hard-sphere interactions. at best, can be a crude approxima- tion to the real phenomena at play. In Sect. 3.4 the lessons learned from the EOS comparisons, were applied to im— provements of the MHD EOS. In particular, the previously published improvements of post—Holtsmark micro—field distributions (Nayfonov (:1 al. 1999) and relativisti— cally degenerate electrons (Gong 61 al. 2001b). were incorporated. together with the two quantum effects; exchange interactions between identical particles and quantum diffraction. For a solar stratification of density and temperature. the largest effect turned out to stem from a proper treatment of higher order Coulomb interactions, beyond the Debye-Hiickel theory. In order to extend the validity of the MHD EOS to the domain of stellar at— mospheres. provisions for including molecules by means of parameterized partition- function, was added. We presently use a data-base containing 315 (Ii-atomic and 99 poly-atomic molecules. A new EOS will also have consequences for opacities, through changed popula— tions of the energy levels, and changes to the dissociation and ionization equilibria. From the OP—team (Seaton 2003, private communications) it has been confirmed that there has been significant improvements to the atomic physics and time is getting ripe for a new opacity calculation. The improvements to the equation of state and the opacities, are accompanied by a comparable improvement to the method for evaluating the radiative transfer. By carefully selecting a small number of wavelength points, it is possible to reproduce 191 the result of the full radiative transfer to within a. percent. A robust algorithm for selecting the wavelength points is presented, and tested on a handful of conventional 1D stellar atmosphere models. The new method is shown to be more stable against the convective fluctuations, i.e., when applied to a vertical slice of a snapshot of a solar simulation, the overall agreement with the full calculation is improved compared to the opacity—binning scheme currently employed. In chapter 5 a set of convection simulations for stars in the (larger) solar neigh- borhood of the HR—diagram, are analyzed in order to improve conventional 1D stellar structure models. It is found that T-T relations can indeed describe stellar atmo— spheres correctly, with just a single (Rosseland) opacity—even in the case of dis— tinctly non-grey atmospheres. It is found that there are significant adiabatic cooling and radiative heating, taking place, when convective overshoot is present (as is always the case in the simulations). The stratification of an atmosphere is not necessarily in radiative equilibrium. Having isolated the radiative effects through a properly defined radiative T- 7' relation we proceed to present a calibration of the main free parameter, oz, of the mixing-length formulation of convection. A small but significant change with Teff and gsurf is found, and for the solar case, the calibration results in a very good agreement with the depth of the solar convection zone—an absolute calibration, with no free parameters to adjust. 192 6.2 Future Work 6.2.0.1 The convection simulations The various improvements out-lined in the present work are now ready for implemen— tation. The new MHD EOS is coded and is ready to run production runs. During the production of EOS tables, 1 will continue the newly opened negotiations with the OP-team and the molecular line database-team, based in Copenhagen and Uppsala. My goal is a compilation and computation of “unified” opacities, that are equally relevant for stellar structure-, and for stellar atmosphere—research. At least in the at— mosphere, this EOS and opacity combination should be valid for stars ranging from cool, late type stars and all the way to white dwarfs. The issue of stability in the convection simulations, as explored in chapter 2.2, will also be addressed in the near future, to facilitate early resumption of production runs. There are a number of potential solutions to explore, of which the piece-wise cubic—spline interpolation will probably be the first one to be tested. A candidate solution will have to prove its worth in a number of tests with analytical solutions, such as simple advection tests of different shapes, as well as sound—waves and shocks of varying amplitude. As soon as a solution has materialized, the simulations can be restarted, and that probably with the new EOS. The next step is implementation of the SOS radiative transfer—scheme and direct comparison between that and the current opacity-binning scheme and various stages in between. When properly tested, this will be applied in the production runs. The last update will probably be the opacity calculations, since the computation 193 of absorption coefficients of the last 11 elements included in the OP and IP projects, is still work in progress. With all the improvements in place, and validated by a handful of convection simulations of individual stars, work on the grid of 3D convection simulations can begin. This will pave the way for many new insights, much improved analysis and interpretation of observations, and the possibility for performing “measurement” on the simulations that are not (yet) possible for real stars. 6.2.0.2 Applications There are a great number of pressing issues, concerning the remaining space-based astero-seismology missions. One issue is the amount of convective “noise” that can be expected from other stars, potentially drowning a p—mode signal. Another issue is the damping of modes, broadening the profiles of the eigen—modes. Preliminary results (Kjeldsen & Bedding, 2003 private comnuinications) points to higher than expected damping, and maybe even to the point of splitting the line into several independent random components. The convection simulations will also be useful for simulating the signal that will be observed, to improve on the understanding and interpretation of exactly what the instruments observe. This also applies to the identity, or global parameters, of the observed stars, as the convection simulations are excellent tools for abundance determinations (Asplund 2000; Asplund 615 (1!. 2000b) and for determining Tefl' and gsurf from spectral line—shapes. When the atmospheric parameters have been determined, the interior of the star 194 can be investigated by means of matching 1D envelope models to the bottom of the convection simulations, and e.g., evaluating the eigen—modes for the combined model, as has previously been done for the Sun (Rosenthal 61. al. 1999). As shown in chapter 5, 1D stellar models can be calibrated against the simulations, 7.6. the atmospheric structure and mixing-length. which can then be fed into stellar evolution calculations, to gain further insights into the lives of stars. The possibilities are as many, and as far—reaching, as our imaginations can take 195 Bibliography Abe, R. 1959, “Giant cluster expansion theory and its application to high temperature plasma”, Prog. Theor. Phys. 22(2), 213 Aguilar, A., West, J. B., Phaneuf, R. A., Brooks. R. L., Folkmann, F., Kjeldsen, H., Bozek, .1. D., Schlachter, A. S., Cisneros, C. 2003, “Photoionization of isoelectronic ions: Mg“r and A12“, Phys. Rev. A 67. 12701 Alastuey, A., Cornu, F., Perez, A. 1994, “Virial expansion for quantum plasmas, diagrammatic resummations”, Phys. Rev. E 49, 1077 Alastuey, A., Cornu, F., Perez, A. 1995, “Virial expansion for quantum plasmas, maxwell-boltzmann statistics”, Phys. Rev. E 51, 1725 Alastuey, A., Perez, A. 1992, “Virial expansion of the equation of state of a quantum plasma”, Europhys. Lett. 20, 19 Alexander, D. R., Ferguson, J. W. 1994, “Low-temperature Rosseland opacities”, ApJ 437, 879 Allard, N. F., Kielkopf, J., Feautrier, N. 1998, “Satellites on the Lyman 13 line of atommic hydrogen due to H—I-I+ collisions”, A&A 330, 782 Allard, N. F., Koester, D., Feautrier, N., Spielfiedel, A. 1994, “Free—free quasi— molecular absorption and satellites in Lyman—alpha due to collisions with H and H+”, A&AS 108,417 Asplund, M. 2000, “Line formation in solar granulation. III. The photospheric Si and meteoritic Fe abundances”, A&A 359, 755 Asplund, M., Gustafsson, B., Kiselman, D., Eriksson, K. 1997, “Line—blanketed model atmospheres for R Coronae Borealis stars and hydrogen—deficient carbon stars”, A&A 318, 521 Asplund, M., Nordlund, A., Trampedach, R.., Allende Prieto, C., Stein, R. F. 2000a, ‘ I o c a I c I ‘Line formation 1n solar granulation. 1. Fe line shapes, shifts and asymmetries”, 196 A&A 35.9. 729 Asplund, M., Nordlund, A., Trampedach, R.., Stein, R. F. 1999, “3D hydrodynamical atmospheres of metal-poor stars. Evidence for a low primordial Li abundance”, A&A 346, L17 Asplund, M., Nordlund, A., Trampedach, R., Stein, R. F. 2000b, “Line formation in solar granulation. II. The photospheric Fe abundance”, A&A 359, 743 Asplund, M., Trampedach, R., Nordlund, A. 1998, “Confrontation of stellar surface convection simulations with stellar spectroscopy”. In (Ciménez et al. 1998), 221 Badnell, N. R., Seaton, M. J. 2003, “On the importance of inner—shell transitions for opacity calculations”, J. Phys. B: At. Mol. Opt. Phys. (submitted) Basu, S., Antia, H. M. 1995, “Helium abundance in the solar envelope”, M.N.R.A.S. 276,1402 Basu, 8., Antia, H. M. 1997, “Seismic measurement of the depth of the solar convec- tion zone”, M.N.R.A.S. 287, 189 Basu, S., Christensen-Dalsgaard, J. 1997, “Equation of state and helioseismic inver- sions”, A&A 322, L5 Basu, S., Dappen, VV., Nayfonov, A. 1999, “Helioseismic analysis of the hydrogen partition function in the solar interior”, ApJ 518, 985 Basu, S., Thompson, M. J. 1996, “On constructing seismic models of the Sun”, A&A 305, 631 Berrington, K. A. (ed.) 1997. The opacity project, Vol. 2. Institute of Physics Publishing Berthomieu, (3., Cooper, A., Cough, D., Osaki, Y., Provost, J., Rocca, A. 1980, “Sensitivity of five minute eigenfrequencies to the structure of the sun”. 111: H. A. Hill and W. Dziembowski (eds), Nonradial and Nonlinear Stellar Pulsation, Vol. 125 of Lecture Notes in Physics, IAU C011. 38, Springer Verlag, Berlin, 307 Bohm—Vitense, E. 1953, “Die Wasserstofikonvektionszone der sonne”, Zs. f. Astroph. 32,135 Bohm-Vitense, E. 1958, “IIber die Wasserstoffkonvektionszone in Sternen ver— schiedener Effektivtemperaturen und Leuchtkrafte”, ZS. f. Astroph. 46, 108 Borysow, A., Frommhold, L. 1989, “Collision—induced infrared spectra of 112-He pairs 197 at temperatures from 18 to 7000 k. 11. Overtone and hot bands”, ApJ 341, 549 Borysow, A., Frommhold, L., Moraldi, M. 1989, “Collision-induced infrared spectra of Hg-He pairs involving 0 H 1 vibrational transitions and temperatures from 18 to 7000 k”, ApJ 336, 495 Canuto, V. M. 1993, “Turbulent convection with overshooting: Reynolds stress ap— proach. ii.”, ApJ 416, 331 Canuto, V. M., Mazzitelli, I. 1991, “Stellar turbulent convection: A new model and applications”, ApJ 370, 295 Canuto, V. M., Mazzitelli, I. 1992, “Further improvements of a new model for tur— bulent convection in stars”, ApJ 389, 724 Carlsson, M., Stein, R. F. 1995, “Does a nonmagnetic solar chromosphere exist?”. A p] 440, L29 Carlsson, M., Stein, R. F. 1997, “Formation of Solar Calcium H and K Bright Grains”, ApJ 481, 500 Castelli, F., Gratton, R. G., Kurucz, R. L. 1997, “Notes on the convection in the ATLAS9 model atmospheres”, A&A 318(3), 841 Cauble. R., Perry, T. S., Bach, D. R.., Budil, K. S., Hammel, B. A., Collins, G. W., Gold, D. M., Dunn, J., Celliers, P., Silva. L. B. D., Foord, M. E., Wallace, R. J., Stewart, R. E., , Woolsey, N. C. 1998, “Absolute equation-of—state data in the 10-40 Mbar (1—4 TPa) regime”, Phys. Rev. Letters 80, 1248 Chabrier, G., Baraffe, I. 1997, “Structure and evolution of low—mass stars”, A&A 327, 1039 Chase, M. W., Curnutt, J. L., Downey, J. R.., McDonald, R. A., Syverud, A. N., Valenzuela, E. A. 1982, “JANAF thermochemical tables, 1982 supplement”, J. Phys. Chem. Ref. Data 11, 695 Christensen—Dalsgaard, J. 1982, “On solar models and their periods of oscillation”, M.N.R.A.S.199, 735 Christensen-Dalsgaard, J. 1991, “Solar oscillations and the physics of the solar inte- rior”. In (Cough & Toomre 1991), 11 Christensen-Dalsgaard, J., Dappen, W. 1992, “Solar oscillations and the equation of state”, A&AR 4(3), 267 198 Christensen-Dalsgaard, J., Dappen, WK, Ajukov, S. V.. Anderson, E. R.. Antia, H. M., Basu, S., Baturin, V. A., Berthomieu, C., Chaboyer, B., Chitre, S. M., Cox, A. N., Demarque, P., Donatowicz, J., Dziembowski, W. A., Gabriel, M., Cough, D. O., Cuenther, D. B., Cuzik, J. A., Harvey. J. W., Hill, F., Houdek, C., Iglesias, C. A., Kosovichev, A. C., Leibacher, J. W., Proffitt, P. M. C. R... Provost, J., Reiter, J., Rhodes M., E. J., Rogers, F. J., Roxburgh. I. W., Thompson, M. J., Ulrich, R. K. 1996, “The current. state of solar modeling”, Science 272(5266), 1286 Christensen—Dalsgaard, J., Dappen, \N., Lebreton, Y. 1938. “Solar oscillation fre- quencies and the equation of state”. Nature 336. 634 Christensen-Dalsgaard. J., Frandsen. S. 1983. “Radiative transfer and solar oscilla— tions”, Sol. Phys. 82, 165 Christensen—Dalsgaard. J.. Cough. D. 0.. Thompson. M. J. 1991, “The depth of the solar convection zone”, ApJ 378. 413 Cooper, M. S., DeWitt, H. E. 1973, “Degeneracy effects in gases in the near—classical limit”, Phys. Rev. A 8(4), 1910 Cox, J. P., Cuili, R. T. 1968, Physical principles. Vol. 1 of Principles of Stellar Structure. Gordon and Breach, Science Publishers Dabrowski, I. 1984-, “The Lyman and Werner bands of H2”, Canadian .1. Phys. 62, 1639 Dappen, W'. 1992. “The equation of state for stellar envelopes: Comparison of theo— retical results”, Rev. Mex. Astron. Astrofis. 23. 141 Dappen, W. 1996, “Helioseismic diagnosis of the equation of state”, Bull. Astr. Soc. India 24, 151 Dappen, W., Anderson, L., Mihalas, D. 1987, “Statistical mechanics of partially ionized stellar plasmas: The Planck-Larkin partition function, polarization shifts, and simulation of optical spectra”, ApJ 319, 195 Dappen, W., Lebreton, Y., Rogers, F. 1990, “The equation of state of the solar interior: A comparison of results from two competing formalisms”, Solar Physics 128, 35 Dappen, W., Mihalas, D., Hummer, D. C., Mihalas, B. W. 1988, “The equation of state for stellar envelopes. III. Thermodynamic quantities”, ApJ 332, 261 de Boor, C. 1978, A practical guide to splines. Springer Verlag 199 Debye, P., Hiickel, E. 1923, “Zur Theorie der Electrolyte”, Physic. Zeit. 24(9). 185 DeWitt, H. E. 1961, “Thermodynamics functions of a partially degenerate, fully ionized gas”, J. Nucl. Energy, Part C: Plasma Phys. 2, 27 DeWitt, H. E. 1969, “Statistical mechanics of dense ionized gases”. In: S. Kamar (ed.), Low Luminosity Stars, Cordon and Breach. New York. 211 DeWitt, H. E., Schlanges, M., Sakakura, A. Y., Kraeft, W. D. 1995, “Eos”, Phys. Lett. A 30. 326 Di Mauro, M. P.. Christensen-Dalsgaard. J., Rabello—Soares. M. C., Basu, S. 2002. “Inferences on the solar envelope with high degree—modes”. A&A 384, 666 Ebeling, W., Férster, A., Fortov, V. E.. Cryaznov. V. K.. Polishchuk, Y. A. 1991, Thermodynamic properties of hot dense plasmas. Teubner, Stuttgart, Germany Ebeling, W., Kraeft, W., Kremp, D. 1976. Theory of bound states and ionization equilibrium in plasmas and solids. Akademie Verlag, Berlin, DDR Eggleton, P. P., Faulkner, J., Flannery, B. P. 1973, “An approximate equation of state for stellar material”, A&A 23, 325 Elliot, .1. R.., Kosovichev, A. C. 1998, “Relativistic effects in the solar equation of state”. In: S. C. Korzennik (ed.), Structure and Dynamics of the Interior of the Sun and Sun-like Stars, SOHO 6/CONC 98 Workshop. ESA, Noordwijk, The Netherlands, 453 Fermi, E. 1924, “I-Iejsa”, Zeits. f. Phys. 26, 54 Fernandes, J., Neuforge, C. 1995. “a Centauri and convection theories”, A&A 295, 678 Fowler, R., Guggenheim, E. A. 1956, Statistical thermodynamics. Cambridge Uni- versity Press Fritsch, F. N., Butland, J. 1984, “A method for constructing local monotone piecewise cubic interpolants”, SIAM J. Sci. Stat. Comput. 5, 300 Ceorgobiani, D., Stein, R. F., Nordlund, A., Trampedach, R. 2003a, “What causes p-mode asymmetry reversal?”, ApJ (in preparation) Ceorgobiani, D., Trampedach, R.., Stein, R. F., Ludwig, H.—C., Nordlund, A. 2003b, “Excitation of stellar p-mode oscillations”, ApJ (in preparation) Ciménez, A., Cuinan, E. F., Montesinos, B. (eds) 1998. Theory and tests of convec- tion in stellar structure, ASP conf. series, First. Cranada Workshop Cingerich, O., Noyes, R. W., Kalkofen, W. 1971. “The Harvard—Smithsonian reference atmosphere”, Sol. Phys. 18(3), 347 Cong, Z., Dappen, W., Nayfonov, A. 2001a, “Effects of heavy elements and excited states in the equation of state of the solar interior”. ApJ 563, 419 Cong, Z., Dappen, W., Zejda, L. 2001b, “Mhd equation of state with relativistic electrons”, ApJ 546, 1178 Cong, Z., Zejda, L., Dappen, \‘V.. Aparicio, J. M. 2001, “Ceneralized Fermi—Dirac functions and derivatives: properties and evaluation”, Comp. Phys. Comm. 136, 294 Cough, D. O. 1977, “hi-fixing-length theory for pulsating stars”. ApJ 214, 196 Cough, D. O., Toomre, J. (eds) 1991. Challenges to theories of the structure of moderate—mass stars. Vol. 388 of Lecture Notes in Physics. Berlin, IAU Coll. 38. Springer Verlag Cough, D. O., \A/eiss, N. O. 1976. “The calibration of stellar convection theories”, M.N.R.A.S. 176. 589 Craboske, H. C., Harwood, D. J., Rogers, F. J. 1969, “Thermodynamic properties of nonideal gases. 1. Free—energy minimization method”, Physical Review 186(1), 210 Craboske, II. C., Olness, R. J., Crossman. A S. 1975. “Thermodynamics of dense hydrogen—helium fluids”, ApJ 199. 255 Crevesse, N., Noels. A. 1992, “Cosmic abundances of the elements”. In: N. Prantzos, E. Vangioni—Flam, and M. Cassé (eds), Origin and Evolution of the Elements, Cambridge University Press, 15 Crevesse, N., Noels, A., Sauval, A. .1. 1996, “Standard abundances”. In: S. S. Holt and C. Sonneborn (eds), Cosmic Abundances, ASP, 117 Custa'fsson, 13., Bell, R. A., Eriksson, K., Nordlund, A. 1.975, “A grid of model atmospheres for metal-deficient giant. stars I”, A&A 42, 407 Harten, A. 1983, “High resolution schemes for hyperbolic conservation laws”, J. Comp. Phys. 49(3), 357 201 Hauschildt, P. H., Allard, F., Baron, E. 1999a, “The NEXTGEN model atmosphere grid for 3000 S ten 3 10,000k”, ApJ 512, 377 Hauschildt, P. H., Allard, F., Ferguson, J., Baron, E., D. R. A. 19991), “The NEXTGEN model atmosphere grid. 11. Spherically symmetric model atmospheres for giant stars with effective temperatures between 3000 and 6800 K”, ApJ 525, 871 Heisenberg, W. 1927, “Uber die grundprinzipien der quantenmechanik”, Forsch. und Fortschr. 3(11), 83 Henyey, L., Vardya, M. S., Bodenheimer, P. 1965, “Studies in stellar evolution. III. The calculation of model envelopes”, ApJ 142(3), 841 Herzberg, C., Howe, L. L. 1959, “The Lyman bands of molecular hydrogen”, Cana- dian J. Phys. 37, 636 Holtsmark, J. 1924, “Hejsa”, Phys. Z8. 25, 73 Holweger, H., Muller, E. A. 1974, “The photospheric barium spectrum: Solar abun— dance and collision broadening of BaII lines by hydrogen”, Sol. Phys. 19, 19 Hooper, C. F. 1966, “Electric microfield distributions in plasmas”, Phys. Rev. 149, 77 Hooper, C. F. 1968, “Low—frequency component electric microfield distributions in plasmas”, Phys. Rev. 165, 215 Hubeny, 1., Mihalas, D., Werner, K. (eds) 2003. Stellar atmosphere modeling, Vol. 288 of ASP conf. series, San Francisco Huebner, W. F. 1986, Physics of the sun, Vol. 1, 33. D. Reidel Publishing Co., Dordrecht Hummer, D. C., Berrington, K. A., Eissner, W., Pradhan, A. K., Saraph, H. E., Tully, J. A. 1993, “Atomic data from the IRON project. 1: Goals and methods”, A&A 279, 298 Hummer, D. C., Mihalas, D. 1988, “The equation of state for stellar envelopes. I. An occupation probability formalism for the truncation of internal partitition functions”, ApJ 331, 794 Hunter, C., Yau, A. W., Pritchard, H. O. 1974, “Rotation-vibration level energies of the hydrogen and deuterium molecule-ions”, Atomic Data and Nuclear Data Tables 14, 11 Iglesias, C. A., Rogers, F. J. 1991, “Opacity tables for Cepheid variables”, ApJ 371, L73 Iglesias, C. A., Rogers, F. J. 1995, “Discrepancies between OPAL and OP opacities at high densities and temperatures”, ApJ 443, 460 Iglesias, C. A., Rogers, F. J. 1996, “Updated opal opacities”, ApJ 464, 943 Iglesias, C. A., Rogers, F. J., Wilson, B. C. 1987, “Reexamination of the metal contribution to astrophysical opacity”, ApJL 322, L45 Iglesias, C. A., Rogers, F. J., Wilson, B. C. 1992, “Spin-orbit interaction effects on the Rosseland mean opacity”, ApJ 397, 717 Irwin, A. W. 1981, “Polynomial partition function approximations of 344 atomic and molecular species”, ApJS 45, 621 Irwin, A. W. 1987, “Refined diatomic partition functions. I. Calculational methods and Hg and CO results”, A&A 182, 348 Jergensen, U. C. 2003, “Molecules in stellar and star—like atmospheres”. In (Hubeny et al. 2003), 303 Kippenhahn, R., VVeigert, A. 1992, Stellar structure and evolution. Springer Verlag. Chp. 18.4 Kjeldsen, H., Hansen, J. E., Folkmann, F., Knudsen, H., West, J. B., Andersen, T. 2001, “The absolute cross section for L-shell photoionization of C+ ions from threshold to 105 eV”, ApJS 135, 285 Kjeldsen, H., Kristensen, B., Brooks, R. L., Folkmann, F., Knudsen, H., Andersen, T. 2002a, “Absolute, state-selective measurements of the photoionization cross sections of N+ and O+ ions”, ApJS 138, 219 Kjeldsen, H., Kristensen, B., Folkmann, F., Andersen, T. 20021), “Measurements of the absolute photoionization cross section of Fe+ ions from 15.8 to 180 eV”, J. Phys. B: At. Mol. Phys. 35, 3655 Koester, D., Finley, D. S., Allard, N. F., Kruk, J. W., Kimble, R. A. 1996, “Quasi— molecular satellites of Lyman beta in the spectrum of the DA white dwarf Wolf 1346”, ApJ 463, L93 Kovetz, A., Lamb, D. Q., Horn, H. M. V. 1972, “Exchange contribution to the thermodynamic potential of a partially degenerate semirelativistic electron gas”, ApJ 174. 109 203 Krae'ft, W., Kremp, D., Ebeling, W., R6pke, C. 1986, Quantum statistics of charged particle systems. Plenum, New York Krishna Swamy, K. S. 1966a. “Profiles of strong lines in C and K dwarfs”, Ph.D. thesis, University of California at Berkeley. Presented in Henyey et al. (1965) Krishna Swamy, K. S. 1966b, “Profiles of strong lines in K-dwarfs”, ApJ 145, 174 Kurucz, R. L. 1992a, “Atomic and molecular data for opacity calculations”, Rev. Mex. Astron. Astrofis. 23, 45 Kurucz, R. L. 1992b, “Finding the ”missing” solar ultraviolet opacity”, Rev. Mex. Astron. Astrofis. 23, 181 Kurucz, R. L. 1992c, “Model atmospheres for population synthesis”. In: B. Barbuy and A. Renzini (eds), The stellar population of galaxies, IAU, 225 Kurucz, R. L. 1992d, “Model Atmospheres for Population Synthesis”. In: B. Bar- buy and A. Renzini (eds), The Stellar Populations of Galaxies, IAU Symp. 149, Springer Verlag, Dordrecht, 225 Kurucz, R. L. 1992e, “Remaining line opacity problems for the solar spectrum”, Rev. Mex. Astron. Astrofis. 23, 187 Kurucz, R. L. 1993. Atlas9. Kurucz CD-ROM No. 13. Stellar atmosphere programs and 2km s"1 grid Kurucz, R. L. 1995, “A new opacity—sampling model atmosphere program for arbi— trary abundances”. In: K. C. Strassmeier and J. L. Linsky (eds), Stellar surface structure, IAU Symp. 176, Kluwer Academic Publishers, 523 Landau, L. D., Lifshitz, E. M. 1989, The classical theory of fields, Vol. 2 of Course of theoretical Physics. Pergamon Press, Oxford, England, 4th edition Lester, .I. B. 1996, “The status of continuous opacities”. In: S. J. Adelman, F. Kupka, and W. W. Weiss (eds), Model atmospheres and spectrum synthesis, ASP, 19. Conf. Series, Vol. 108 Lewis, J. L. 1957, You win again. Sun Records Ludwig, H.—C., Freytag, B., Steffen, M. 1999, “A calibration of the mixing-length for solar—type stars based on hydrodynamical simulations. 1. Methodical aspects and results for solar metallicity”, A&A 346, 111 Ludwig, H.-C., Jordan, 8., Steffen, M. 1994, “Numerical simulations of convection at 204 the surface of a. ZZ Ceti white dwarf”, A&A 284, 105 Lydon, T. J., Fox, P. A., Sofia, S. 1992, “A formulation of convection for stellar struc- ture and evolution calculations without the mixing-length theory approximations. I. Application to the sun”, ApJ 397, 701 Lydon, T. J., Fox, P. A., Sofia, S. 1993, “A formulation of convection for stellar struc- ture and evolution calculations without the mixing-length theory approximations. II. Application to a Centauri A and B”, ApJ 413, 390 Mihalas, D. 1978, Stellar atmospheres W. H. Freeman and Company, 2nd edition Mihalas, D., Dappen, W., Hummer, D. G. 1988. “The equation of state for stellar envelopes. II. Algorithm and selected results”, ApJ 331, 815 Morel, P., Provost, J., Lebreton, Y., Thévenin. F., Berthomieu, C. 2000, “Calibra— tions of alpha. Centauri A & B”, A&A 363, 675 Nahar, S. N. 2003, “Photoionization cross sections of 011, 0111, OIV and OV: benchmarking R—matrix theory and experiments”. Phys. Rev. A (submitted) Nayfonov, A., Dappen, W'. 1998, “The signature of the internal partition function in the thermodynamical quantities of the solar interior”, ApJ 499, 489 Nayfonov, A., Dappen, W., Hummer, D. G., Mihalas, D. 1999, “The MHD equation of state with post—Holtzmark microfield distributions”, ApJ 526, 451 Nordlund, A. 1982, “Numerical simulations ofthe solar granulation. I. Basic equations and methods”, A&A 107, 1 Nordlund, A. 1985, “Solar convection”, Sol. Phys. 100, 209 Nordlund, A., Dravins, D. 1990, “Stellar granulation. III. Hydrodynamic model at- mospheres”, A&A 228, 155 Nordlund, A., Spruit, H. C., Ludwig, H.—G., Trampedach, R. 1997, “Is stellar granu— lation turbulence?”, A&A 328, 229 Nordlund, A., Stein, R. F. 1990, “3—D simulations of solar and stellar convection and magnetoconvection”, Comput. Phys. Commun. 59, 119 Nordlund, A., Stein, R. F. 1991, “Granulation: Non-adiabatic patterns and shocks”. In (Cough & Toomre 1991), 141 Nordlund, A., Stein, R. F. 1997, “Stellar convection; General properties”. In (Pijpers 205 et al. 1997), 79 Nordlund, A., Stein, R. F. 2000, “3—(1 convection models: Are they compatible with I-d models?”. In: L. Szabados and D. Kurtz. (eds), ASP Conf. Ser. 203: The Impact of Large-Scale Surveys on Pulsating Star Research, IAU Coll. 176, Springer Verlag, Berlin, 362 Parkinson, W. H. 1992, “Summary of current molecular databases”. In: P. L. Smith and W. L. VViese (eds), Atomic and Molecular Data for Space Astronomy: Needs, Analysis and Availability, Vol. 407 of Lecture Notes in Physics, 2lst IAU general assembly, Springer Verlag, Berlin, 149 Pauli, WI. 1925, “Uber den zusammenhang des abschlusses der elektronengruppen im atom mit der komplexstruktur der spektren”, Zeits. f. Physik 31, 765 Pijpers, F. P., Christensen-Dalsgaard. J., Rosenthal, C. (eds) 1997. Solar convection and oscillations and their relationship, Dordrecht, Kluwer Pillet, P., van Linden van den Heuvell, H. B., Smith, W. W., Kachru, R., Tran, N. H., Gallagher, T. F. 1984, “Microwave ionization of Na Rydberg atoms”, Phys. Rev. A 30(1), 280 Plez, B., Brett, J. M., Nordlund, A. 1992, “Spherical opacity sampling model atmo— spheres for M—giants. 1. Techniques, data. and discussion”, A&A 256, 551 Riemann, J., Schlanges, M., DeWitt, H. E., Kraeft, W. D. 1995, “Equation of state of the weakly degenerate one—component plasma”, Physica A 219, 423 Rogers, F. 1981a, “Analytic electron—ion effective potentials for Z S 55”, Phys. Rev. A 23(3), 1008 Rogers, F. 1994, “Stellar plasmas”. In: G. Chabrier and E. Schatzman (eds), The Equation of State in Astrophysics, IAU Coll. 147, Cambridge University Press, 16 Rogers, F. .1. 1977, “On the compensation of bound and scattering state contributions to the partition function”, Phys. Lett. 61A, 358 Rogers, F. .1. 1981b, “Equation of state of dense, partially degenerate, reacting plas— mas”, Phys. Rev. A 24, 1531 Rogers, F. J. 1986, “Occupation numbers for reacting plasmas — The role of the Planck-Larkin partition function”, ApJ 310, 723 Rogers, F. J., Iglesias, C. A. 1992a, “Radiative atomic Rosseland mean opacity ta— bles”, ApJS 79, 507 206 Rogers, F. J., Iglesias, C. A. 1992b, “Rosseland mean opacities for variable composi- tions”, ApJS 401, 361 Rogers, F. J., Nayfonov, A. 2002, “Updated and expanded OPAL equation of state tables: Implications for helioseismology”, ApJ 576, 1064 Rogers, F. J., Swenson, F. J., Iglesias, C. A. 1996, “OPAL equation-of—state tables for astrophysical applications”, ApJ 456, 902 Rosenthal, C., Christensen~Dalsgaard, J., Nordlund, A., Stein, R. F., Trampedach, R. 1999, “Convective contributions to the frequencies of solar oscillations”, A&A 351, 689 Saha, M. N. 1921, “On a physical theory of stellar spectra”, Proc. Royal Soc. London 99(A 697), 135 Saumon, D., Chabrier, G. 1989, “A new hydrogen equation of state for low mass stars”. In: G. VVegner (ed.), \Nhite Dwarfs, IAU Coll. 114, Springer Verlag, Berlin, 300 Saumon, D., Chabrier, G. 1991, “Fluid hydrogen at high density: Pressure dissocia- tion”, Phys. Rev. A. 44, 5122 Saumon, D., Chabrier, C., Horn, H. M. V. 1995, “An equation of state for low—mass stars and giant planets”, ApJS 99, 713 Sauval, A. J., Tatum, J. B. 1984, “A set of partition functions and equilibrium constants for 300 diatomic molecules of astrophysical interest”, ApJS 56, 193 Seaton, M. 1987, “Atomic data for opacity calculations: 1. General description”, J. Phys. B 20, 6363 Seaton, M. J. (ed) 1995, The opacity project, Vol. 1. Institute of Physics Publishing Seaton, M. J., Zeippen, C. J., Tully, J. A., Pradhan, A. K., Mendoza, C., Hibbert, A., Berrington, K. A. 1992, “The opacity project — computation of atomic data”, Rev. Mex. Astron. Astrofis. 23, 19 Shibahashi, H., Noels, A., Gabriel, M. 1983, “Influence of the equations of state and of the 2 value on the solar five-minute oscillation”, A&A 123, 283 Shibahashi, H., Noels, A., Gabriel, M. 1984, “Influence of the equations of state and of the 2 value on the solar oscillations”, Mem. Soc. Astron. Ital. 55, 163 Skartlien, R. 2000, “A Multigroup Method for Radiation with Scattering in Three— 207 Dimensional Hydrodynamic Simulations”, ApJ536, 465 Slattery, W. L., Doolen, C. D., DeWitt, H. E., Slattery, W. L. 1982, “Equation of state of the one-component plasma derived from precision Monte Carlo calcula— tions”, Phys. Rev. A 26, 2255 Sneden, C., Johnson, H. R., Krupp, B. M. 1976, “A statistical method for treating molecular line opacities”, ApJ 204, 281 Stein, R. F. 1989, “Convection and waves”. In: R. J. Rutten and G. Severino (eds), Solar and stellar granulation, Kluwer Academic Publishers, 381 Stein, R. F., Nordlund, A. 1989, “Topology of convection beneath the solar surface”, ApJ 342, L95 Stein, R. F., Nordlund, A. 1998, “Simulations of solar granulation. I. General prop— erties”, ApJ 499, 914 Stein, R. F., Nordlund, A. 2003, “Radiation transfer in 3D numerical simulations”. In (Hubeny et al. 2003), 519 Stringfellow, G. S., DeWitt, H. E., Slattery, W. L. 1990, “Equation of state of the one—component plasma derived from precision Monte Carlo calculations”, Phys. Rev. A 41(2), 1105 Trac, H., Pen, U. 2003, “A primer on eulerian computational fluid dynamics for astrophysics”, PASP115, 303 Trampedach, R. 1997. “Convection in stellar atmospheres”, Master’s thesis, Aarhus University, Arhus, Denmark Trampedach. R., Asplund, M. 2003, “Radiative transfer with very few wavelengths”. In: ASP Conf. Ser. 293: 3D Stellar Evolution, 209 Trampedach, R.., Asplund, M. 2004, “Opacity sampling for 3D convection simula— tions”, A&A (in preparation) Trampedach, R.., Christensen—Dalsgaard, J., Nordlund, A., Stein, R. F. 1997, “Near— surface constraints on the structure of stellar convection zones”. In (Pijpers et al. 1997),?3 Trampedach, R.., Christensen—Dalsgaard, J., Nordlund, A., Stein, R. F. 2004a, “Im- provements to stellar structure models, based on 3D convection simulations. II. Calibrating the mixing—length”, A&A (in preparation) 208 Trampedach, R., Christensen—Dalsgaard, J., Nordlund, A., Stein, R. F. 2004b, “Im- provements to stellar structure models, based on 3D convection simulations. III. Stellar evolution with a varying mixing-length parameter”, A&A (in preparation) Trampedach, R., Dappen, W. 2004a, “MHD 2000”, ApJ (in preparation) Trampedach, R.., Dappen, W. 2004b, “Pressure ionization in the MHD equation of state”, ApJ (in preparation) Trampedach, R., Dappen, W., Baturin, V. A. 2004c, “A synoptic comparison of the MHD and the OPAL equations of state”, ApJ (accepted) Trampedach, R., Stein, R. F., Christensen-Dalsgaard, J., Nordlund,.A.1998, “Stellar evolution with a variable mixing—length parameter”. In (Giménez et al. 1998), 233 Trampedach, R., Stein, R. F., Christensen-Dalsgaard, J., Nordlund, A. 2004d, “Im— provements to stellar structure models, based on 3D convection simulations. 1. T-T relations”, A&A (in preparation) Tsuji, T. 1973, “Molecular abundances in stellar atmospheres. 11.”, A&A 23, 411 Ulrich, R. 1982, “The influence of partial ionization and scattering states on the solar interior structure”, ApJ 258, 404 Ulrich, R., Rhodes, E. 1983, “Testing solar models with global solar oscillations in the 5—minute band”, ApJ 265, 551 Vardya, M. S. 1966, “Pressure dissociation and the hydrogen molecular ion”, M.N.R.A.S. 134, 183 Vernezza, J. E., Avrett, E. H., Loeser, R. 1981, “Structure ofthe solar chromosphere. 111. Models of the EUV brightness components of the quiet sun”, ApJS 45, 635 Waech, T. C., Bernstein, R. B. 1967, “I-Iejsa”, J. Chem. Phys. 46, 4905 Willson, R. C., Hudson, H. S. 1988, “Solar luminosity variations in solar cycle 21”, Nature 332, 810 209 .. 1. ...II ...I ll. .17 . .l. l vflIbix . ...nfifi; .311. )Iubh‘v .. i \ :31 ”Ft... 5 .2... r . .5... . tween...“ 4).. . . . . , . a . . .. . . stew. .. 1 [Tu . .p . . L... 3.... : aw.