THE UNIVERSITY OF MICHIGAN INDUSTRY PROGRAM OF THE COLLEGE OF ENGINEERING THE ONSET OF THERMAL CONVECTION IN A RADIATING ATMOSPHERE Edward A. Spiegel A dissertation submitted in partial fulfillment of the requirements for the degree of Doctor of Philosophy in the University of Michigan 1958 October, 1958 IP-323

ACKNOWLEDGMENTS I should like to express my deep appreciation to Dr. M. So Uberoi whose discussions and criticisms during the whole course of this investigation have provided both guidance and incentive. I am also indebted to Dr. Leo Goldberg for considerable and thoughtful advice, and to him and the other members of the doctoral committee, Dr. L. H. Aller, Dr. A, K. Pierce, and Dr. G. E, Uhlenbeck for suggesting a number of improvements in this dissertation. There are, in addition, many others with whom I have discussed the work reported here. Of these, I should like to mention particularly Dr. G. K. Batchelor for some stimulating comments at the outset of the work; Mr. R, C, Bless for his useful discussions of the material of the last chapter; Dr. K. M. Case for some helpful remarks on variational methods; and Dr. W, Unno who has provided a number of enlightening discussions of convection theory and made many suggestions. I am grateful to the National Science Foundation for a predoctoral fellowship during the tenure of which the investigation was undertaken. I should also like to thank the Office of Naval Research for providing partial financial support during a major portion of the work, Finally, I am indebted to the Industry Program of the College of Engineering for the care with which they undertook the preparation of the final manuscript and copies of this thesis, ii

TABLE OF CONTENTS Page ACKNOWLEDGMENTS......................... ii LIST OF TABLES............................................... LIST OF FIGURES................0...... o................... vi I. INTRODUCTION........................................ 1 A. The Aim of This Investigation.............. o.. 1 B. Convection in Stars................................. 2 C. Convective Instability in a Layer of Varying Density. 4 D. The Rayleigh Theory.......................,........ 7 E. The Stability of the Radiative Gradient.............. 12 II. THE SMOOTHING OF TEMPERATURE FLUCTUATIONS BY RADIATIVE TRANSFER..,O, o...o.................. 14 A. Preliminary Remarks..........o................... 14 B. The Radiative Heat Equation.................... 14 C. The Linearized Equation....0.......................... 18 D. The Homogeneous Medium................ o............ 21 E. Spherically Symmetric Disturbances............o.. 6 25 F. Some Numerical Values.,,...........o................. 28 III. THE EQUATIONS OF THE PROBLEM........o............... 30 A. The Basic Equations......................o 30 B. The Equilibrium-State Equations............o... 31 C. The Perturbation Equations........e............o o32 D, The Equations in G and w Alone.................. o 35 E. The One-Dimensional Equations...................... 37 F. The Boundary Conditions............................. 38 IV. THE ONSET OF CONVECTION.............................. 40 A. The Equation of Marginal Stability................... 40 B. A Variational Principle............................ 41 C. The Constant Gradient Approximation.......o......... 44 D. Calculation of Xcrit............................... 46 E. Physical Interpretation....... o................... 54 V. CONVECTIVE INSTABILITY IN EARLY STARS................... 60 A. Plan of This Chapter............................. 60 B. Expression of X for Application to Stellar Atmospheres.............................. 61 C. Convective Instability in the B Stars................ 63 iii

TABLE OF CONTENTS (CONTtD) Page Do The Effectiveness of Convection in the Early B Stars................................4... 69 E. Instability in Early A Stars................. 76 F. Summary of the Investigation........................ 87 APPENDIX I THE PRINCIPLE OF THE EXCHANGE OF STABILITIES....... 90 APPENDIX II THE KERNEL E o....................... 95 APPENDIX III REDUCTION TO THE CLASSICAL THEORY*....*... *........ 98 APPENDIX IV THE CALCULATION OF k1....... eoo4..o*......... 103 APPENDIX V THE CALCULATION OF X2.....1....................... 109 APPENDIX VI TABULATION OF U(a,.) AS CALCULATED BY A FIVE-POINT RADAU FORMULAo o... o*. o... o 4. o..4. 114 BIBLIOGRAPHY............................ 115 iv

LIST OF TABLES Page I Critical Rayleigh Numbers and Modes of Maximum Instability............................................ 11 II Characteristic Times of the Solar Photosphere............ 28 III Scale of Disturbance with 4-Minute Lifetime.............. 29 IV Values of Xcrit and acrit...................... 52 V Stability Criterion Quantities for Saito's Model with lI' = 15,500~K..................................... 67 VI d and T for the Helium-Free Models....................... 80 v

LIST OF FIGURES Figure Page 1. Minimum Eigenvalue of R Vs. k (Schematic)............. 10 2. Decay of Central Temperature in a Spherically Symmetric Disturbance..................... 27 3. Illustrating the Trial Functions..................... 49 4. Xcrit Vs. T..... o o........... o o......... 55 5. acrit Vs. eT o*eoe*......o*............................. 56 6. A Vs. T in Underhill's Model.............. o....... o... 65 7. A Vs. z for Saitots Model with log g = 3.8 and and.ie = 15,500~.........o.o......o....e............. 68 8. A Vs. T for TC = 10,700~0..... oo..o.... o 77 9. A Vs. z for I', = 10,7000.~ o.................. 78 10. A Vs. T for Osawats Models (Te = 8,900 ).............. 81 11. A Vs. z for Osawats Models ('i-f= 8,900)............. 82 vi

I. INTRODUCTION A. The Aim of This Investigation If viscosity and heat conduction are negligible in a stellar atmosphere, then whenever the temperature gradient exceeds the adiabatic gradient, convection can occur. The character of the flow pattern, and particularly the scale of motion, will depend on the superadiabatic excess and its space variation. For reasons we will discuss later, the convective pattern will tend to develop in cells of zero horizontal dimension so that the motion would be entirely vertical. However, when the effects of viscosity and heat conduction are introduced, the situation is altered. A certain minimum excess over the adiabatic temperature gradient is required before convective instability becomes possible, and even when this minimum excess is achieved, the pattern of motion is not the same as in the inviscid, nonconducting case. The most striking feature of the altered convective instability is that, although cellular motion is again expected, there is one particular cell dimension, or mode of maximum instability, which is preferred. These qualitative statements describe experimental and theoretical results of long standing. The occurence of cellular convection when a certain critical temperature gradient is exceeded was detected in the experiments of Benard (1900). In a famous paper stimulated by Benard's experiments, Rayleigh (1916) showed that a mode of maximum instability and a critical gradient were to be expected on theoretical grounds. Even earlier, Rayleigh (1883) and Lamb (1890) had discussed the convective instability of inviscid, nonconducting media. Yet, the -1

,2astrophysical importance of these studies was not appreciated until Unsold (1930) showed that temperature gradients in stars can be superadiabatic in regions where hydrogen or other abundant elements are being ionized. Since that time, much attention has focussed on problems of convective instability, especially those involving constraints found in stars and planets. In spite of this increased interest in convective instability, the nature of the heat transport mechanism has not been widely discussed, and most of the work has followed Rayleigh's (1916) in including only thermal conduction. But in stellar atmospheres, radiative transfer is by far the dominant mechanism of heat transport both because it is a volumetric effect and because the radiative conductivity considerably exceeds the thermal conductivity [Eddington, (930), Chapter V; Edmonds, (1957)]. The aim of the present work is, therefore, to study the modifications of the classical stability theory that are required by a radiating medium. The usual classical approximations are retained; most notably it is assumed that the superadiabatic layer does not encompass large variations in temperature and pressure. This restricts the direct application of the results to early-type stars, though, it is hoped that the insights gained may lead to approximations which will permit the extension of the work to include the conditions demanded by the later spectral types. B. Convection in Stars Observational evidence for the presence of convection in stars comes mainly from the absorption-line profiles of stellar spectra.

-3Doppler broadening, too great to be explained by thermal motions of atoms, indicates the existence of large-scale motions in a variety of stars. These motions have generally been referred to as "turbulence" in the astronomical literature, and probably the convection is turbulent in the late spectral types. For example, in the sun, the r.m.s. vertical velocity as found from Doppler shifts of "granules" is about.4 km/sec [Richardson and Schwarzschild, (1950)], and with reasonable estimates of viscosity [Edmonds (1957)] and scale height, this implies a Reynolds number of 1011. Attempts to produce models of stars with convection have generally been based on some idealized picture of the convective processes. No complete theory has been produced and each model has left at least one parameter to be fitted to the observations. The most elaborate approach has been the application to the solar convection zone of the mixing-length theory of turbulence initiated by Biermann (1937, 1948) and Siedentopf (1932). Nowhere in the theory is the mixing length specified, though in her model of the solar convection zone, Vitense (1953) uses the scale height as a mixing length. Perhaps the most promising alternative approach lies in the study of nonlinear convective instability begun by Malkus (1954 a) and developed further by Stuart (1958) and Malkus and Veronis (1958). Malkus' (1954 a) theory has succeeded in predicting what seem to be transitions from laminar to turbulent convection. Such transitions have been observed by Schmidt and Saunders (1938) and further investigated experimentally by Malkus (1954 b).

-4The starting point of the nonlinear theories is a complete understanding of the problem of instability with respect to perturbations of small amplitude. The main concern of the work done so far has been with problems in terrestrial contexts, and they naturally develop from the Rayleigh-type convection. If the techniques used are to be extended to stellar problems, a complete solution of the linear stability problem under astrophysical conditions must be found. These conditions include (1) large variations in density and temperature with height in the atmosphere, (2) the bounding of an unstable layer by stably-stratified gases which however may be penetrated by motions, and (3) the smoothing of temperature fluctuations by radiative transfer. Studies of the effect of these conditions on the motion have been attempted where only one of the conditions is introduced at a time, but the complete solution of even the linear problem with. all three effects present is still awaited. In the next section we make reference to the work on the first two effects. C. Convective Instability in a Layer of Varying Density The study of convective instability began with the work of Rayleigh (1883) and Lamb (1890)o The interest at the time was primarily in water and atmospheric waves, and as we have mentioned, the relevance of such studies to problems of stellar atmospheres was not appreciated until Unsoldts pioneering work. In this early work the effects of viscosity and heat conduction were ignored. A static fluid, stratified in horizontal plane-parallel

-5layers was subjected to a uniform, vertical gravity field. Small, arbitrary perturbations were introduced in the state variables and these were coupled through the equations of motion to a flow pattern of small amplitude. The critical situations of zero density gradient (imcompressible fluid) and zero superadiabatic gradient (compressible) were assumed to be exceeded. All disturbances grew in time since no stabilizing influences were included. The aim of the investigations was to see which disturbances developed most rapidly. Only disturbances periodic in horizontal planes were considered. Thus, if Cartesian coordinates are used, and x and y are the horizontal axes, the disturbances depend on x and y like exp(ikxx + ikyy) where the kIs are wave numbers. A given mode of disturbance is characterized by the horizontal wave number k = (k + k2) / It was assumed that the time dependence of a disturbance is like ent, where n must depend on k. In all cases, for both compressible and incompressible flow, and for a number of initial configurations, the rate of development, n, is zero for k = 0 and increases monotonically with increasing k. Thus, disturbances of small horizontal dimension develop most quickly. The increase of n with increasing k can be understood intuitively; we will here follow a brief discussion given by Hide (1955). Hide first recalls that a given static density distribution is unstable if it is not a minimum potential energy configuration. The tendency of an unstable configuration, he argues, will be to convert potential energy

-6into kinetic energy in the most efficient way possible. Since only vertical motions are effective in this conversion, a disturbance consisting entirely of vertical motions should develop most quickly. This is the periodic disturbance of zero horizontal wave length (infinite k) and it is, therefore, the preferred mode of instability. Similarly, the disturbance of infinite horizontal wave length consists entirely of horizontal motion and should be entirely stabilized, as the theory predicts. The details of the theory are to be found in Lambbs book (1931; see especially page 541 ff); Lamb (1890) considered the compressible flow and Rayleigh (1883) the incompressible. Some numerical results, especially a graph of n versus k, were recently derived by Skumanich (1955) for a polytropic atmosphere. Skumanich stressed the result that, in a polytropic atmosphere, n - oo as k - oo. Another aspect of convective instability in atmospheres of varying density is the penetration of the motion into the ambient stable layers. Unno (1957) has considered this effect in an atmosphere of exponentially varying density in which the (constant) scale height changes discontinuously at one depth in the atmosphere. Here again n increases with k. Further work on this aspect of the problem has been done by Gribov and Gurevich (1957)*, but apparently the only readily accessible reference is an abstract. * I am indebted to Dr. W. V. R. Malkus for drawing this work to my attention.

-7The situation changes entirely when the effects of viscosity and heat transport are admitted. Each of these favors stability, and each increases in effectiveness with increasing k. The case in which viscosity alone enters the picture has been considered in connection with variable density and incompressible flow by several people at Yerkes Observatory [Chandrasekhar (1955), Hide (1955), Fan (1955)]. The situation is that infinite wave number implies infinite viscous stress, and the result is found that the greatest rate of development occurs for a finite k, depending on the initial configuration and the parameters of the problem. Even in this case the critical gradients for the onset of convection remain the adiabatic and zero gradients for compressible and incompressible flow respectively. It is not until both viscosity and some form of heat transport act that an adverse density gradient can be stable. D. The Rayleigh Theory If a density fluctuation appears in a medium with an adverse density gradient, the disturbed region of fluid will be driven vertically by the resulting net buoyancy force. Viscous forces will oppose the motion, but so long as nothing acts to wipe out the fluctuation, the buoyancy and viscous forces will balance and a sustained motion will result. This situation, which was discussed at the end of the preceding section, bears a similarity to the case of a sphere falling through a viscous fluid with the Stokes terminal velocity. However, in real fluids heat conduction will destroy the concomitant temperature fluctuation and the longevity of the buoyancy force is curtailed. The question

-8- in any given case is then, does a, sustained motion arise or are the fluctuations destroyed so quickly that convection cannot be maintained. The answer, of course, must depend on the initial temperature gradient and the nature of the fluctuation, as well as on the hydrodynamic parameters of the fluid. The familiar heat transport mechanism, thermal conductivity, was first introduced into stability problems by Rayleigh (1916). This complication of the theory required the introduction of equations other than the equations of motion, since an additional variable, the temperature, now appearedo A heat conductivity equation and an appropriate equation of state were therefore introduced into the problem by Rayleigh (l916). [The set of equations was not original with Rayleigh, and he credits them to Boussinesq (1903).] To compensate for the new difficulties, Rayleigh made use of certain very natural approximations, which he attributed to Boussinesq (1903), of which only one is a severe limitation. This approximation treats all variations of density as negligible except when they appear in conjunction with gravity. Thus, any derivative of the static density distribution is neglected when multiplied by a perturbation quantity if the term involved does not contain a gravity factor. Hence, Rayleigh's work (and the present one as well) does not apply to layers of fluid with large density variations from top to bottom. Rayleigh also treated the motion as incompressible so that both the equations of continuity and state become simplified. However, Jeffreys (1930) has shown that so long as the density does not vary greatly across the layer, the results obtained for incompressible

-9motion can be transformed into exact results for the compressible case by everywhere replacing the temperature gradient by its excess over the adiabatic gradient. As in the previous investigations, Rayleigh sought solutions of the perturbation equations which depend on x, y and t like exp(ikxx + ikyy + nt). He assumed n to be real and set n = 0 to obtain the critical equation of marginal stability. Thereby, he found a single equation for the vertical velocity component alone. The equation is a sixth order ordinary differential equation involving the horizontal wave lumber k = (kx + ky) number k = (kx + ky)!/2 and a non-dimensional parameter known as the Rayleigh number, Here, g is the acceleration of gravity, a the compressibility, B is the static temperature gradient, d the layer thickness, v is the kinematic viscosity, and where K is the thermal conductivity, p is the density, and c is the specific heat per gram. All these quantities are constants across the layer, The problem is to find the minimum value of R, for given k, which will permit solutions of the differential equation. The parameter R, it may be added, is interpretable as a ratio of buoyancy force to viscous force. Since the work of Rayleigh, a number of other papers have appeared, extending Rayleigh's results to a greater variety of boundary conditions [Jeffreys (1926), Low (1929)], and introducing some refinements,

-10We will summarize briefly the outcome of these papers which culminated in the comprehensive discussion by Pellew and Southwell (1940). The mathematical aspects are omitted here since they are essentially contained in the later chapters and are also concisely summarized by Prandtl (1952). For a fixed horizontal wave number it is possible to find solutions of the equation of marginal stability. Solutions do not exist for all values of R; there is a minimum allowed value of R for each horizontal wave number depending on whether the boundaries are free surfaces or rigid walls. Figure 1 shows schematically the sort of result obtained for a particular type of boundary; it shows the locus of minimum eigenvalues of R. The curve has a minimum value, Rcrit, which occurs for a value of k, kcrit, known as the mode of maximum instability. 1rit - ~- I - T / Rcrit H -i - - crit Figure 1. Minimum Eigenvalue of R Vs. k (Schematic)

-11The solutions corresponding to the minimum eigenvalue are the so-called fundamental modes; they are the ones which do not have zeros inside the fluid layer. The results of the various investigations are summarized in the table, Wave numbers are expressed in the dimensionless form a = kd, and acrit gives the mode of maximum instability. TABLE I CRITICAL RAYLEIGH NTMBERS AND MODES OF MAXIMUM INSTABILITY Nature of Boundaries acrit Rcrit Two Free Surfaces 2.22 657o7 Two Rigid Surfaces 3513 1707.8 One Free, One Rigid 5o36 1100o7 The conclusion is that a layer of fluid is stable against thermal convection unless R > Rcrito However, if the critical Rayleigh number is exceeded, disturbances in a band of wave numbers are unstable, and the most unstable mode is that corresponding to acrito Pellew and Southwell (1940) showed that a variational principle can be found for the problem which gives accurate approximations to the minimum allowed values of R. Their method gave impetus to later researches on problems of convective instability. They also demonstrated that, if n is complex, its real part is negative. Hence, overstability cannot occur and setting n = 0 to find the condition for marginal stability is valid, Previously the justification for setting n = 0 in this way, the so-called principle of the exchange of stabilities, was appeal to experiment, We will see in Appendix I that the principle is also valid in the radiative case,

-12E. The Stability of the Radiative Gradient We have mentioned that radiative transfer can play a role in stabilizing a stellar atmosphere by effectively equalizing temperature differences. Another aspect of convective instability in a radiating medium is that the initial static temperature gradient is not a constant, as in the case when the static gradient is produced by the action of thermal diffusion alone. Thus the problem arises, how is the Rayleigh theory to be modified to allow for a variable steady-state temperature gradient. The problem of a variable temperature gradient has been discussed by Malkus (1954 a) who noted that the form of Rayleigh's differential equation is essentially unaltered. The Rayleigh number enters as before, except that the temperature gradient is replaced by its average value. Malkus has shown that, for the case of free bounding surfaces, a variational principle exists for the calculation of minimum allowed values of the mean Rayleigh number. The effects of variation in the temperature gradient are thus included in the mean measure of stability. Using Malkus variational principle, Goody (1956) has investigated the stability of static gradients arising frg9 tt ombfeinte action of radiative transfer and thermal diffusion. The static gradients are calculated on the basis of the Eddington approximation. To describe the flow of heat by radiative transfer Goody uses two alternative approximations: (a) the opaque approximation in which radiative transfer is treated as a diffusion process; (b) the transparent approximation in which essentially Newton s law of radiation is employed and transfer

>13 effects are neglected. For several values of the ratio of radiative to thermal conductivity, numerical calculations of the mean critical Rayleigh number are performed. The results depend upon the conductivity ratio and the optical thickness of the layer. Goody's results apply strictly only to very large and very small optical thicknesses, though he graphically interpolates to obtain estimates for the intermediate cases. No attempt is included to subtract the adiabatic gradient out, and application to a compressible medium would require a separate numerical calculation.

II. THE SMOOTHING OF TEMPERATURE FLUCTUATIONS BY RADIATIVE TRANSFER* A. Preliminary Remarks In the last chapter we saw that an atmosphere which is able to wipe out temperature differences effectively will tend to be stable against thermal convection. The equation describing the smoothing of temperature fluctuations by radiative transfer is developed in this chapter. One of the by-products of the development is an expression for lifetimes of such fluctuations under the action of radiative transfer alone. Radiative lifetimes are of interest in a variety of problems of stellar hydrodynamics, and the processes considered have received attention in the astrophysical literature. Some numerical estimates are given by Whitney (1955) for radiative lifetimes in stellar pulsation and by Oster (1958) for radiative relaxation times in solar flares. Approximate formulae were developed by Unsold [cf. Unsold (1955)] and by Vitense (1953) for studies of the solar convection, but these are for optically thick and thin perturbations. An equation for the smoothing of fluctuations in the chromosphere was developed by Parker (1953) which also neglected transfer effects. The trend of these investigations indicates a need for a formalism to incorporate effects of radiative transfer into problems of stellar hydrodynamics. This need is partly met in the work of the present chapter. B. The Radiative Heat Equation We seek to express an equation for the time variation of temperature in a radiating medium. The basis for this equation will be * Much of this chapter has appeared elsewhere [Spiegel (1957)]. -14

-15the statement of the conservation of energy for the medium. Changes will be supposed to proceed quasi-statically, so that a temperature distribution, source function, and mean absorption coefficient can be defined for the medium at any instant. Let t be a position vector, t be the time, and let e(3;t) = internal energy per unit mass, p(i;;t) = matter density, T(5; t) = temperature, p (i;t) = pressure, Q(3-;t) = radiative heat exchange rate per unit volume, K(];t) = thermal conductivity, and (5t;t) = rate of viscous dissipation of energy per unit volume. Then we have [e.g., see Goldstein (1938); Chapter XIV] o + t) = (KT) Q as the expression of conservation of energy where &-t at x8= at The symbol ui denotes the velocity vector with components (u,v,w). The viscous dissipation term is of second order in the velocity and will be negligible in all problems considered herein. Likewise, thermal conduction is not important in most stellar problems [Eddington (1930), cf. Chapter V], and though the retention of the thermal conduction term would entail no special difficulty, we neglect it. For

-16example, Edmonds (1957) has noted that the thermal conductivity is always less than or equal to 5% of the radiative conductivity for the physical conditions encountered in Vitense's (1953) model of the solar convection zone. For the internal energy, we have the relation de = c4rT (2) where cv is the specific heat per gram at constant volume. The net rate of radiative energy exchanged per unit volume with the surroundings is rate of absorption minus rate of emission. If 3B(;t) = source function integrated over frequency, J(;t) = mean intensity integrated over frequency, ^(2;t) = mean (grey-body) absorption coefficient per cm., where K is an inverse length, and J and B are in units of intensity (i.e., ergs/cm2/sec/sterad). Then we may write @Q = (J -B)J (3) The mean intensity is an average of contributions from over the entire medium, suitably weighted according to the extinction. The extinction of a beam between two points depends exponentially on the optical separation between them. We define the optical separation between two points x and J as ( X; t) -J f( ( t) &A (4) Jo(h

-17where X and. Then the mean intensity can be expressed J I! tl= Kf tB t i t) 4( where d & is a volume element and the integration is over all space. The transit time of light has been neglected. To proceed further we will assume local thermodynamic equilibrium and thus take the Planck function for the source function. Hence B 7iC T_ r(6) where a is the Stefan-Boltzmann constant. Finally, if we assume incompressibility we may neglect the PP d (-) term, though in the absence of motion the effect of compressidt p bility can be included by replacing cv by c. Now on combining Equations (1), (2), (3), (5) and (6) we obtain CV IA:t - 4- 0- K_ which will be one of the fundamental equations of this investigation.

-18C. The Linearized Equation Let the temperature be expressed as T(,t), = T t), e (Tt't) (8) where To is some equilibrium distribution and Q is the temperature fluctuation. In the convective instability problem which we will consider in later chapters, the assumption is made that any perturbations in equilibrium quantities are small. Actually, for the case of stellar temperatures this assumption still allows reasonably large fluctuations in temperature, so that even in the relaxation-time study the approximation 1| < < To is of interest. We can then neglect terms of second order in @. Since we have assumed incompressibility, one variable suffices to specify the state of the system. Thus, if we choose the temperature as the state variable, any other physical parameter of interest can be expressed in terms of the temperature explicitly, and in terms of x and t implicitly. If we call f a typical quantity, we may write the Taylor expansion, to first order in Q, as F(7T) 7 C(To) +( )T e (9) In the sequel, we will use the notation fo = f(To) and f = (r)T o 0 UTT=To and write the t and t dependence explicitly in the following way: f(~;t) = fo(3) + f$(t) Q (f,t) (10)

-19Then, if there are no motions in the equilibrium state, Equation (7) reduces to To ( ) = To( ) 4- ) () which is the condition of radiative equilibrium. More generally, on combining Equations (7), (8), (10) and (ll), and neglecting o(2) we get \P" Lt tcL @(V ~t)4\ (T5) ) — {-)-.\(- ytk'V) eo L + ) Ltr7[, ( x) 9 (C;t) t (12) 4t/' t, ( X )r W(xt ) \9)(x^n?;t)c ]Th)?, where - n the last integration, let us replace the where In the last integration, let us replace theby variableR by a= h/n

-20Then Equation (12) becomes,'- 0',-t + u.v (-( 9 ( l J4te + ( X \a l /o l X In this linear formulation only the effects of perturbations on the (13) 4Toi )e^ 9 source function and absorption coefficient appear. The first term on e terms on the right hand side of Equation (13) represent the change in energyious physical origins of the radiative smoothing of temperature fluctuations. In this linear formulation only the effects of perturbations on the source function and absorption coefficient appear. The first term on the right hand side of Equation (13) represents the change in energy absorbed at x as a result of the alteration of source function throughout the atmosphere. In the second term is contained the alteration in the rate of emission at 3? which results from the disturbance in the source function there. The third term arises from the perturbation in the absorption coefficient throughout the atmosphere and gives the resultant change in energy absorbed at 3t. The absorption-coefficient perturbation term consists of two opposing tendencies, since a change in absorption coefficient produces a change in both the absorption and

-21emission. We will see that these opposing tendencies cancel exactly in the isotropic radiation field of a homogeneous medium. In order to obtain the relaxation time for the radiative smoothing of the fluctuation 0 we will discard the convective terms in Equation (13). Even the remaining terms are too complicated to deal with in the general case, and we turn to an examination of the equation for homogeneous media. D. The Homogeneous Medium For an infinite homogeneous medium with no motions, Equation (13) becomes C -;/Ff r e'.) */{r e((2tt)) + fat L1O(Rj itJ f i J e? 13 ()d el' r -j-K (14) where - _ l^ Tl (15) 1 f cv is an inverse time, and the subscript denoting equilibrium values has been suppressed. Let us first consider the absorption-coefficient perturbation term p - JcJ[en(X+^ t)- n (X +,l e A? (16) e IkQ ec- ir -- ('-rr

-22In the second term of ~, we will use a as an integration variable instead of r. However, we will call the new integration variable r to retain consistency with the first term. Then f J *HL'l _ 1 - — *" *:,-,t l,IP-41(17) ij - ( i ia;t) L- i Lb If we now interchange the order of integration and perform the integration over A, we find that (P) ~ 03 (18) Hence, the perturbation of absorption coefficient has no net effect on the smoothing process in a homogeneous medium. This is true only in this linear approximation; the linear theory neglects the interaction of absorption-coefficient perturbation at a point with the perturbation in source function at other points. With the vanishing of it, Equation (17) reduces to If we take te F e,trnsor -of thi 1el9) If we take the Fourier transform of this equation, we obtain at ~}(t;t}= - (v ) < t t)) (20) D t k )~~~~~~~~)"

-23= where ~4)t1) 7(WQ31(x ){)c (21) and -X m / i ( fit e - it d —) p(22) Then ( decays exponentially in time according to the equation -,ot) = i(t,);o) e (23) where the initial disturbance specified by ((kS;O) is of arbitrary form. The integral term in n, -can be evaluated using spherical coordintes. If i, we ca can be evaluated using spherical coordinates. If L = cos( ), we can write After integrating over l, we have i-t) = (26) ~() \ E- cT~ ~ Soi

-24Hence, (I1i (~0 i -Kt (8 O' ^ t(27) Then, n becomes (28) (10 = X ( f -_'< (8 The quantity 1/n is the approximate decay time for a perturbation of characteristic length 1/k. The ratio K/k is a measure of the optical thickness of the perturbation. As K/k varies from zero to infinity, n(k) varies monotonically from 7 to 0. In the former limit, when the perturbation is optically thin, only local emission operates to smooth out the fluctuations. Wherever the temperature is raised, the rate of emission is increased; where it is lowered, there is a decreased emission rate. But no self-reversal occurs, and any given point does not, in a manner of speaking, see neighboring points. In this way, the smoothing rate becomes independent of k. At the other extreme, large optical thickness, K/k - co and n(k) -> (- )k2. The dependence of n on k is just as in ordinary thermal 3K2 conductivity and 7/3K2 has the same dimensions as the conduction coefficient in the usual heat equation. In this limit, a photon travels several mean free paths before escaping from a perturbed region, and the smoothing rate is no longer volumetric. In general, if the initial perturbation is Q(x;0), (P(I o) — e(., )o( ) e 0) - (29)

-25and the solution of Equation (19) is obtained by taking the inverse transform of Equation (23). We then find that *. d (f r p ^ -^'*?'': ^':) (30) is the general solution for small disturbances in a homogeneous medium. E. Spherically Symmetric Disturbances As an example of the sort of calculation that might arise, let us consider initial disturbances which are symmetric about the origin. The general solution, Equation (30), becomes ~14*K^A )-, ^ -' an V'1 0a P. ^ (31) where, |. We can integrate over all angles in the i integral and obtain n-^ f. a x GC)t) ^.J e'' (;, [ )-,o).": (32) In turn, the ki integration can be performed over all angles and we find that the distribution remains spherically symmetric at all later times.

-26Thus, e tcX ) -) j A kJ; \ Ra4 4',Xo0) e-' O i',/'Ib'qV ",'"' (33) where x = x | In particular, let us consider ('),0) =-.5 e~ where Q0 = G(0;0) and a are constants. We find, then, that uee 1 _k ^kx ^(-^y (34) The integral in Expression (34) is difficult to perform in general and must be done numerically. For the present discussion we will restrict ourselves to an examination of the time dependence of the temperature at the center of a disturbance with K2 = 1. Then, if we set x = 0 and K1 = 1 in the integral of Equation (34), we obtain',, J +, - -)2- (3 l) ^ \ ~ Riyr M

-27With the transformation ^ ^= co)0 /' 8(36) Equation (35) becomes u o "it _<t(l-^ ) 9o )' ) e Oi C,;3 0 r (37) This integral expression may be readily evaluated by Gaussian quadratures; a four-point quadrature leads to the result presented in Q(O;t) Figure 2. The solid curve is a plot of -;0 against the dimensionless time, yt. We see that the decay is somewhat slower than would result from the simple exponential decay curve, e-t, represented by the broken line. 1.0 \ \.8.6 I4- \ Q(o;t) ~~~~\''^^*-~ — -_ @~Q(o;o).2 e-7t 1 2 3 4 5 6 7 8 9 10 At Figure 2. Decay of Central Temperature in a Spherically Symmetric Disturbance

-28It is interesting to note that the rate of cooling at the center of the disturbance decreases in time. This is probably due in part to the spreading out of the disturbance as time passes, with the result that the center of the disturbance becomes optically shielded from undisturbed regions. F. Some Numerical Values To obtain some idea of the relative importance of radiative effects, let us briefly examine the situation in the solar photosphere. Using the model given by Allen (1955), and cv from Biermann's (1942) tables, we find the results in Table II. TABLE II CHARACTERISTIC TIMES OF THE SOLAR PHOTOSPHERE Mean optical depth 0.5 1.0 1/Y (seconds) 2.5 1.7 The value 1/7 is the limiting lifetime of optically thin perturbations and represents the lower limit of possible lifetimes. These are small compared to observed lifetimes of about 4 minutes for features in the solar photosphere. However, with these values of 1/y, we can ask what optical thickness of perturbation would lead to a 4-minute lifetime in Equation (28). For the uppermost visible layers, the following results are found.

-29TABLE III SCALE OF DISTURBANCE WITH 4-MINUTE LIFETIME Mean optical depth 0.5 1.0 K/k 5.5 7.0 1/k (km) 500 350 These values of 1/k are of the order of magnitude of the so-called granule size [cf. Rogerson (1958)] although a number of geometrical factors are omitted, and we may expect radiative effects to be of importance in stellar hydrodynamics.

III. THE EQUATIONS OF THE PROBLEM A. The Basic Equations Consider an atmosphere in a uniform vertical field of gravitational acceleration, g. Let us use a Cartesian coordinate system with z as the vertical coordinate measured positively upwards. The equations of motion are "fd D — ^ - -k 5t j, (38) at = - 4 fii d~r (39) jcA^ _^ 4; -. va _ (40) where v is the kinematic viscosity and u, v, and w are the velocity components. We will restrict ourselves to the Boussinesq approximations discussed in Section D of Chapter I; in this way we can treat the motion as incompressible, and still have the possibility of converting to compressible results by means of the Jeffreys theorem discussed in that section. Then the equation of continuity is?u ~'anr- 3-ur,- -- + -" = 0 (41) -30

-31Since the fluid is treated as incompressible, we require only one variable to specify the state of the system. We may, therefore, take as an equation of state, As = [I-~<(T-T J] (42) where a is the compressibility and p0 is the density corresponding to some reference temperature To. Our final equation is the radiative heat Equation (7) developed in the preceding chapter,.Cdt = A JET (7) - Here we have written c for the specific heat with the understanding that when the motion is incompressible cv is implied, but cp should be used when extending the results to compressible fluids by the Jeffreys theorem. The six equations, Equations (38) - (42) and (7), for the six unknowns p,T,p,u,v,w form the basis of this investigation. B. The Equilibrium-State Equations Let us assume that the basic equations admit an equilibrium solution with no motions. If we then set all time derivatives equal to zero and let the velocity vanish, we obtain from Equations (38), (39) and (40), 9; %~ (43) - = - (44)'p -() (45)

-32Equations (43) and (44) imply that the steady-state pressure depends on z alone; i.e., the atmosphere must be plane-parallel in the equilibrium state. Equation (45) is the well-known condition of hydrostatic equilibrium. The radiative equilibrium condition obtained from Equation (7) we found to be Expression (11) T4 = J KT Lc A(11) We will denote solutions of the equilibrium equations by pO, po, and Toy where these three quantities depend only on z. It is taken for granted that the equilibrium solution can be obtained. Equations (45) and (11) are the basic equations of model atmosphere theory for grey atmospheres in local thermodynamic equilibrium. C. The Perturbation Equations'The question to be investigated is whether an atmosphere described by the equilibrium solution po, To, po will be restored to equilibrium if the medium is slightly displaced from the equilibrium configuration. To answer this question, we consider a configuration described by Tr=~~~~~ ~~T0,~~+6 (46) I = T + e (46) and ask whether the initially small quantities 5p, Q, 6p will grow or diminish in time. The velocity in this disturbed state, (u, v, w), is

-33also supposed to be small. We will then introduce solutions (46) into the basic equations and neglect all second order terms in small quantities. The equations of motion, Equations (38) (40), become LLI aat P x tZ, (47) J^Set - -t -ft <1- P'^7 Y (48) j~4, =... _ f v r - iSt pQ?7?iT C~ ~ f $a~-sp(49) tin which we have made use of te em e, Es in which we have made use of the equilibrium equations, Equations (43), (44), (45), and (11). The equation of continuity, Equation (41), remains unaltered: gt +be g "- 0 -t~i a-rt -— O (41) The equation of state, Equation (42), becomes i\ =z - x< Po 6 (50) In Chapter II, we derived the linearized heat equation, Equation (12). In that equation, let us introduce the notation a & T~7 (51) cvzX

-34where it is understood that when we speak of the applications to compressible media, this definition is to be replaced by lUTo (UTn dTo cT)D (52) in which ()A is the adiabatic gradient. In either case, we know that dz AD an atmosphere is completely stable unless P<0 for at least some values of z. We will here consider only atmospheres having some horizontal layer, of thickness d, in which P may be negative. We assume that ( vanishes on the bounding planes of this layer and that P is positive definite outside the layer. The origin of coordinates will be taken such that z = 0 on the plane midway between the two bounding planes. Whether any particular combination of negative P, and d will be adequate to produce convection is still a matter to be investigated. In any case, for IzL>d there will be stabilization and the disturbance quantities will go to zero as Izl ->. We have already assumed that the density changes across the layer Izl ~ are not appreciable. Consistent with this approximation, we will assume that To and Ko also do not vary much in Iz14. Then since the integrands of the terms on the right hand side of Equation (12) become small for Izl|>, variation in To and Ko outside the (possibly) unstable layer will not cause much change in the values of the integrals of Equation (12). Hence, Equation (12) becomes -_ i ( t ) (53),^1 ry[ 9 -Q) al'i-YK VJn'

-35where, as in Section D of Chapter II, the absorption-coefficient perturbation term has vanished identically. Equations (47), (48), (49), (41), (50), and (53) form the set of six equations governing the six disturbance quantities 5p,G,6p,u,v, and w. It will be our next step to reduce the number of equations. D. The Equations in Q and w Alone Ultimately, we will reduce the six equations derived in the previous section to a single equation in one unknown. As we shall see, it is also convenient to carry out the reduction to two equations in two unknowns as an intermediate step. Let us differentiate Equations (47) and (48) with respect to x and y respectively. In these equations, we may treat p0 as constant, as implied in the Boussinesq approximation introduced above. On adding the results of the differentiation we obtain a (1^ ^ ^.}- - f i " (54) f^( ai a o ^ where 7. - aX " (55) and the equilibrium subscripts are hereafter omitted. We can then combine Equations (54) and (41) to obtain Vt ~. ^ +? (56) Pa~~~~~~~~~~~~z~~C t~r.

636If we eliminate bp between Equations (49) and (50) and operate on the result with A7 we find that l'jt-,r= -rev b -VI>r + t me We can eliminate 6p by differentiating Equation (56) with respect to z and adding the result to Equation (57). Then ~-tV. TV gw)+ 1%6 ~ <09 (58) together with Equation (53), — "~L-P"UIU~e ~ o j e(53) { r er-, d) provide two equations in the two unknowns ~ and w. At this point, the contrast between the present work and the classical theory is made clear. In the classical case, in the reduction to two equations in Q and w, one of the equations is our present Equation (58). The other equation would differ from Equation (53) primarily in having the Laplacian of ~ on the right hand side instead of the terms which are there. Also, in the spirit of Malkus' (1954 a) extension of the classical work, we may mention that 3 in Equation (53) may vary with z according to the equations of radiative equilibrium. It should be emphasized that the use of the Boussinesq approximations appeared frequently in the operations leading to Equation (58).

-37E. The One-Dimensional Equations The first step in dealing with Equations (53) and (58) will be to reduce the partial derivatives to total derivatives by seeking separable solutions. In particular, let us take Xir(X'e) = e W% ( ) (59) and (et ).XX "11li^@~ ( ) (60) where kx, ky and n are constants. We see that in the space dependence of these solutions, we have merely restricted attention to particular harmonic components in the horizontal plane, Since the equations are linear, any linear combination of these will also be a solution. With these forms of the solutions, Equations (53) and (58) become the one-dimensional equations ('. 7 )-)V (- )vT - (61) and n ) ( ) (r) (62) R X - V ^ i d/ 1r n" i_'-,: >)

-38where c - i. a )' d~ ~~~~~~~~ ^ ~~(63) and G. -Xvs-)/ J J.4 () J J (64) The two coupled equations, Equations (61) and (62), in E and W, together with the boundary conditions stated in the next section, will be used in Appendix I to demonstrate the principle of the exchange of stabilities. F. The Boundary Conditions In the regions of stable stratification the perturbations will be damped out, and for great enough distances from z = 0 they must vanish entirely. In principle, all that we can say is that as z approaches the limits of the atmosphere (possibly at infinity) the perturbation must go to zero. Actually, as far as practical calculation is concerned, it is likely that the perturbations can be taken as zero outside some finite z interval containing z = 0. If hi and h2 are positive, real constants, we can then state that u = v = w = = 6p = p = 0 for z > h1, z < -h2 (65)

-39d d where the h's will have to be specified for any given situation, and may be infinite. Of course, h1i and h2 2From the vanishing of u and v, we conclude that =-V= o -0 Z>1, i XI (66) Also since u and v vanish, we deduce from Equation (41) that __ zQ and, hence, that ^D)V- o J=O~V tZ' f,,z (67) The vanishing of bp and 5p yields no additional information. The question arises whether these conditions suffice to specify the problem uniquely. When we combine Equations (61) and (62) into a single equation in W (say), we will have an integro-differential equation with fourth order derivatives. We expect to require at least 4 boundary conditions. Indeed, we do have 4 conditions on W and its first derivative. We have also included a pair of conditions on e which can be looked on as integrability conditions, though rather strong ones. We will see that these conditions suffice to determine the critical stability condition unambiguously.

IVo THE ONSET OF CONVECTION A. The Equation of Marginal Stability We have mentioned that the principle of exchange of stabilities is that the case n = 0 separates the stable and unstable regimes. The validity of this notion for the radiative case is demonstrated in Appendix I. If the equations admit non-trivial solutions in the case of n = O we say that the conditions under which this may happen, are the conditions of marginal, or neutral stability. We then seek the condition of marginal stability by setting n = 0 in Equations (61) and (67). Then, in marginal stability,, _ G. {s'} (i)- ) - (7) k Q() (68) and A/ v ) TJ - J\ 9() (T) (69) Limits of integration in Equation (69) and in the sequal are omitted with the understanding that the range of integration is over values of for which perturbations do not vanish, ie.,- L3 do d A combination of Equations (68) and (69) yields the following equation in W alone. J( I, l );,)X(X) d; -N-J, a f W. —- a X -4o r -40

-41where (<) is the average 3 in -i < ^ and -- - X (rl) Equation (70) then governs the marginal stability of a medium of given Po Subject to the boundary conditions, there are only certain values of X for which solutions of Equation (70) can be found. For a given value of a, the minimum value of X that admits solutions gives the critical condition for the onset of convection in the mode described by a. It is our aim to find these minimum allowed values of X. B. A Variational Principle We can rewrite Equation (70) as i.!V)r (' ) -(, P ~)(0 ) (72) where F} = (JD- a) W (73) and A K) i - ( (74) where 6 is the Dirac function. Let W be a solution of the adjoint equation Jw~...... \,A,,! I X i) (75)

42where nox/'^. /99mt-rIFl /) )^c ^-\ QC~ T(76) The boundary conditions on W and F are, from Equations (66) and (67), /77/ NJ4 ThA 3 X ) i (77) and on W and T D - T-o -_ (! ) iKT7-VV7T I4 - ~-I)5 0 d(78) Then, as we shall now show, the expression ~~~X ~~~~~lr~~~~~~ I A Scot (3,(79) - J Ei- (W nroi) (A JC I!,J — is stationary with respect to small, arbitrary variations in W and W which satisfy the boundary conditions. For suppose we know that there are functions W,F,W and F satisfying the equations. Let us make the guesses q and Fi r at F and F, where \f and vanish on the boundaries. By means of Equations (73) and (76) we can then find correspondingW*JandV+ T i and make them satisfy the boundary conditions.

-43Now, if our guesses are reasonably good, s rand p F will be small quantities. In that case, the error we make in computing X by Expression (79) will be, to first order in small quantities, CL I( V() ck JJrJ/al('F()Vij)(Sc (LIF- O8 With Equation (79) this reduces to (81) Since W and F satisfy Equation (72), the foregoing equation becomes - (,,. i,., - v\ v 1,,.',,';(82) The second term in Equation (82) may be transformed by partial integration; we find that

-44In this partial integration, the boundary conditions on W and J\Vwere required. Now, since the Fts must vanish at the integration limits, we have by the definition of adjoint integral operators [cf. Morse and Feshbach (1953), p. 877], r r L/ Ya J ai'a rXXt 4 s-> (84) So that on introducing Equations (83), (84), and (75) into the Expression (82) we find s(l0)- =o0 (85) Hence, the Expression (79) for X is an extremum, and moreover, in order to specify the extremum expression the boundary conditions stated in Equation (77) are both necessary and sufficient. C. The Constant Gradient Approximation The variational Expression (79) can in principle be evaluated for any particular equilibrium situation specified by a B/<(3. In particular, by careful calculation based on comparison with radiativeequilibrium model atmospheres, it should be possible to determine which spectral types are stable against thermal convection. But it would be rather a long process to do carefully for all spectral types. Indeed, one might wonder whether the stabilizing effects of radiative transfer are important in any stars which satisfy the approximations introduced. It would, therefore, be of some value to produce by less elaborate calculation than that called for in general, a survey of possible

i45instability conditions. To do this, we will make some idealizations of the situation. First, let us suppose that the stabilization outside the superadiabatic zone is very effective. Then, the perturbations would die off very quickly for I | > o (i.e., for j 1 )~ This suggests d setting hi = h2 -, or in other words treating the boundary conditions v r - A T -' oT -- ^ G *. i.'o Y ( (86) Then, as a further approximation we will neglect the variation of P across the superadiabatic layer, taking =,' <^ (87) With Conditions (86) and (87), the equation of marginal stability becomes {lOH'onincivA1-Fc - isw(3, (88) )TrcsA)a? I ^ (88) where, again, F =(] - Q)~)N (89) The boundary conditions and the form of p now resemble those in the classical theory. In Chapter II we saw how the form of the radiative heat equation reduced to that of a thermal conduction equation when -'T -^. We would, therefore, anticipate that Equation (88) should reduce to Rayleights differential equation when — 9. That this is the case is shown in Appendix III.

46In this constant-gradient approximation the extrema of X may be computed approximately by the variational principle, which in the present approximation becomes A DF') JJ- tipJ J-(__________ 2 tI F(J)T (90) It may be seen [for example, in the discussion of the eigenvalues of symmetric kernels by Courant and Hilbert, (1953); Chapter III] that the extrema are minima for positive X. D, Calculation of Xarit To approximate the minimum allowed value of X for a given choice of a, we begin by selecting an J satisfying Fct) F(-&=o (91) which is the third of Conditions (86). Also, since t and 8 are proportionalo We may be guided in our choice of P by any intuition we have about the temperature distribution. Having chosen F, we may find W by solving Equation (89) subject to the remaining boundary conditions. In general, the original choice of Fwould involve one or more variational parameters, and we would select those values of the parameters which give minimum values for X in Expression (90).

_47The situation exists in the present problem that the exact eigenvalues are known for the classical theory (i.e., the case t- o). 0L There, the function F(0 = (c^ Y%, m=l 3 (92) produces values of X in good agreement with those resulting from the exact solutions. The f\,is essentially a variational parameter (though not all values are allowed by the boundary conditions) and the value, - 1 minimizes X as a function of lfor all a. With this choice of,i) a value of about 1715 is found for the critical Rayleigh number when the boundaries are rigid; the exact answer j as seen in Table I of Chapter I, is 1708 to four significant figures Hence, in the present calculation, we expect that the trial function El (X) = ^ T? (93) should give reliable results for large e However, for small /', it is possible that altogether different trial functions would be preferable. But the inClu;ion of variational parameters would be somewhat involved here in that minimizing the expression for X with respect to these parameters involves a long numerical process in general We willa therefore, follow the alternate procedure of computing X for two quite different trial functions which involve no variational parameters, and selecting the one of these which gives the minimum values of Xo The justification for this procedure is the a osteriori one that where the procedure is needed ( r/ not large) the different trial functions lead to nearly

-48the same eigenvalues, so that X is fairly insensitive to the choice of trial function. Our first choice of trial function is T of Expression (93); this is a natural choice because we know it to be quite adequate for large /C. On the other hand, we might expect that for transparent layers (small 4I ) the temperature disturbance would go to zero more steeply than in the opaque case. At least we are familiar with this sort of behavior from ordinary radiative transfer problems. Another motivation in our choice of an alternative trial function is the formal equivalent of the preceding reason. Namely, since the derivative of is left unspecified by the boundary conditions, we would like to vary this property of F as much as is possible. A function with vanishing derivatives on the boundaries seems appropriate. Therefore, we are prompted to adopt as a second trial function ok ( j) -+ =Wt i3K (94) The amplitude factors in both trial functions are arbitrary, and have been left unspecified since they have no effect on the results. The two trial functions are exhibited in Figure 3; the normalization of both in the figure is such that jT c) =( I

-4977(N) \Ti~s IC~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ \ 1S —- -------— [Figure 3. Illustrating the Trial Functions Even though we have chosen trial functions containing no variational parameters, the calculations of X by Expression (90) are rather lengthy. The procedures are carried out in Appendix IV for ~ and in Appendix V for We summarize the results here. Using as trial function, we obtain from Expression (90) the result A-x - -^ - t (b *ic k QC (95) ^+'n-H' -+ I/\'b" I^ "s L ^I "~ c' 1T

-50where Ca x2TC (96) and. +1 U^1-)- = -^ — 3: —--- ^ ^-; ~ C^ (97) + I with _. ^ - (98) The integral U(C)j has been computed by a 5-point Radau formula, and the term involving U is significant in the result to about only one part in ten thousand. We have tabulated U in Appendix VI. A twopoint Radau formula results in the approximate expression UVaZ)" - ^^ (99)

-51For A we obtain the formula 4AN - - CCS l 3 - C ( A itt +-11 lf t-+o1 | q t /, a q -I C OLI T(100) tt1-HXJ —Rt^a-f^ lO a a t(az |) where.gal = G+-^ (ioi) and kx Via:~r)_ 7 (102) Y(at) A A, wlx —- -—, ( 02) The two-point Radau approximation gives =-e y Oa1^ = -t (103) The term in X2 which involves V makes an even smaller contribution than the counterpart U-term in X1 and we may neglect it.

-52Given these formulae, the procedure for computing critical values of X can be carried out graphically. For a fixed, we merely plot X versus a, and locate the minimum. The minimum value of X is then Xcritj and we know that a radiating layer of gas is stable against convection unless XI i k\crt. (104) The values of Xcrit and acrit for the two different trial functions are given in Table IV. TABLE IV VALUES OF Xcrit AND acrit T al Xa2X2 -crit lcrit crit 2crit Diff. 0 4.77 70.10 4.8 69.03 1.5.1 4.75 68.09 4.72 67.13 1.4.5 4.55 6o042 4.55 59.82 1.0 1 4.32 51.63 4.4 51.51 0.2 3.56 25.30 3.7 26.42 4.3 5 3.27 14.41 3.5 15.62 8.1 10 3.14 4.80 1000 3.11 5.27xl0-4 3.1 6.52x10-4 13.5 Our first concern is whether the results obtained are at all accurate. The transition of our Equation (88) into the Rayleigh differential equation at large shows (cf. Appendix III) that 3r2X goes into a Rayleigh number. Thus, our present results at T = 1000

-53give radiative-Rayleigh numbers 371 )tcr = 1716 and 3L cL[ -- 1956. The former result compares favorably with the classical 1708. Hence, the difference of 13% in Xlcrit and X2crit are of no concern since we know that Xl is quite accurate. As T becomes smaller, the crit two computations give more similar results. The fact that two quite different trial functions give comparable eigenvalues lends confidence to the accuracy of the results. We will at each T adopt the minimum of the two results. Thus we take Xcrit = Xlcrit for T?1 and Xcrit = X2crit for T$1. The nature of the differences in Xlcrit and X2crit are quite reasonable on physical grounds. Our choice of f was motivated by a feeling that, for transparent layers, the drop to the boundary value is quite rapid. In this way we might expect 1 to be more nearly a solution than 1 for small T. The results bear this out. However, the differences in Xl t and X2cr for small T, are small. This is crit crit because ~ and L do not differ so much in actual value as in derivative. Since at small values of T, the smoothing of temperature fluctuations is primarily through emission losses and depends mainly on the value of O (or F), the two trial functions give similar results. On the other hand, for large T radiative transfer behaves like a diffusion, and diffusion is determined primarily by the derivative of ~ (or F). Hence, we find the relatively large difference between Xlcrit and \2crit at large T plausible. To summarize the calculation, we plot Xcrit and acrit as functions of T, drawing the curve through the preferred values as discussed above for the case of Xcrit. The values for acrit can be

-54considered identical for the purposes of our accuracy. Figure 4 shows Xcrit as a function of T, and acrit versus T is plotted in Figure 5. E. Physical Interpretation As is the case with most of the dimensionless parameters describing hydrodynamic stability, the parameter X can be arrived at as a ratio of certain forces from purely qualitative arguments. We offer here a semi-qualitative derivation of the instability criterion of the last section in order to clarify this remark and some of the physical significance of X as an instability parameter. However, when we speak of this as a physical interpretation it must be made clear that we have succeeded only in making plausible that X is indeed the relevant stability parameter in this problem. The question of why the critical values of X are what they are is in no way answered by these crude arguments. Thus,X, just as the classical Rayleigh number, can be thought of as a ratio of buoyancy force to fiscous force. In a rough way, we can take as the requirement for instability the condition buoyancy force ~ viscous force. For a perturbed fluid element of dimension I, we would have' e', \2 l(105) for instability, where r is the coefficient of viscosity. Suppose that in the lifetime, t, of the disturbance w represents a typical speed and V' is the distance traversed. Then Lv A) v -I Al&~~~~ = I ~ Lt~ t (106)

70 s60 -^e ICRIT \2 CRIT 50 40 50VJ1 30 0.,.!,,,,!! t. z I z,, i I!. I,' 0.02 0.05 0.1 0.2 05 1.0 2.0 5.0 10.0 20.0 50.0 T Figure 4. Xcrit Vs.

5.0 4.0 7_o~~~~~~~~~~~~~~~~~ _\x~~~~~~~~~~~~~ O 0 3.0 I 1 I I 11 1111 __I_ I I 11 I I 1 I __11 I I 1 0.04 0.1 0.2 0.5 0.7 1.0 2.0 5.0 10 20 50 T Figure 5. acrit Vs. T

-57Also, as in Equation (50) g = e-g4e0 (50) Let us approximate G by 0-~ o. Q (ia) Then, with Equations (106), (50), (107) and the definition (lo8) 9 = a 1 the requirement for instability becomes PW~~~~~l\ > I|~~ ~(109) We see that the greater is i, the more unstable the perturbation is, and in particular i d is the most favorable possible case for instability. Then, the condition for any instability becomes -;I (110) V In general, we should expect the lifetime t to be set by the dominant dissipation mechanism. If, for example, this is thermal diffusion w- - t- c Te (111) where X is the coniucticn coefficient. Then the instability condition becomes ~ 1 ~ i (112)

-58Of course, as we have seen, the critical Rayleigh number is not unity, but more like 10. The variational procedure in the classical theory with the trial function f gives sTiGa I+clar 0o. - at | X: +:a - (113) as the criterion for instability to wave number a. This criterion is the particular one for the boundary conditions we have adopted, the so-called rigid boundary conditionso Here, as in Equation (96), q is the three-dimensional, non-dimensional wave number CL - (a,a. )/L (96) For the case of a radiating medium, the lifetime as given by Equation (28) is +- ^ I___________ (114) The instability Condition (110) becomes \LB = X( \ _\- Cort, ~ci _) (115) This, like Condition (112), is very crude. But we know that the radiative instability condition should reduce to the diffusive one in the optically thick case, and X in that case should be replaced by 7/3K2 as in Section D of Chapter IIo With this expectation, and Condition

-59(113), we can make a reasonable guess at the missing factors in Condition (115). Using the simplxt form that will produce the correct limit at large T/a, we obtain as the approximate criterion for instability, J- = rI q [I [ Vf,- + a tc- ~ -(116) The critical values of X in this approximate criterion differ only by about 20% from those of the more precise Condition (95).

V. CONVECTIVE INSTABILITY IN EARLY STARS A. Plan of This Chapter In Section C of Chapter I, we mentioned the need for a theory of convective processes in stars and suggested the possibility of extending current work in non-linear hydrodynamics [e.g., Malkus and Veronis (1958)] to application in stellar problems. The study of the linear problem under stellar conditions was proposed as a necessary first step in such a development. In Chapters III and IV, a restricted linear problem of this sort was considered; it was undertaken to investigate the convective stability of a radiating layer of fluid with (nearly) constant mean properties. We saw that the stability of such a layer is described by two non-dimensional parameters, a radiative analog of the Rayleigh number and the optical thickness of the layer. The natural next step in the development of the convection theory is the extension of the linear study to include the effects of variable density. However, this will not be undertaken here, and we will conclude the present investigation by examining the stability of those stellar atmospheres most suited to the approximation of the present theory —the stars of early spectral type. In this chapter, we will evaluate the stability parameter X for a number of model atmospheres of class B and early A. We will see that one can expect convective instability in the atmospheres of most of these stars. In the case of early B stars, we will make a hypothesis about the relative importance of convective and radiative energy transport, while in the case of the early A stars we will see that -60

-61there is some evidence that the convection is stabilized by rotation in at least large regions of the stellar atmospheres. Finally, we will close the chapter with a summary of the results of this dissertation. B. Expression of X for Application to Stellar Atmospheres To test for the possibility of convection, we need to know a number of properties of the equilibrium state for the evaluation of the instability parameter, X. Many of these are tabulated in the published model atmospheres, but the practice is not uniform, and supplementary calculations are often required. We may minimize such calculations by re-expressing X in terms of quantities most usually tabulated in the model atmosphere work. First, we may use the hydrostatic condition, Equation (45), to write lT _T CQT %T s - -" ~ "e T ^(117) <^Z^ a' The quantity -- oT is usually tabulated in preference to dT ckLat d I in model atmospheres, though even the former is not always presented. We may then introduce the ideal gas equation = k.T ((118) wherek is the Boltzmann constant, ~ is the mean atomic weight, andT.

-62is the mass of unit atomic weight. We find that 4,__m con wAT (119) ars " k Age so that A (S)T T S ^_ /G A (120) where A = /(/>a T (w) Also, in an ideal gas, the expansion coefficient is OC = 7, (122) and if n is the coefficient of viscosity,./v= \ (123) Hence, Equation (71) for X becomes - (\Ib ( A) (K Cp) (124) i-~ [^) \^} ^A K

-63This expression for X is then to be evaluated and compared to the appropriate value of Xcrit for the model. The quantities on the right hand side of the expression have been taken as constant in the development and the choice of average values must somehow be made, though this should not be an important consideration in cases where the approximations of the theory apply. We may also note that most of the quantities on the right hand side of Formula (124) are tabulated as part of the results of the model atmosphere calculations, except for principally c and rb. Rosa P and UnsBld (1948) have tabulated... as a function of T and for stellar material with 15% helium by numbers of atoms. Edmonds (1957) has similarly tabulated r for the 15% helium abundance. Let us then turn to an examination of the convective stability of stellar atmospheres. C. Convective Instability in the B Stars Hydrogen and helium are the most abundant elements in stellar atmospheres and, wherever these elements are in an intermediate stage of ionization, the temperature gradients may become superadiabatic [Uns6ld (1930)]. Depending on the effective temperature, Te, the surface gravity, g, and the chemical composition, a star may have ionization zones of H, He I or He II near its surface. A summary of the various possibilities for different Te and g and for a 15% numerical abundance of helium is given by Mme. Pecker (1953) in her extensive study of stellar convection zones. Her considerations show that the possibility of convection, even due to He II ionization, does

-64not arise if the effective temperature is above 35,000~K, and the surface gravity is below 10' cgsu, since the ionization of these elements is essentially complete beyond this limit. The limit seems to exclude almost no known stars, and the possibility of convection arises in most stellar atmospheres. 1. The BO Stars (Te = 36,800~)* In the BO stars, temperatures are so high that hydrogen is almost completely ionized all the way to the surface of the star. Even the first ionization of helium is reasonably complete, and the only possibility for an ionization instability lies in the second ionization of helium. Underhill (1950) and Rudkjbing (1947) have computed models of such stars and found that convective instability is indicated by the Schwarzschild criterion. Traving (1955) finds a similar result in his discussion of the BO star, T Scorpii, but he concludes that convection cannot develop because of the stabilization by radiative transfer. Traving's argument concerns properties of the fully developed motion, and we will summarize it briefly in the next section. For the purpose of discussing the stability of the BO stars, Underhill's tabulations are the most convenient of those mentioned. The parameters in this model are Te = 36,8000 log g = 4.2 ________H/He = 85/15 * In the assignment of spectral types to various models, we will always adopt the type mentioned by the computer of the model. This may mean deviations of a few tenths in the type assigned to a given effective temperature. However, this has no bearing on our considerations, and we use the spectral designations only as nomenclature. The meaningful quantity here is the effective temperature.

-65where H/He denotes the ratio of numbers of hydrogen and helium atoms. Figure 6 shows the adiabatic excess, A, as a function of the mean optical depth, T, in Underhill's model, where ft - -iKCd (125) and K, it will be recalled, is the absorption coefficient per cm. 0-,2 T A 0 I / * I I \ -- -.02 -.04 - 1 2 3 4 Figure 6. A Vs. T in Underhill's Model The points plotted are those obtained directly from Underhill's paper; the solid curve is used to extrapolate these computations. The dashed curve is the form of A called for in the constant-P approximation. We see that the value of A drops fairly rapidly just outside the superadiabatic layer. This means that the gases just outside the layer are quite stably stratified, as is required by our boundary conditions for the constant-n approximation. Also, the temperature

-66and density respectively vary by 7% and 30% across the superadiabatic layer, and this seems to fulfill the Boussinesq approximation, i.e. that the mean properties of the fluid do not change appreciably across the superadiabatic layer. The optical thickness of the layer is T = 1, and using a mean value of K across the layer we find d = T/ = 1.2 x 108 cm. With these data and the tables of Rosa and Unsbld, and of Edmonds, we find X = 500. The critical value of X for r = 1 is Xcrit = 50 (cf. Table IV). Thus, X exceeds the critical value by an order of magnitude and there seems to be ample reason to expect thermal instability in Underhill's model. Of course, the question raised by Traving as to the later development of the motion cannot be answered by stability considerations alone, and we will defer discussion of this matter to the next section. 2. B2 Stars (Te= 20,400~) As we decrease the effective temperature from BO, the zone of He II ionization will move deeper into the star. Saito (1954) has computed models which show this effect. Though these models do not closely satisfy the radiative equilibrium condition (i.e., constancy of flux conditions), they probably provide an adequate description of the temperature gradients. We find that for Saito's hottest models, i.e., for Te = 20,400~K log g = 3.8, 4.2 H/He = 85/15, the gradient does not become superadiabatic until a mean optical depth of about 8 is reached. Saito's calculations do not extend beyond optical depth 10 and, hence, do not provide adequate information for

-67examining the stability. But even if convection does occur, it will probably be too deep to produce easily observable effects. 3. B5 Stars (TR = 15,500~) If we turn next to Saito's models withTe = 15,500~ and the same chemical composition and gravities as above, we find that He I now produces a superadiabatic gradient and no effect of the He II ionization is detected above mean optical depth of 10. Figure 7 shows the variation of A with physical depth, A, in the atmosphere of the log g = 3.8 model. The origin of 7 is set arbitrarily at mean optical depth 1. Thus, if we use T to denote mean optical depth, {? - d c.) ( ) (126) as may be seen on integrating Equation (125). A summary of values for the stability is shown in Table V. TABLE V STABILITY CRITERION QUANTITIES FOR SAITO'S MODEL WITH = 15,500~K g d p r X Xcrit 6.31 x 103 3.6 x 108 1.6 x 10-9 1.5 7.6 x 105 43 1.59 x 104 1.4 x 108 3.0 x 10-9 1.8 1.5 x 105 39 We see that instability condition is easily met in B5 stars. However, the Boussinesq approximation for density is not well satisfied; a factor of about two in density occurs across the superadiabatic layer.

-68The temperature, on the other hand varies by only about 20%. Also, A drops off quickly outside the superadiabatic layer, as in the case of Underhill's BO model. It seems reasonable to accept the result that there is marked instability in B5 stars. 0.3 0.5 0.8 1.2 1.8 2.4 3.0 7.064,04.02. I I I I I I - z x 10-8 (cm) 3 2 1 0 -2 -3 Figure 7. AVs. z for Saito's Model with log g = 3.8 and i = 15,500~ The situation with B stars thus is that convection due to He I ionization disappears as we go from B5 to B2. However, at the disappearance of He I convection, the Te is not yet high enough for He II ionization to appear near the surface of the star. It is not until we reach temperatures over 30,000~ that He II produces effects near the surface. Hence, there would seem to be a small gap of spectral

-69types that do not show the effects of convection. As we have mentioned, it is not clear whether convective effects become detectable even in BO stars. The remarks of the next section are directed at this problem. D. The Effectiveness of Convection in the Early B Stars We mentioned in the preceding section a suggestion advanced by Traving (1955) that the development of convection in the BO stars is inhibited by radiative transfer. Traving argues in terms of a mixing length model for the developed motion. He equates the mixing length to the scale height and computes the ratio of energy excess in a turbulent element to total energy radiated in the lifetime of the element. Vitense (1953) has given an approximate formula for this ratio. At a mean optical depth of 1 in his own BO star model, Traving -4 finds this ratio to be 10. Traving, therefore, concludes that convection will not be important in BO stars. Our own considerations have shown that in Underhill's BO model there is convective instability with - = 10. In this section, we will attempt to find an approximate criterion for the effectiveness 1\ of convection as indicated by the degree of instability, It should be emphasized, however, that we will utilize the results of experiments whose application andf interpretation in the present context is not unambiguous and our conclusions must be viewed with caution. In his discussion of convection in stars, Unsold (1955) cites the experimental work of Schmidt and Saunders (1938) on Rayleigh-type convection. Schmidt and Saunders observed a transition from laminar to turbulent convection at a Rayleigh number of 47,000. [The critical

.70number for the onset of convection was found to be 1700 by Schmidt and Milverton (1935)]o The transition to turbulent motion found by Schmidt and Saunders is accompanied by a discrete transition in slope of the curve of heat transport versus Rayleigh number. Unsold remarks that, since the Rayleigh number in stars usually exceeds the value of 47,000, we may regard all stellar convection as turbulent. While we will not take exception to the viewpoint that stellar convection is normally turbulent, we may express concern at the application of the laboratory results to stars in general. The experiments of Schmidt and Saunders, and of most investigators, provide solid boundaries for the convecting layer. Moreover, the conditions of the Boussinesq approximation are well met in the laboratory. We will see in the next section that these restrictions are already violated in early A stars; they are even less closely approximated in later spectral types. On the other hand, we have seen that the case is somewhat different for the early B stars. There the Boussinesq approximation is not a bad one, and the gas surrounding the superadiabatic layer is stably stratified so that penetration beyond the unstable layer is likely to be small. Hence, any application of laboratory results to stars should be restricted to these early B stars. Another major difference between the laboratory situation and the stellar conditions is that radiative transfer acts in the latter case. That, of course, has been the concern of this work and we have learned that in the constant-gradient approximation the classical and the radiative theories are formally equivalent with a differential operator replaced by an integral operatoro Also, we saw that X

-71is physically analogous to the Rayleigh number, the difference being in the critical value of each for the onset of convection. We will then make the hypothesis that we may apply experimental results to early B stars by replacing R/Rc by X/Xc in these results, where RC and XAc represent the critical values of the Rayleigh number and of X for the onset of convection. The experimental results we will use are those of Malkus (1954b) who has studied the heat transport in layers of distilled water and acetone up to Rayleigh numbers of 1010. Malkus has found a number of discrete transitions in the heat transport as a function of Rayleigh number, some (including apparently the one noted by Schmidt and Saunders) are merely transitions from one type of turbulent flow to another, Malkus points out that it is possible to represent the total heat flux~j) across the convecting layer by wH =3 p ^ f) (127) where'iand a are constants, K is the conductivity, and f is the (constant) static gradient. The experimental values of a and kodiffer on different sides of one of the aforementioned discrete transitions, but this one occurs at a high enough K to be of no interest here. Also, there is some question about the agreement with other experiments, but for the cases of interest to us, the agreement with other work is good. Let us now assume that even in the convective case, the heat flux across the layer due to thermal conduction is KS. Then, the

-72average fraction of convective energy transport is r -C (-(128),' ~,," This neglects the distortion of f by the motion. Then, let us introduce the quantity x such that R = x K RcC (129) so that - othi m er w/(130) By the hypothesis made earlier, we may replace / by X/Xc in these expressions. Then, in a radiating layer fl = (131) where A and a are constants. One objection to this procedure which may be raised is that in the experiments convective transport competes with heat conduction, but in stars convective and radiative transport compete. Since radiative transfer is more efficient than heat conduction, we expect that the estimate of r made in Equation (131) may be too large. Another difficulty to be noted is that P will not generally be constant in the equilibrium configuration of a star and the expression for | is, at best, a gross estimate for that reason as well. To adequate accuracy for our purposes, Malkus' determination of the constants in the heat transport law leads to A = 1 and

-73" a = 0.3. Thus 0,3 r -'1-(') (132) In Underhill's model with X,/X =.1, we find r =.5. Hence, our hypothesis leads us to suspect that an appreciable fraction of the heat flux across the unstable layer is carried by convection. However, since the unstable layer is not extensive in optical depth, the consequences of convection for the structure of the star may still be small. In order to see what the effect on the structure of the star may be, let us adopt a model in which a fraction,, of the total flux is transported across the superadiabatic layer by convection. Elsewhere the flux is supposed to be entirely radiative. If _ is the net radiative flux, T1 is the optical depth at the top of the superadiabatic layer, and T is again the optical thickness of the layer, we have Of-sTe ~ For F(r, = (, —T. Cot t,<T<t (133) rTl Far t t,+t~ This approximation for Y is meant only for early B stars in which the penetration of motions into the subadiabatic layers is expected to be small. Also, the question of how the mean absorption coefficient is to be defined is left open. Presumably, the choice of mean will be that which gives the best approximation to the radiative equilibrium condition.

-74For a grey atmosphere in local thermodynamic equilibrium, but not necessarily in radiative equilibrium, the Eddington approximation leads to the expression [cf. Unsold (1955), p. 142] <rTF t(O)+ + TF (134) With the assumed form of, we then find T ()= [ Te [t + 7- (P )] (135):where., %p v, ) - 3r(t-), r<' ( (136) For a model in radiative equilibrium, we would have ) (7) O. Considering Underhill's model once more, we have T 2 T = 1, and from our semi-empirical discussion of her model, r = 1/2. Then C($r= q ^(-X) < 3 (137) ( t 3

-75and we see that the discussion of convection does not lead to drastic revisions of the temperature distributions. At most, a 3% change in is required at r = 3. The conclusion drawn by Traving, that convection will not develop, is borne out by these considerations to the extent that no great changes in structure are produced. The reason seems to be more that the unstable layer is optically thin rather than that convective transport does not occur. The questionable nature of the calculations of this section makes it desirable to look for an observational check on their usefulness. We should like to compute spectral changes resulting from the changes in temperature distribution. The hydrogen discontinuities would probably be suitable spectral features for this purpose. Unfortunately, it seems that the uncertainties in the abundance of helium in B stars [Jugaku (1958)] would mask any effects of the magnitude we would be interested in. Hence, any further progress in this problem must come from improvement in the hydrodynamical theory. It is hopeful that some definite progress has been made in understanding the discrete transitions in convection [cf. Malkus and Veronis (1958)].

-76E. Instability in Early A Stars In studying the stability of B star atmospheres, we began at le= 36,800~ and moved to successively cooler stars. We noticed that in B5 stars, as contrasted to the BO stars, deviations from the Boussinesq approximation appeared. It should come as no surprise then, that deviations from the conditions of this approximation are even larger in the A stars. In spite of this, it seems worthwhile to carry the stability discussion to early A stars where we encounter instability due to hydrogen ionization. In this way, though the conclusions drawn must remain qualitative, we will at least gain some idea of the restrictions to be anticipated in later work directed at studying hydrogen convection zones. 1. AO Stars (Te= 10,700~) For the purpose of examining the instability of AO stars, two of the models calculated by Saito (1954) are convenient. These are for Te. = 10,700, H/He = 85/15, and log g = 3.8 and 4.2. Saito tabulates the adiabatic and radiative gradients and the geometrical depth as functions of mean optical depth. A similar tabulation has been given by Ueno and Matsushima (1950) for log g = 4,09 and the same effective temperature, ie = 10,700~; their model differs from Saito's in excluding helium and in having a slightly different temperature distribution. Figures 8 and 9 show A as a function of T and X respectively as given in the three models. As previously in this chapter, i = 0 at T =1. (In the Ueno-Matsushima model, the tabulated X is too large by a factor of 10 for T < 1. This has been corrected in Figure 9.)

0.18 A: UENO a MATSUSHIMA, LOG g 4.09 Q14 ~ \ A B: SAITO, LOG 942 0.12 CL: SAITO, LOG 9g3.8 008 0.06 -C 0.04 - \\l f \ 0.04 - 0.02 -0.02 -0.04 0 1 2 3 4 5 6 7 8 9 10 11 12 13 T Figure 8. A Vs. 7 for T — = 10,700~

g —-- 0.16 A: UENO a MATSUSHIMA, LOG g =4.09 | \B sA, ~o B: SAITO,LOG g =4.2 12 - C: SAITO, LOG = 3.8 0.08 ~~~~~~~~~0~~~~~~~~~~~~~~~~~~~~~~~~~~ 0.04 2 01 0 -- -2 -3 -4 -5 -6 -7 -8 Z xlO (CM) Figure 9. AVs. z forTe 10,700~

-79In each of the Saito models, there are two potential instabilities. One is due to the ionization of hydrogen, and the other deeper one is due to the first ionization of helium. Comparison of the Saito models with the Ueno-Matsushima model indicates the diminution of the hydrogen ionization zone by the presence of helium. It would seem that the introduction of the helium ionization zone would more than compensate for this weakening of the instability, and it is plausible that two zones would merge once convection developso This may not be the case in cooler stars where helium would curtail the hydrogen convection without producing an unstable layer near the surface of the star. For the various superadiabatic layers in the figures, the Boussinesq approximation is poor; the zones are each two or more scale heights in geometrical thickness. For general information we may, nevertheless, mention that X/Xc is a factor 2 larger for Saito's log g 4.2 model than for his log g = 3.8 model with the great density of the former outweighing the effect of large d in the latter. We also notice a tendency in these models for layers just beneath the deepest superadiabatic layers to become rapidly subadiabatic, hence, stable. 2. A3 Stars ( Te = 8900~) As a last example, let us consider the models computed by Owawa (1956) for Te = 8900~, log g = 3.5, 4.0, 4.5 and no helium present. Osawa has iterated his model calculation to achieve the condition of radiative equilibrium quite closely. This has the advantage, for us, of providing a rough check on the adequacy of the non-iterated models at J^ = 10,700. Osawa has not tabulated the gradients nor z, and these must be found by supplementary calculation. The adiabatic and radiative

-80gradients may be computed with formulae given, for example, in Unsold's (1955) book and z may be found by numerical integration from Equation (126). The results of such calculations are shown in Figures 10 and 11 giving A respectively as function of T and of Z. Only hydrogen ionization zones appear, though we would surmise that the presence of helium would produce other, deeper superadiabatic zones. The optical and geometrical thicknesses of the hydrogen ionization zones for Osawa s three models and the Ueno-Matslushima model are shown in Table VI. For the Osawa models with log g = 4.0 and 4.5, we have had to extrapolate the plots of A versus T and z beyond the point to which the model has been computed. TABLE VI d AND T FOR THE HELIUM-FREE MODELS Model Te(~K) log g d(km) T Osawa #1 8,900 3.5 4.5 x 103 8.6 Osawa #2 8,900 4.0 1.5 x 103 12.1 Osawa #3 8,900 4o5 4.0 x 102 15.5 Ueno & Matsushima 10,700 4.09 2.4 x 103 2.4 We see that the geometrical thickness of the hydrogen ionization zones does not change appreciably from X = 10,700~ to 8,900~, but the optical thickness does. Some of the difference in optical thicknesses is due to the differences in absorption coefficient tables used by the different authors and also to the inclusion of 4I in the absorption coefficient by Osawa. But some of the increase must be real,and probably

2.6 2.4 22 2.0 1 LOG g 3.5 1.8 L* \ \ 2 LOG g =4.0 3 LOG g =4.5 1.6 1.4 1.2_ 1.0_ 0.8 0.6 3 / 0.4 0.2 0 -0.2 1 2 3 4 5 6 7 8 9 10 11II 12 13 14 15 16 17 T Figure 10. A Vs. T for Osawa's Models. (Te = 8,900~)

2.6 2.4 2.2 2.0 1.8 4.5 * 1.6'LOG g- 4.0 x r'*6 - l \ 3.5 e 1.4 1.2 1.0 0.8 0.6 _ 0.4 0.2 0 -0.2'2 +1 0 -I -2 -3 -4 -5 Z x I0- (CM) Figure 11. AVs. z for Osawa's Models. (Te= 8,900~) eT 8,900O~)

-83arises from the increased effectiveness of HE absorption at the lower temperature. The implication of this increase of T with decreasing would seem to be an increase of the observational (i.e., optical) importance of the convective zones as well as a diminishing importance of radiative effects on the large scale motions. In this connection, we may mention that for Osawa's models #1 and #2 X/Xcrit = 5 x 103 and for his model #3 X/Xcrit = 10, though the precise values depend on the method of averaging such variable quantities as the density. An interesting feature of Osawa's models is the slow rate of decrease of A with z, just below the superadiabatic layers. It seems that this would permit penetration of the motions into the subadiabatic regions. We might, therefore, expect the convection zones in stars of this sort to be of great vertical extent, and convection to play an important role in the determination of the stars' structures. The opposite point of view has been taken by Hunger (1955) who has computed model atmospheres in the temperature range Te = 8,660~ to 9,500, and assumed that convection effects are unimportant. Hunger's conclusion is based on the same line of reasoning which we described in Section C of this chapter in connection with Traving's (1955) work. That is, Hunger computes the ratio of energy excess in a turbulent element to the total energy radiated in the lifetime of the element at T = 1 in one of his models. He finds a value of 10-3 and concludes that no appreciable amount of energy will be transported by convection in the atmosphere under consideration. However, the dependence of A on z displayed in Osawa's models would have led us to suspect that if convective heat transport is important in these stars, it becomes so only at optical depths greater than 2 or 3. Hence, it is not unexpected

-84that the ratio computed by Hunger, is small at optical depths where the material is quite transparent. On the other hand, Osawa has found reasonable agreement between his computed Balmer discontinuities, and those observed by Chalonge and Divan (1952). More recently, Bless (1958) has had remarkable success in matching his observed continuous spectra of A stars with the spectra predicted with Osawa. We must then conclude that either convection is not extensive in A star atmospheres or that if it is, the effects on the structure of the stars is too deep optically to be detected. It would be difficult,, at present, to distinguish between these two possibilities. The possible effects of convection on the observed continua have been examined by various authors who compute two models for each star considered; one in convective equilibrium, and one in radiative equilibrium. For example, Hack (1956) has done this for A stars and has found considerable difference between convective and radiative models. However, Hack assumes that when convective equilibrium occurs, it does so as soon as the radiative gradient exceeds the adiabatic, The result of such an assumption is to include convective effects too high in the atmosphere, and thus, to overestimate the effect of convection. The possibility that convection is stabilized by other factors must also be examined. The rotational velocity of A stars is often quite high, and a typical equatorial velocity is 100 km/sec, [Struve (1950)]. Also, magnetic fields are quite common among A stars, [Babcock (1958)]. It seems likely that coriolis and electromagnetic forces play an important part in convection in A stars.

-853. The Effect of Rotation In general, rotation seems to stabilize convection in a fluid, but Cowling (1951) has shown how non-uniform rotation may lead to shear instability. If we neglect this latter possibility, we can in general expect that A star rotation has an inhibiting effect on the atmospheric convection. The influence of coriolis forces on the convective instability of a horizontal fluid layer with an adverse temperature gradient has been studied by Chandrasekhar (1951). The framework of his investigation is that of the classical Rayleigh-type convection, and it includes the Boussinesq approximation.;Chandrasekhar shows that the critical Rayleigh number for the onset of convection in a rotating layer is dependent on the non-dimensional parameter |- L4 - C where ~ is the angular velocity about a direction making an angle with the vertical. The ratio 0~, where Hi is the critical Rayleigh number in the absence of rotation, has been tabulated as a function of by Chandrasekhar. Because of the formal equivalence between the classical and the radiative Rayleigh theories, in at least the constant-s approximation, we can expect that the values of \c/X~ are the same as fo for givenT, where XO is the critical value of X when = o. For an A star, we may take -- - 10-6 sec-1, and for Osawa's models, d = 108cm and 105 cm2/lc. Then T- 4 x 1010 cos2, At the pole of the star, 0 = 0= and we find from the asymptotic form given by Chandrasekhar for large T, that?c/Xc = 5 x 104. For sawa's models we found X/X~ - 5 x l03. Hence, at the poles of a typical A star,

-86the rotation would suffice to stabilize ordinary convection. However, as Chandrakhar shows, instability under most astrophysical situations is likely to arise as oscillations of increasing amplitude (overstability). Since the situations at the poles of the stars we are discussing are just stable to disturbances of the steady type, we may anticipate that overstability will arise. At somewhat lower latitudes the situation is more complicated, but in the equatorial regions the possibility of steady convective instability is the likely one. At this point, we must stress the qualitative nature of this discussion. All the stability criteria employed are derived with the use of the Boussinesq approximation, and we have already seen that the application of such criteria to A stars is of questionable validity. Moreover, we are actually dealing with the stability of a spherical shell by treating various latitude regions independently. These are but a few of the limitations which must be emphasized here, the qualitative features of the arguments may be correct and we are confronted with some interesting conclusions. Instability seems to arise as overstability at the polar regions of these rotating stars, and as steady convective motion at the equatorial regions. This would lead to latitude variations other than those due to flattening of the star by rotation, and we would expect a latitude variation in certain spectral features. The Balmer jump, for example, should be weaker in stars seen equatorially than in stars seen pole on. [Of course, as Bless (1958) has discussed, other aspects of rotation might lead to latitude variations of spectral features.] We may also wonder about the evolutionary significance of these phenomena. As a star evolves from the main sequence, it expands

887and its angular velocity decreases. Then T and Xc/XO must also decrease, and convective instability spreads to all latitudes. At the same time, the expansion of the star leads to changes in the general instability situation by the increase of d, and the decrease of g and the density in the outer layers. We are confronted with a number of interacting physical processes and perhaps the origin of stellar pulsation is among them. F. Summary of the Investigation The basic motivation in this work is a need for an understanding of the convective processes in stars. This problem involves the solution of non-linear equations for which not even the corresponding linear problem has been solved in general. The specific aim of this investigation is then, to study the linearized convection problem, i.e., the problem of the onset of convection, and to include the effects of radiative transfer as exactly as possible. The spirit of the work has been to make whatever approximations are necessary to facilitate the inclusion of the radiative effects. Thus, the Boussinesq (1903) approximation introduced into this type of problem by Rayleigh (1916) has been retained; that is, throughout the fluid layer whose stability is studied it is assumed that all the equilibrium quantities are nearly uniform. The present investigation begins with a discussion of the heat equation appropriate to a radiating medium. This medium is assumed to be grey and in local thermodynamic equilibrium at all times. For the purpose of exemplifying the nature of the radiative terms, as opposed to the more familiar conductive term, the decay of temperature fluctuations in the absence of fluid motions is studiedo The fluctuations are

-88taken as departures of small amplitude from the temperature distribution in radiative equilibrium. In the case of an infinite, homogeneous medium, the equation simplifies considerably, and the problem may be reduced to quadratures by Fourier transform methods. In particular, an expression for the lifetime of any Fourier harmonic of the temperature perturbation is obtained. It is found that in the solar photosphere the observed *granule1 lifetime of four minutes corresponds to a half wavelength of about 400 km. In the next phase of the work, the radiative heat equation is introduced into the study of the convective instability of a horizontal, radiating layer of fluid. The development generally parallels the classical work of Rayleigh (1916) and its extensions by Jeffreys (1926, 1928, and 1930) and by Pellew and Southwell, and includes the approximations just discussed, An equation of marginal stability is derived, and it is seen that convection may arise when a non-dimensional analog of the Rayleigh number exceeds a certain critical value. It is indicated that a variational principle exists for the calculation of this critical number, and that the effects of variable temperature gradient on the critical value may be admitted, as in the variational principle derived by Malkus (1954 a) in another connection. The critical value of the instability parameter is seen to depend on the optical thickness of the layer as well as the nature of the boundaries. In the case where the temperature gradient is treated as constant and the boundaries are treated as rigid walls, critical values of the instability parameter and corresponding values of the wave number of the most unstable periodic disturbance are computed. The results are shown for several optical thicknesses of the layerO This portion of the work concludes with a

-89, phenomenological derivation of an approximate form for the critical value of the instability parametero The final part of the investigation deals with the convective instability in the atmospheres of stars of early spectral type. The ionization zones in these atmospheres, in which the equilibrium temperature gradients (radiative gradients) exceed the adiabatic gradients are tested for instability in terms of the criterion derived in the main body of the investigation. A number of model atmospheres are examined, and it is noted that the Boussinesq approximation is well satisfied for the superadiabatic layer in the early B stars, but that it becomes less valid as one goes toward later spectral types. In general, application of the stability criterion shows that convective instability is to be expected in the atmospheres of most stars, though there are some early B stars in which the instability occurs quite deep, optically.'The next step in this line of inquiry should be the extension to the case where the layer may have large density variations from top to bottom. The writer plans to undertake this extension of the work.

APPENDIX I The Principle of the Exchange of Stabilities We have Equations (61),MV(, )W = 0(i - T~D (61) and Equation (62) ~-W NJ= v'C [(l iH T)~ (n') X -e( yi3 (62) where 3 G,,, such that W = T J = > o 4 ~ T,.L In general, n may be complex, and 8 and W will be too. The complex conjugate of Equation (61) is i d (D <~) W^ - i (g- o')v^ - ^ ^ ~V o c (AI.1) The quantities 6' and W* satisfy the same boundary conditions as 8 and W. On multiplying Equations (62) and (AI.1) together, we obtain IweCiYtjtui(lYI-) er (a) i cryi- o-'c3-a m-e^ ^^Wd-; - fidL (d ) m CV ^-&) (AI 2) If we integrate Equation (AIo2) over J we obtain, on carrying out the various partial integrations, g1., + dIz+ c%'aC(C3^-( 4- o l-I2) 0 (AI.3) -90

-91where I,. SJ DT J- [(W )-w1 ir (AI.4) I J [('DW ) (DW) + V V Jj1 cA4 (AI1 5) 13 - 0* 9 cK5 (AI.6) I J 8 t) dljE(i- J]0 (J)d 3 (AI. 7) We see at once that I1, I2, and 13 are real and positive. Also 14 is real, since <s_ j( G 0dJ )) g(! r = (AI.8) We can also see that 13>I4. Let 9 (as) = J~(-o IO:H 1) *C "(AI.9) and define the operator by 40e i., eii~ o'~ c~ X ( dt~3 ~(AI.0lo) Now, since ( is a convolution of ~ andK, s-W}~ _~ y[6]S~ <~{ (AI.11)

-92The definition of I is Expression (64) /.)e re i e a+ tr(64) I j@ ) d +- 4 L+7 ho (64) Then ^ rt~ r ^^O }c =t6 e% +L e (AI.12) - o On comparison with Expression (24) for we see that f[K] =g ^ c= t } __R (AI.13) and we may conclude that ^it{ 1C} | < 1 (AI.14) Hence:I~c~i f(g~~ f j t } ~~(AI.15) and by Parseval's theorem, Fr) om | th| <(AI.16) From this we see that -T ~|^ A^Il~l (AI.1)

-93whence JL^ ~3 ~L >fd } * (AI.18) and T > T J(AI.19) If we now let Is = L-LT (AI.20) Eauation (AI.3) becomes t$I,+^^A^+^d A ~ + + 0=dalyIc = o (AI.21) where all the I's are positive. We let n = nr + ini (AI.22) where nr and ni are real. Then, on equating the real and imaginary parts of Equation (AI.21) to zero, we obtain the conditions lJic gCw yI5 +(t&cI+ ^Cd IC) o (AI.23 (tf IL- y2cl I3) Ot o (AI.24) All the quantities in these equations are positive with the possible exception of P. If P is positive, then we see from Equation (AI.23) that nr is negative. If B is negative, then nr may be positive, but ni

-9must be zero. Thus, whenever n is complex, its real part is negative, and overstability is ruled out.

APPENDIX II The Kernel E The kernel in the equation of marginal stability is defined in Equation (64) as E(~&l ) J(@o~ Tt%( } (64) Consider first the integral iL;J p^ + e a, (AII.1) We have dL _ e X, e 1 (AII.2) But this last integral is a well-known Fourier transform [see Magnus and Oberhettinger, (1954); p. 118], and we have 1a I.L, (MIJL+ CA X) (AII.3) where K is the special Hankel function of zeroth order. We may then integrate Equation (AII.3); we find tL(P) L ( o) - | ~jK4() &to (AII.4) -95

-96From Equation (AII.1) we find that' = c(~2i. ~ (i_,2) $(AII.5) whence [Magnus and Oberhettinger, (1954); p. 116) Le(o) - e (AII.6) Then _-1 (T)-^ h7'T- e^^'JKJ-F3 t\&(AII.7) and i ( l& ( 4R |y' )4 T <V (AII.8) The first integral in Equation (AII.8) we have already encountered in Equation (AII.1); thus, we have r(.X\\ - }(<, i t ( I I, (J Q' ) ol; d u (AII. 9),t^; J^ J^ Q,s),-': a

-97We can perform the Tq integration in the foregoing equation since we see that it is the inverse of the Fourier transform encountered in Equation (AII.2). Then we can write (1 Ko O - L 0^ — ^ ~ (, )~, t(AII.10) In the second integral, let -S[ = -^( t? /(All.11) Then tz2-' ~~s = a~ " gS~ = ktT (AII.12) Hence, ^ ('aTl|a =.. T- (AII.13) F'inally, we may introduce an integral representation forJX, namely, 12, /-CS 3l ^(AII.14) [see Magnus and Oberhettinger, (1954), p. 271. With this expression, we conclude that I e~ 7j c ^(iTO^^J~~isr~

APPENDIX III Reduction to the Classical Theory When the perturbation dies out at the condition of marginal stability is Equation (88), (hl~~l)?^)^ - F~~~~n -- - X^V~~~n (88) 4-iwhere F Q= CD&)w t,7(89) and, as we saw in Appendix II, -o e-IT-r'l.3 ^ -1) ji-. ^ s ~ (AIIl15) 4^ Let LL "~~ =^ +0^-~ ~(AIII.l) Then IL( -T'i) I (AIII.2) Also, let us write$ as the operator for theAiL-transform, so that A.-8(AI.3) ~98

-99" Since J-"J in F is always less than unity, we may write Q= ~ ((T) o 4) ~o ( 5 F~ ( J) f'-C(AIII o4) where we assume that all derivatives exist We have then that OiV;t R = a! T ) (AIII 5) where M uTQjI ( J ( ~-j ) (Ir) d/ J )(AIII 6) We can write (-QJ (1-i) )KJ 3/) L/ (AIII.7) or fJQ(d) = z XX XQ F(X)d cx (AIII.8) v 0 If we define A then (AIII.9) then NQ (%) = Ai - + (-I)aQ( ) (AII.I10)

-100Now, using Expression (AIII.2), we may write AQas AQ( C) - JI q X( e d X (AIII 11) r \f -2~ j -J(X By repeated partial integration, we find that gb e! - [I ~-o (AIII.12) whence 0 (lB R!r~l -J AE5 [ b (AIII.13) Let us now investigate the transition to the classical limit as T-* I. In that case, the AQ become AQY (^ b -) l -t (AIII.14) where we have neglected only terms of negative exponential order. For N we have N,^ ( -A) PR Jl~a"'^ _;a~= for y even (AIII.15) and N(0) = forQ odd.

-101By integration by parts, we can show that S* _! r t+ (a- Q- i) ^ ^ - (AIII.16) We further find that 4^0 = a it (AIII.17) No s c C5C ^ " ^W (AIII.18) (AIII.18) [lJ t -+ — ) - + C NL (AIII.19) Also, to 0 ( X (AIIII 20) ^ ^^ = 0^"i: (AIII.21) (AIII.22) Hence, {rf_3o3rf } ( e ()+.(AII.23)

-102and Equation (88) reduces to (re-pac -(y 3i aW which is Rayleights differential equation with the Rayleigh number replaced by 3T2X.

APPEKNDIX IV The Calculation of Xl We wish to evaluate _ cVS ffcki -'K(lSJ31)TI7~' O~ ) (90) where Tl,(l = ^^ Y, (93) (J - jV\T, (AIV.l) and, as we saw in Appendix II, eO ^ds;L 47-7T (AnI.l5) Equation (AIV.l) must be solved for W1, subject to VTJ, 5V, L & S - Iz (AIV.2) The general solution of (AIV.l) is "W, = cJ^O i 1)^Sal + C-+ +Ma + pS - i ~ -t -I (ZT^ (AIv.3) where A, B, C, D are arbitrary constants and q2= a2.2 The equation of marginal stability and the boundary conditions are even: -103

-104let us take even solutions. Then -vJ; = 0", 0L Ja - q (AIV)4) Also, T)'W Q=l (cxAfd^ c^Aa^ - aCb ^^A^T - ^ (AIV.5) Applying the boundary conditions we find A eh ^ ^^.Ao. h(AIv.6) and ^.,_ q (AIv.7) Hence W1 is specified for any a. We may now carry out the integrations. First, /t~~~~~~ "F'^S'= t~ fcwL~~io(AIv.8) -'-t j i The more difficult part is J.= dr Ti J/ IFWw 17(K)^I'ea5^ (AIV.9) or le~~~~~~~~ (Aveo r IC~~~~~~~~~~~~~~~~ Joos9d2 Jct csS0~YQ&T i L~~~~~~~i

-105Consider first the integral J, = cs 1 e dAs (AIV.11) i We can write Jo, t J-, r"- - ),^ + ia As 2 (AIV.12) and integrate. We obtain J, = act^t [-o-M <^ 5 -+*He (eas + e (AIV.13) Consider next, the integral J^. = w T3P J. d 5(AIV.14).J This may be readily integrated; it turns out to be ^? X+t [S at r;" ( 1 J) (AIV.15) Then Jo becomes C _1 NO __ —-+^ - (i- d~ _ JG-i7 s

-106The first integral in (AIV.16) can be performed by using s2 as integration variable. Then QO OtTt 7S 30 I st c3~I+ t The second term in JO can be integrated by using _ as inteJS gration variable. Thus, r -F 2 J o Bt ^-o: i ad I te(r + L + Q "OS (AIV.18) The final integral in Jo cannot be done explicitly, but if we transform it so that the integration limits are -1 and +1, we may apply quadrature formulae quite readily. If we let h4. =; ( > (AIV. 19) where;}

-107the transformation may be carried out. Then T I f t -1 t _ t ( + _ _ _ _ _)^8 1v AtI'; tT T (C- t ~ ~(A IV.20) where rJ( {t) _ i e >-' = O; r H 2 n (AIV.21) The integral J(j ),t) is tabulated in Appendix VI according to the results of a five-point Radau formula. Finally, we have to consider the integral Th This last integration involves no special difficulty and we obtain [nrT1~.. _____ ek I, s (AIV.22) &f- ~$, I —T ~+ o,'3.

-108Then, on gathering the results of the integration from Equations (AIV.10), (AIV.20) and (AIV.22), we have 7' I I -T ~ I+ 0qj 2- I r+ t - ~ \| / it c| _A i 1 1 "52t^'; A -;| -t < t — Ad?^ (Gin

APPENDIX V The Calculation of X2 We have again to evaluate a number of integrals to find \^- _~^__.- ___ ~,,.................... oi L2W, dS (AV.1) where r^:(I) = It L;XC (94) and -(0-02.)-,\ T= (AV.2) We repeat the result of Appendix II, namely, w -&1'-Is ~IP~'I gd^-^ - i l) k- ds (AII.15) First, we solve Equation (AV.2) for W2. We find Vw = Aw$aS + 6B f + C O 6a + I +o3) -109t-1 CA5.4, I

-110where A, B, C, D are as yet arbitrary and q2= a2 + (2it)2. The conditions to be satisfied by W2 are -= _ -~ —. ~ (AV.4) If we take only even solutions and use these conditions, we learn that A = - G < i i ( -) Aac4 (a+ & ) (Av.5) and h> q ( __Qc_ ) (AV.6) B and C are zero. The first integral is'i.L "^~J Jl'-t >T|) i1 - 3 (AV.7) Then, we must consider I L Z d

-11Il which may be rewritten aIt < ^ 3 T (t 15 + This is equivalent to ~__f ds | i-f-i+ Qt'sxL'+,::Asi \ -CLS I *o&^^^ o&^ltn^ ) Ed (AV.10) Consider (AV.12) t+1&s-' which may be rewritten as = ~ tS Cser

-112Thus sJ (2,u-L ( 1 Fe ( (AV.13) With Formula (AV.13) and some manipulations, we find Z~-=.1~~~~~ 3 -(a(SL,cS (- + ) C2 C) ) e TC. F tS)'t4 (~,Tt] After considerable calculation, this becomes ^ / ^ \Cr- It- Try 4I Tk caa t) ^ (Av.15) W1L- _c,,:

-113where +| a 41~ CL~~~~~t Vr(4) ckx~G~ - ~~l't721Th VTO'J --— f/ \ J^^I" (Av.16) Also, fT1T -crg - 1 ^^^ S~lr'tca~"^01 (Av.17) Hence, we obtain \, 2^f %1 Rc-' 3~ Jt-^ ^ 1 ^^ ^ k — + l - Lt~ ~ ~ ~ ~ C IT + 127 \4(:'n7 Y -1 ^irLq-^ L <x^(^Aoi~~~~i~a)~ ) 7' I

-114APPENDIX VI Tabulation of TJ (Qt)) As Calculated _ By a Five-Point Radau Formula a \ 0.1 0.5 1.0 5.0 1.0 3.656 (-3) 2.898 (-3) 1.818 (-3) 5.99 (-6) 2.0 1.204 (-3) 9.27 (-4) 6.28 (-4) 1.84 2.5 5.83 (-4) 4.71 2.28 (-4) 1.38 3.0 283 2.33 1.69 1.00 (-6) 3.1 2.44 2.01 1.48 9.3 (-7) 3.3 1.82 1.51 1.24 (-4) 8.1 3.5 1.35 1.13 (-4) 8.5 (-5) 6.8 3.7 1.02 (-4) 8.5 (-5) 6.4 5.8 4.0 6.5 (-5) "5.4 4.2 4.6 4.5. 3.1 1.9 1.7 3.0 5.0 1.5 (-5) 1.3 (-5) 1.0 (-5) 1.9 (-7) 6.0 4 (-6) 3 (-6) 2 (-6) 1 (-8) Numbers in parentheses indicate powers of ten to be included as factors.

BIBLIOGRAPHY Allen, C. W. 1955 Astrophysical Quantites, London, University of London, Athlone Press. Aller, L. H. 1953 Astro hysics, New York, Ronald Press. Babcock, H. W. 1958 A. J., III, 141. Benard, H. 1900 Rev Geno Sci. Pur. App., 12, 1261, 1309. Biermann, L. 1937 Ast. Nach., 264, 395. 1942 Zts. f. Ap., 21, 320. 1948 Zts. f. Ap., 25, 135. Bless, R. C. 1958 University of Michigan Thesis. Boussinesq, J. 1903 Theorie Analytiue de la Chaleur, Tome II, Paris. Chalonge, D. and 1952 Ann. d1 Ap., 15 201. Divan, L. Chandrasekhar, S. 1952 Proc. Roy. Soc., A, 217, 306. 1955 Proc. Camb. Phil. Soc., 51, 162. Courant, R. and 1953 Methods of Theoretical Physics, Vol. I, Hilbert, D. New York, Interscience. Cowling, T. G. 1951 Ap. J., 114, 272. Eddington, A. S. 1930 The Internal Constitution of the Stars, Cambridge, University Press. Edmonds, F. N. 1957 A.., 125, 535. Fan, T. Y. T. 1955 Ap. J., 121, 508. Goldstein, S. 1938 Modern Developments in Fluid Mechanics, Oxford University Press. Goody, R. M. 1956 Jo Fluid Mech., 1, 424. Gribov, V. N. and 1957 Phys. Abstracts, No. 3045. Gurevich, L. E. Hack, M. 1956 BAN, 13, No. 466. -15

-116BIBLIOGRAPHY (CONT' D) Hide, R. 1955 Proc, Camb, Phil, Soc., 51, 179* Hunger, K. 1955 Zts, f Ap,, 36, 42, Jeffreys, H. 1926 Phil. Mag,, 2, 833* 1928 Proc. Roy, Soc., A, 118, 195. 1930 Proc. Camb. Phil. Soc., 26, 170. Jugaku, J. 1958 University of Michigan Thesis. Lamb, H. 1890 Proc. Roy. Soc,. A, 84, 551. 1931 Hydrodynamics, 6th Ed., Cambridge, Low, A, R. 1928 Proc. Roy. Soc., A, 125, 180o Magnus, W. and 1954 Formulas and Theorems for the Functions of Oberhettinger, F. Mathematical Physics, New York, Chelsea Pub. Co. Malkus, W. V. R, 1954 a Proc. Roy. Soc., A, 225, 196. 1954 b Proc. Roy, Soc., A, 225, 185. Malkus, W. V, R. 1958 Jo Fluid Mech, 4, 225. and Veronis, G, Morse, P. M. and 1953 Methods of Theoretical Physics, Vol. I, Feshbach, H. New York, McGraw-Hill* Osawa, K, 1956 Ap, J., 1235 513* Oster, L. 1958 Zts, f. Ap.* 44, 26, Parker, E. N, 1953 Ap, J^, 117, 431. Pecker, C. 1953 Ann. dA'Ap., 16, 321. Pellew, A, and 1940 Proc, Roy. Soc., A, 176, 312. Southwell, R, V. Prandtl, L. 1952 Essentials of Fluid Dynamics, Chap. V, New York, Hafner Pub. Co. Rayleigh, Lord 1883 Proc. London Math, Soc., 14, 170 (also Scientific Papers, Vol. 2, 200). 1916 Phil, Mag., (6)532, 529.

-117BIBLIOGRAPHY (CONT D) Richardson, R. S. & 1950 Ap.o 111, 351. Schwarzschild, M. Rogerson, Jr., J. B. 1958 Sky and Tel., 52, 112. Rosa, A. and 1948 Zts. f. Ap., 25, 20. Uns old, A. Rudkjobing, M. 1947 Pub, Cop Obs., No. 145. Saito, S. 1954 Contributions from Inst. Ast,, University of Kyoto, No. 48. Schmidt, R. J. and 1935 Proc. Roy. Sc., A, 152, 586. Milverton, S. W. Schmidt, R. J. and 1938 Proc. Roy. Soc., A, 165, 216. Saunders, O. A. Schwarzschild, K. 1906 Gott. Nach, 41. Siedentopf, H. 1932 Astr. Nach., 247, 297. Skumanich, A. 1955 Ap. Jo, 121, 408. Spiegel, E. A. 1957 Apo J., 126, 202. Struve, 0. 1950 Stellar Evolution, Princeton University Press, Princeton. Stuart, J. T. 1958 J. Fluid Mech., in press. Traving, G. 1955 Zts. f. Ap., 36, 1. Ueno, S. and 1950 Pub. Ast. SocO Japan, 2, 32. Matsushima, S. Underhill, A. B. 1950 Pub. Cop. Obs., No. 151. Unno, W. 1957 Ap. J., 126, 259. Unsold, A. 1930 Zts. f. Ap., 1, 138. 1948 Zts. f. ApQ, 25, 11. 1955 PhYSik der Sternatmospharen, Berlin, J. Springero Vitense, E. 1953 ZLtso f Apo, 32, 135. Whitney, C. A. 1955 Stellar Pulsation, Harvard University Thesis.