"Science, Faculty of"@en . "Physics and Astronomy, Department of"@en . "DSpace"@en . "UBCV"@en . "Carlberg, Raymond G."@en . "2010-03-04T23:42:11Z"@en . "1978"@en . "Doctor of Philosophy - PhD"@en . "University of British Columbia"@en . "This thesis investigates the quantitative nature of the variability which is present in the stellar winds of high luminosity early type stars. A program of optical observations with high time and spectral resolution was designed to provide quantitative information on the nature of the fluctuations. These observations found no optical variability over a time period of six hours and hence restrict the variability over this period to size scales of less than 5x1011 cm, but do confirm the variations on time scales exceeding one day. A class of X-ray sources comprised of a neutron star orbiting a star with a strong stellar wind provides another source of information on the variability of stellar winds. A theory of accretion onto a neutron star was developed which is used with X-ray intensity data to derive estimates of the density and velocity of the stellar wind. This analysis performed on Cen X-3 suggests that the velocity in the stellar wind increases as the wind density increases. A theoretical analysis of the stability of a stellar wind is made to determine whether the variability may originate in the wind itself. Two types of instability are founds those that amplify pre-existing disturbances, and absolute instabilities which can grow from random motions within the gas. It is found that short wavelength disturbances (<10\u00E2\u0081\u00B4 cm) are always strongly damped by conduction, and long wavelength ones (>10\u00C2\u00B9\u00C2\u00B9 cm) are damped by radiation if the gas is thermally stable, that is if the net radiative energy loss increases with temperature. Intermediate wavelengths of about 10\u00E2\u0081\u00B8\u00E2\u0081\u00BB\u00E2\u0081\u00B9 cm are usually subject to an amplification due to the density gradient of the wind. The radiation acceleration amplifies disturbances of scales 10\u00E2\u0081\u00B7 to 10\u00C2\u00B9\u00C2\u00B9 ca. Absolute instabilities are present if the gas is thermally unstable, if the flow is deccelerating, or if the gas has a temperature of several million degrees. On the basis of the information derived on stellar wind stability it is proposed that a complete theory should be based on the assumption that the wind is a nonstationary flow."@en . "https://circle.library.ubc.ca/rest/handle/2429/21496?expand=metadata"@en . "RADIATION DRIVEN INSTABILITIES IN STELLAR WINDS by RAYMOND GARY CARLBERG M.Sc. , University of B r i t i s h Columbia, 1975 B.Sc., University of Saskatchewan, 1972 A THESIS SUBMITTED IN PARTIAL FULFILLMENT OF THE REQUIREMENTS FOB THE DEGREE OF DOCTOE OF PHILOSOPHY in THE FACULTY OF GRADUATE STUDIES* Department of Geophysics and Astronomy We accept t h i s thesis as conforming to the reguired standard THE UNIVERSITY OF BRITISH COLUMBIA October, 1978 (c) Baymond Gary Carlberg, 1978 In presenting th i s thes is in pa r t i a l fu l f i lment of the requirements for an advanced degree at the Univers i ty of B r i t i s h Columbia, I agree that the L ibrary sha l l make it f ree ly ava i lab le for reference and study. I fur ther agree that permission for extensive copying of th i s thesis for scho lar ly purposes may be granted by the Head of my Department or by his representat ives. It is understood that copying or pub l i ca t ion of th i s thes is fo r f inanc ia l gain sha l l not be al1 owed without my wr i t ten permission. Department of The Univers i ty of B r i t i s h Columbia 2075 Wesbrook P l a c e Vancouver, Canada V6T 1W5 Date Odf *? jl% i i ABSTRACT This t h e s i s investigates the quantitative nature of the v a r i a b i l i t y which i s present i n the s t e l l a r winds of high l u -minosity early type st a r s . A program of o p t i c a l observations with high time and spectral resolution was designed to provide quantitative information on the nature of the fluctuations. These observations found no op t i c a l v a r i a b i l i t y over a time period of s i x hours and hence r e s t r i c t the v a r i a b i l i t y over t h i s pericd to size scales of less than 5x10 1 1 cm, but do confirm the variations on time scales exceeding one day. A class of X-ray sources comprised of a neutron star o r b i t i n g a star with a strong s t e l l a r wind provides another source of information on the v a r i a b i l i t y of s t e l l a r winds. A theory of accretion onto a neutron star was developed which i s used with X-ray i n t e n s i t y data to derive estimates of the density and vel o c i t y of the s t e l l a r wind. This analysis performed on Cen X-3 suggests that the velocity i n the s t e l l a r wind increases as the wind density increases. A t h e o r e t i c a l analysis of the s t a b i l i t y of a s t e l l a r wind i s made to determine whether the v a r i a b i l i t y may originate i n the wind i t s e l f . Two types of i n s t a b i l i t y are founds those that amplify pre-existing disturbances, and absolute i n s t a b i l i t i e s which can grow from random motions within the gas. I t i s found that short wavelength disturbances <<10* cm) are always strongly damped by conduction, and leng wavelength ones (>10 1 1 cm) are damped by radiatio n i f the gas i s thermally stable, that i s i f the net radiative energy loss increases with temperature. In-termediate wavelengths of about 10 8~ 9 cm are usually subject to i i i an amplification due to the density gradient of the wind. The radiation acceleration amplifies disturbances of scales 10 7 to 10*1 ca. Absolute i n s t a b i l i t i e s are present i f the gas i s ther-mally unstable, i f the flow i s deccelerating, or i f the gas has a temperature of several m i l l i o n degrees. On the basis of the information derived on s t e l l a r wind s t a b i l i t y i t i s proposed that a complete theory should be based on the assumption that the wind i s a nonstationary flow. CONTENTS Chapter 1: Introduction .................................... Chapter 2: Optical Observations Chapter 3: Supersonic Accretion Chapter 4: Physical Description of The Gas ................. Chapter 5: The S t a b i l i t y analysis .......................... Chapter 6: Conclusions ..................................... Bibliography Appendix 1: Radiative Effects In Supersonic Accretion ...... Appendix 2: Gas Physics .................................... Appendix 3: The Dispersion Relation ........................ Appendix 4: The Major Computer Prcgrams .................... V FIGURES 1. Lambda C e p h e i : The E f f e c t Of R e s o l u t i o n ................ 10 2. lambda Cephei Time S e r i e s .............................. 14 3. lambda Cephei Day Tc Bay ...15 4. Alpha C a m e l o p a r d a l i s Day To Day ........................ 16 5. D e l t a O r i o n i s Day To Day ...17 6. S u p e r s o n i c A c c r e t i o n Schematic .........................22 7. Schematic X-ray I n t e n s i t y V a r i a t i o n of Cen X-3 .........27 8. The D e n s i t y V e l o c i t y V a r i a t i o n o f Can X-3 .............. 28 9. I o n i z a t i o n 8 a l a nc \u00E2\u0082\u00AC ..........\u00C2\u00BB\u00C2\u00AB .... ..........\u00C2\u00AB\u00E2\u0080\u00A2........ .35 10. Standard H e a t i n g And C o o l i n g Rate ..................... 36 11. CNO Abundances 10 Times S o l a r H e a t i n g And C o d i n g .....38 12. L o s s e s I n An O p t i c a l l y T h i c k Medium ................... 40 13. R a d i a t i o n F o r c e As A f u n c t i o n Of Temperature ..........42 14. Momentum B a l a n c e ...................................... 46 15. Roots For A S t a t i c Pseudo I s c t h e r m a l Atmosphere ....... 55 16. H i t h H e a t i n g And C o o l i n g , V, dv/dz Nonzero ........... .60 17. 8 i t h The R a d i a t i o n F o r c e ..............................62 18 \u00E2\u0080\u00A2 No C o o X x n c j \u00E2\u0080\u00A2 \u00E2\u0080\u00A2 \u00E2\u0080\u00A2 \u00E2\u0080\u00A2 \u00E2\u0080\u00A2 \u00E2\u0080\u00A2 \u00E2\u0080\u00A2 \u00E2\u0080\u00A2 \u00E2\u0080\u00A2 * \u00E2\u0080\u00A2 * \u00E2\u0080\u00A2 \u00E2\u0080\u00A2 \u00E2\u0080\u00A2 \u00E2\u0080\u00A2 * \u00E2\u0080\u00A2 \u00E2\u0080\u00A2 \u00E2\u0080\u00A2 \u00E2\u0080\u00A2 * \u00E2\u0080\u00A2 \u00E2\u0080\u00A2 \u00E2\u0080\u00A2 \u00E2\u0080\u00A2 * \u00E2\u0080\u00A2 \u00C2\u00AB \u00E2\u0080\u00A2 \u00E2\u0080\u00A2 \u00E2\u0080\u00A2 \u00E2\u0080\u00A2 * \u00E2\u0080\u00A2 \u00E2\u0080\u00A2 \u00E2\u0080\u00A2 \u00E2\u0080\u00A2 63 1. Thermal I n s t a b i l i t y T=5x10s K ......................... .65 19. High Temperature I n s t a b i l i t y .......................... 66 20. D e c c e l e r a t i o n I n s t a b i l i t y .............................67 v i TABLES 1. Catalogue Of Observations 9 2. S c a l e s In Supersonic A c c r e t i o n ........................ .23 3. Atomic Abundances ......................................32 4. Maximum V e l o c i t y For An A c c e l e r a t i n g S o l u t i o n ..........47 5. The D i s p e r s i o n R e l a t i o n s P l o t t e d .......................68 6. P h o t o i o n i z a t i o n Cross S e c t i o n Parameters ............... 90 7. Gaunt F a c t o r For Hydrogen 92 8., The Recombination Hate Constants .......................93 9. I o n i z a t i o n P o t e n t i a l s And S u b s h e l l P o p u l a t i o n s .........97 10. Gaunt Factor Constants ................................ 102 11. The Resonance L i n e s ................................... 103 v i i ACKNOWLEDGEMENTS The completion of t h i s thesis owes a great deal to the teaching, assistance, and encouragement that I have received. My supervisor. Dr. Greg Fahlman, always provided an interested ear and a careful c r i t i c i s m of my ideas. He was also a scurce of support when needed. Drs. Gordon Walker and Jason Auman made a number of suggestions which substantially improved the presentation of the thesis. My induction into the world of observational astronomy came with three very pleasant summers spent in V i c t o r i a at the Dominion Astrophysical Observatory. . The fr u s t r a t i o n s of doing spectrophotometry with photographic plates were largely elimin-ated by Dr. Gordon Walker's encouragement to use his fieticon detector system. The successful observations obtained with t h i s system owe a great deal to the advice and technical support of Dr. Walker, Dr. Bruce Campbell, Dr. Chris P r i t c h e t t , Tim Lester, and Mike Creswell. The whcle data reduction process has been reduced to a straightforward and enjoyable a c t i v i t y by the use of Chris P r i t c h e t t * s data handling package, RETICENT. A large part of t h i s thesis r e l i e d on the computer for i t s completion. The OBC Computing Centre has been a great a s s i s -tance and deserves much praise for i t s f a c i l i t i e s , and the pro-vision of a large number of well documented u t i l i t y programs. Pa r t i c u l a r thanks go to A. C. Hearn, o r i g i n a l l y at the University of Utah, who wrote the program REDUCE. Financial support was provided by the National Research Council of Canada and by a MacMillan Family Fellowship. 1 CHAPT IB 1. INTE0DUC1I0N The o p t i c a l spectra of many hot stars were noted to have emission components on the Harvard objective prism plates {see the Henry Draper Catalogue, Cannon and Pickering 1918). Many of these stars were studied in d e t a i l , leading to Beal's (1929) proposal that the l i n e p r o f i l e s could be explained by emission from gas being ejected from the star. A great deal of informa-tion was accumulated on the o p t i c a l spectra of emission l i n e 0 and E stars i n the following years, which i s summarized i n Eeals (1951) and Dnderhill (1960). A new observational window^ was opened by Morton (1967) using a rocket borne spectrograph. He found that 4 of the Orion s t a r s , S , 6, and \u00C2\u00A3., had absorption l i n e s blue s h i f t e d to v e l o c i t i e s of order 2000 Km s _ l . These v e l o c i t i e s exceeded the t y p i c a l escape speed of 300 Km s ~ 1 by such a large factor that there was no question that the stars were losing mass at a large rate. The Snow and Morton u l t r a v i o l e t survey (1976) showed that a l l stars with an e f f e c t i v e temperature greater than about 3x10* K, and luminosities greater than a bolometric magnitude of -6, have a detectable supersonic wind which c a r r i e s away a s i g n i f i -cant amount of the star's mass during i t s l i f e t i m e . Reliable mass loss rates have been obtained from o p t i c a l observations (see Hutchings 1976) and OV observations (Snow and Morton 1976) which now have been extended to infrared (Barlow and Cohen 1977) and radio wavelengths (Bright and Barlow 1975), a l l of which provide confirming and complementary data on the magnitude of the mass loss . The derived mass loss rates for OB stars l i e i n the range of 10~ 9 to 10~ s M0/year, with terminal v e l o c i t i e s 2 ranging from 1000 to 3000 km s-\u00C2\u00BB. The s t e l l a r wind phenomenon poses several questions: what i s the physical mechanism driving the mass l o s s ? how does the mass loss e f f e c t the star's evolution? and how do these hot, luminous stars - % f feet the i n t e r s t e l l a r medium and the evolution of a galaxy? The answers to a l l cf these questions hinge on a thorough understanding of the physics of the s t e l l a r wind. This thesis i s a contribution to the understanding of the basic nature of the s t e l l a r wind. The physical problem i s to describe the dynamics of a qas moving i n an intense radiation f i e l d ; a situation which occurs i n a number of astrophysical situations, including quasars and active g a l a c t i c n u c l e i , (Mushotzky, et a l . 1972 and Kippenhahn et a l . 1975). The l i n e p r o f i l e s cf s t e l l a r wind stars, especially ones with high mass loss rates have long been known to show some v a r i a b i l i t y over one day (see for instance Beals 1951, Underbill 1960, Conti and Frost 1974, Leep and Conti 1978, Brucato 1971, Snow 1S77, and Hosendahl 1S73). For ground based observation t h i s time period makes i t d i f f i c u l t to resolve the time evolu-tion of the vari a t i o n . An observational program was i n i t i a t e d to more cl o s e l y define the nature of these reported variations, using a modern detector capable of measuring very small changes. This w i l l be discussed i n Chapter 2. Besides o p t i c a l evidence of v a r i a b i l i t y , there are a number of X-ray binary stars where the X-ray source i s a neutron star accreting mass from the s t e l l a r wind fCcnti 1978). These sources show a number of scales of variation of the i r i n t e n s i t y which can be ascribed to variations of the s t e l l a r wind. A 3 theoretic a l analysis of the accretion process and how i t e f f e c t s the ctserved intensity was made i n order that the X-ray data could be used to derive the prevailing density and vel o c i t y of the s t e l l a r wind at the location of the neutron s t a r . This an-a l y s i s performed on X-ray data for the source Cen X-3 indicates a c o r r e l a t i o n between the wind velocity and density. This w i l l be described in Chapter 3. The basic formulation of the theory of a s t e l l a r wind from a hot, luminous star was i n i t i a l l y put forth by Lucy and Solomon (1970), who proposed that the acceleration was produced by the scattering of photons with wavelengths that f e l l within a few resonance l i n e s . This was a generalization of Milne*s (1926) idea that momentum transfer from photons could s e l e c t i v e l y acce-lerate certain ions. This was l a t e r extended by Castor, Abbott, and Klein (1975, referred to as CAK) to include the force on many l i n e s of many ions. The theory provided an encouraqinq agreement with the limited data available on the velocity as a function of radius and mass loss rates. Recently the u l t r a v i o l e t s a t e l l i t e observations have re-vealed that some highly ionized species, i n par t i c u l a r 0 VI and N V, are present i n the wind. These ions would not be expected to be ionized i n any observable quantity by the radiation f i e l d appropriate to these stars. York, j t a l . (1977) have observed variations i n the 0 VI line i n three stars over a time periods as short as s i x hours. This observation suggests a \"slab\" moving outwards at an increasing v e l o c i t y . . The presence of 0 VI in the s t e l l a r wind presents a puzzle as to the source of i t s excitation. At the present time there are three proposals. 4 F i r s t , Castor (1978) has modified his radiation driven wind to an a r b i t r a r i l y s p e c i f i e d temperature higher than radiative equi-librium, which provides a suitable abundance of C VI. Second, Lamers and Snow (1978) have an empirical \"warm radiation pressure\" model, in which they show that the ions can be provi-ded i f the s t e l l a r wind i s at a temperature of aJsout 2x10 s K* Neither of these models specify the source of the additional heating. Third, Hearn (1975) has proposed that s t e l l a r winds are i n i t i a l l y accelerated i n a hot corona with a temperature cf several m i l l i o n degrees. Pursuing t h i s idea Clson and C a s s i n e l l i (1978) have shown that a small corona, about 10% of a s t e l l a r radius, generates enough thermal X-rays to produce the reguired i o n i z a t i o n r a t i o s . To provide a heating mechanism for a corona, Hearn (1972) showed that radiation driven sound waves could be amplified while propogating outward i n the atmosphere. The waves grow to a saturated amplitude s u f f i c i e n t to provide enough shock heating to maintain a corona (Hearn 1973). There are two d i f f i c u l t i e s with t h i s analysis. Berthomieu et a l . (1975) have pointed out that Hearn*s simplifying assumptions resu l t i n a scale length for the wave amplification which i s the same as the atmospheric scale length. Therefore s i g n i f i c a n t amplification only occurs over lengths which invalidate the as-sumption cf small variations of the zero order quantities over the length for amplification. In addition, the unstable waves that he finds are amplifying i n s t a b i l i t i e s (Castor 1977), and require some o s c i l l a t o r to i n i t i a t e the wave motion. Motivated by t h e o r e t i c a l arguments and the observations of fluctuations I have performed a s t a b i l i t y analysis on the egua-5 tions governing the moving gas i n the s t e l l a r radiation f i e l d . The complete set of eguaticns governing the motion with no a \u00C2\u00A3\u00C2\u00A3ipri s i m p l i f i c a t i o n s were used. An accurate description of the gas physics was developed using an approximate treatment of radiation transfer dependent ccly on l o c a l guantities. As a result the state of the gas can be completely s p e c i f i e d by the l o c a l radiation f i e l d , the gas velocity and i t s gradient, and the density and temperature, as described in chapter 4. The s t a b i l i t y of the gas against v e r t i c a l disturbances was i n v e s t i -gated with the aid of a computer to provide the numerical solu-tions to the dispersion r e l a t i o n . Chapter 5 i s comprised of t h i s discussion. The purpose of t h i s thesis i s to investigate the quantita-tiv e nature of i n s t a b i l i t i e s i n s t e l l a r winds and r e l a t e i t to the observational and t h e o r e t i c a l problems which have been out-l i n e d . This i s not an attempt to create a unified theory of a s t e l l a r wind. Bather i t i s a detailed investigation of ce r t a i n areas of the question i n order to illuminate some of the physi-c a l mechanisms which are important i n a s t e l l a r wind. This i s required because not much i s known about the basic physical pro-cesses which dominate the observed v a r i a b i l i t y of the s t e l l a r wind. The investigation i s confined to the s t e l l a r wind i t s e l f , which i s loosely defined as the region where the o p t i c a l depth in tbe continuum i s les s than one and the gas i s moving with greater than scnic v e l o c i t i e s . As has been emphasized by Cannon and Thomas (1978), i t i s possible that some of the driving force for the the wind and hence seme of the wind i n s t a b i l i t i e s may 6 originate within deeper layers of the star. It i s assumed that there i s no magnetic f i e l d . This i s done mostly because of the tremendous s i m p l i f i c a t i o n of the problem which r e s u l t s . But there i s no observational evidence for a magnetic f i e l d , although i f the wind i s as chaotic as t h i s thesis suggests, a magnetic f i e l d would be d i f f i c u l t tc detect. In summary t h i s thesis i s motivated by observations of s t e l l a r wind v a r i a b i l i t y , and suggestions by ether authors that i n s t a b i l i t i e s do exis t which may be responsible for the creation of a high temperature corona. The investigations described are carried out i n twe parts. observational evidence of the v a r i a -b i l i t y i s acquired which suggests length and times scales cf the fluctuations which are present, and a c o r r e l a t i o n between the wind ve l o c i t y and density. The t h e o r e t i c a l analysis provides physical sources of several i n s t a b i l i t i e s which can exist i n the s t e l l a r wind. From t h i s information I suggest that the s t e l l a r wind i s an extremely chaotic medium i n which the i n s t a b i l i t i e s not only provide the source for the observed v a r i a b i l i t y , but also can be used to provide an ion i z a t i o n source for the 0 VI ion and the corona as postulated by Hearn. The presence of these i n s t a b i l i t i e s means that a model for a s t e l l a r wind should be i n the form of mean flow guantities and associated f l u c t u -ating guantities. 7 CHAPTER 2. OPTICAL OBSERVATIONS For many years several of the l i n e s i n the o p t i c a l spectrum of several early type stars have been reported as varying (see references c i t e d in the Introduction). P a r t i c u l a r attention has been paid to the star Lambda Cephei, because i t i s a bright 06f star i n the northern sky. The H \u00C2\u00A9cline has been reported to vary on time scales of one day (Leep and Conti 1978) and longer, with no apparent systematic variation. The amplitude of the v a r i a -tion i s t y p i c a l l y 10% of the in t e n s i t y . This behaviour i s f a i r l y t y p i c a l of the more luminous mass loss s t a r s . The shor-test period variations with a high confidence l e v e l are the UV observations made by the s a t e l l i t e Copernicus of S o r i A, T j Q r i , and \u00C2\u00A3 Pup (York et a l . 1977), where a small feature of width about 150 Km s - 1 was seen to \"move\" i n the 0 VI l i n e between two observations spaced about 6 hours apart. Host of the observations at o p t i c a l wavelengths have been made with photographic plates, which have a photometric accuracy barely able to reveal the presence of the v a r i a t i o n , l e t alone reveal much information as to i t s character. In fact Lacy (1977) made scanner observations of some of the li n e s i n stars that were reported as varying and concluded on the basis cf a s t a t i s t i c a l analysis of the errors present i n the eguivalent width that any v a r i a b i l i t y present was less than the expected randcm error. However the eguivalent width of a l i n e averages together a l l material emitting at that l i n e freguency,, Observa-tions which resolve the line can provide much more information, but at the cost of longer exposure times. The c l a s s i c a l description of l i n e formation i n a s t e l l a r 8 wind was given by Beals (1951). The observed l i n e can be consi-dered to be made up of three almost independent parts; an under-lying abscrpti.cn l i n e formed in the photosphere of the star, a superposed emission l i n e with i t s centroid at zero velocity pro-duced by emission of photons in the s t e l l a r wind, and a blue shifted absorption l i n e which i s formed i n the portion of the s t e l l a r wind which i s silhouetted against the star. The analysis of l i n e formation within the wind i s si m p l i -fied by the Sobolev approximation (Sobolev 1960), which says that the emission and absorption of photons i n a given narrow wavelength i n t e r v a l , outside of the doppler core, i s determined by the amount of gas moving at a velocity such that the l i n e of sight velocity of the gas f a l l s within the wavelength i n t e r v a l . This approximation i s va l i d i f the gas speed i s supersonic. The assumption i s supported by the observations of Hutchings (1976 and references therein) who has shown that the wind has a velo-c i t y exceeding the sound speed for distances greater than 10$ of the s t e l l a r radius. The observations were undertaken to confirm the reported v a r i a b i l i t y and were to be made with s u f f i c i e n t l y high signal tc noise, spectral resolution and time resolution to c l e a r l y r e -solve the variations, as they developed. In p a r t i c u l a r i t was thought that there might be evidence f o r the nature of the mech-anism of the variation, for instance, a spot rotating with the star cr a \"blob\" moving cut through the wind. fill observations were carried out with the 1.2 meter t e l e s -cope of the Dominion Astrophysical Observatory, V i c t o r i a , E.C. The 2.4 meter camera i n the ccude spectrograph was used with a 9 red coated image s l i c e r giving a projected s l i t width cf 6C mi-crons. The spectrum was detected with a 10 24 element array of 25.4 micron diodes (a Beticon RL1024/C17) cooled to a tempera-ture of -80\u00C2\u00B0 C (Walker et a l . 1S76). The image s l i c e r and de-tector pixel size combination were chosen to give a properly oversampled spectrum. , A l l observations were centred on the Ho< l i n e in the f i r s t order r e s u l t i n g i n a dispersion of .125A/diode. Observations were made in September and October of 1977, and are tabulated below and shown i n the accompanying f i g -ures. TABLE 1: CATALOGUE OF OESIBVATIONS # Star Date Time Expos 1977 PSI secon 1 lambda Cep Sept 11/12 22: 13 2250 2 i i \u00E2\u0080\u00A2 i 22:55 u 3 \u00C2\u00AB it 23:41 tt 4 t i t i 00:23 i i 5 n \u00E2\u0080\u00A2 i 01:05 i i 6 t i n 01:48 t i 7 II i i 03:10 t i 8 n tt 03:54 i i 9 \u00C2\u00BB II 04: 36 i i 10 i t Oct 11/12 23:00 3000 11 t i n 01:08 30 00 12 t i Oct 12/13 22:15 3000 13 H Oct 16/17 21:05 3000 14 \u00E2\u0080\u00A2t II 23:38 30 00 15 t i tt 04:28 3000 16 Alpha Cam Oct 11/12 00:27 1500 17 i i Oct 12/13 01:35 2002 18 II Oct 16/17 22:58 1800 19 Delta Ori Oct 12/13 02:31 600 20 \u00E2\u0080\u00A2 i tt 03:01 600 21 t i Oct 16/17 23:00 2641 22 i t it 03:41 2830 The l i n e s present i n the 100 A region examined are i d e n t i -f i e d i n Figure 1., They include the s t e l l a r and He I I 6527, the i n t e r s t e l l a r 6614 A feature, and a multitude of narrow. 10 Fig 1: Lambda Cephei: The Effect of Resolution 11 weak, t e l l u r i c water vapour l i n e s . The t e l l u r i c water l i n e s and t h e i r r e l a t i v e equivalent widths are indicated above the spec-tra. This data was taken from Moore et a l . (1966), and may not give the exact r e l a t i v e i n t e n s i t i e s for these observations. A l l cf the spectra have been f i l t e r e d by a Fourier transform techni-que to U0% of the Nyquist frequency, which i s roughly the true resolution of the spectra,, A l l spectra had considerable (10$) response changes along the array, due mostly to a l i g h t f r o s t on the window of the detector. This was removed by dividing by the spectrum of a lamp which was taken immediately before cr af t e r the observation. For the time series spectra the underlyinq shape of the spectrum did not vary within the error (0. }%) of the lamp c a l i b r a t i o n . A l l spectra were r e c t i f i e d to a l i n e a r continuum, Fiqure 1 shows the absolute necessity to resolve the t e l -l u r i c water l i n e s . As the Fiqure 2 time series of Lamda Cep over 6 hours shows, the water l i n e s vary s i g n i f i c a n t l y over one hour. In Figure 1 the top spectrum shows the mean of the time series of high resolution spectra. Below i t are two spectra recorded at KPNO in June, 1978 (cocrtesy of G. G, Fahlman and G. A. H. Salker) using a lower resolution spectrograph. The bottom spectrum in Figure 1, i s the top spectrum but convolved with a Gaussian to give approximately the same instrumental re-solution as the KfHO spectra. I t i s evident that the variation in an H*, p r o f i l e can be e n t i r e l y due to t e l l u r i c water l i n e variations, i f the instrumental resolution i s inadequate to c l e a r l y separate these variations out. The time series of Lambda Cephei (spectral type 06f) shown 12 in Figure 2 covers 6.5 hours. The average of t h i s time series of observations i s shown at the top of Figure 1. The l i n e s below are the i n d i v i d u a l spectra divided by the mean, then nor-malized. Although there are suggestions cf underlying broad (say about 10 A) changes, these are less than the noise l e v e l . In the day tc day observations shown i n Figure 3, there i s clear evidence of a variation at the Hotline of the emission feature at v e l o c i t i e s near 200 Km s _ l , and on the absorption side at v e l o c i t i e s near -300 Km s~*. The time series difference spectra, number 1 to 9 of Figure 2, can be analyzed to determine the s t a t i s t i c a l significance of any variations. The report by York et a l . (1977) of a feature of FwHM 150 Km s~ 1 changing over a period of 6 hours i s only s l i g h t l y wider than seme of the t e l l u r i c water features, and leads to some d i f f i c u l t l y i n interpreting changes. The standard deviation of the spectra i s i n the range of 0.6 to 0.8% of the mean, other than for spectrum 9. Assuming the noise to have a normal d i s t r i b u t i o n with t h i s variance, the fluctuations must have an amplitude exceeding 2.57 standard deviations to have a less than 1% probability of chance occurence. This amplitude i s indicated in Figure 2., The smoothed serie s of plots i n Figure 2 are the same spectra as those on the l e f t but averaged over 11 diodes. This reduces the variance by a factor of the square root of 11. The l i n e s for a s t a t i s t i c a l s i g n i f i c a n c e cf 99% are again drawn on the plot. I t can be seen that there are many features which do vary s i g n i f i c a n t l y . But the features that are varying a l l correspond to the wavelengths of t e l l u r i c l i n e s , except for the feature at a velocity with respect to the Hot, l i n e 13 of +200 Km s- 1 ( l e f t dotted l i n e ) . This exceeds the 99% s i g n i -ficance l e v e l i n records 1,2,3,6,7, and 8, going from an excess to a deficiency with respect to the mean. The variation occurs i n the l i n e (see Figure 1) near the top of the emission feature. The subtraction would be very sensitive to very small s h i f t s of the l i n e in thi s region. There are two reasons to think that this feature may not be s t e l l a r i n o r i g i n . F i r s t , i t s variation correlates very well with the water li n e s at a velocity with respect to H of -700 Km s~ 1 (dotted l i n e on r i g h t ) . And sec-ond, the feature shows no velocity s h i f t over t h i s time period, which might be expected i n a wind. I conclude that r e a l v a r i a -tions are present i n the time series, but they are most l i k e l y due to t e l l u r i c features. Alpha Cam (spectral type 09.51a) was chosen because cf spe-c t r a l type, and the presence of the emission l i n e at H$(. Of a l l the stars examined i t seems to have the most s i g n i f i c a n t v a r i a -t i o n s , see Figure 4. Delta Orionis (spectral type 09.511) was observed because of the variation reported by York, et a l . , (1977). Observations made within one night, Oct 16/17, have no r e a l i n d i c a t i o n cf a change. There i s only weak evidence for a p r o f i l e change i n fi v e days, because of the confusion created by the different strength of the t e l l u r i c l i n e s . This star i s a spectroscopic binary of period of 5 days, which produces velocity s h i f t s , but probably not p r o f i l e changes. This i s shown i n Figure 5. The Sobclev approximation allows an estimate of the s i z e of the region producing the photons i n a given wavelength i n t e r v a l . Although the thickness of the s h e l l s of equal l i n e of sight ve-JiJ_JJ.j_ 1000 4?0 -1000 I A ' I I I D i f f e r e n c e Spectra Smoothed over 11 p o i n t s 10 Oct 11/12 23:00 11 Oct 11/12 01:08 21:05 1000 500 0 -500 -1000 Km s' 15 12 Oct 12/13 '\Jty&l 14 Oct 16/17 23:38 ytfl 15 Oct 16/17 04:28 1 ni f i q . 3: Lambda Cephei day t c day 16 Fig. 4: Alpha Cam day to day 17 Fig. 5: D e l t a Ori day t o day 18 l o c i t y w i l l vary w i t h the v e l o c i t y and d i s t a n c e from the s t a r , the c o n t r i b u t i o n t o t h e p r o f i l e w i l l be weighted towards r e g i o n s of h i g h e r d e n s i t y . To e s t i m a t e the t o t a l i n t e n s i t y i n a wave-l e n g t h i n t e r v a l , , a l l i n t e g r a l s o f the e m i s s i o n over the s h e l l w i l l be r e p l a c e d by average q u a n t i t i e s . I f the s h e l l has an average l i n e o f s i g h t t h i c k n e s s . As, t h e n t h e t o t a l volume e m i t t i n g i n t h e wavelength i n t e r v a l i s a p p r o x i m a t e l y 4 t r r ^ 2 A s f , where f i s a f a c t o r c o n t a i n i n g t h e d i f f e r e n c e between the t r u e e m i t t i n g a r e a and the assumed d i s k o f two s t e l l a r r a d i i , 2z$. The c o n s t a n t f w i l l be assumed t o be 1. The t h i c k n e s s o f t h e c o n s t a n t l i n e o f s i g h t v e l o c i t y s h e l l s i s a p p r o x i m a t e l y c o n s t a n t over the r e g i o n of dominant e m i s s i o n . These assumptions a r e j u s t i f i e d by t h e a c t u a l c a l c u l a t i o n s o f l i n e p r o f i l e s as done by C a s s i n e l l i , e t a l . (1978). : The e r r o r of these a p p r o x i m a t i o n s may be as l a r g e as a f a c t o r o f f i v e , but depends on t h e l i n e c o n s i d e r e d . The average s h e l l t h i c k n e s s , A s , w i l l be e s t i m a t e d frcm t h e l i n e c f s i g h t v e l o c i t y g r a d i e n t , dv/ds (=cos 2# dv/dr + s i n 2 ^ v / r , where $ i s the a n g l e between t h e l i n e o f s i g h t and t h e s t a r ) as A f l u c t u a t i o n i n t h e wind which changes the e m i s s i o n r a t e w i l l be observed as some f r a c t i o n a l change of t h e f l u x o b t a i n e d by i n t e g r a t i n g o ver a l l r e g i o n s e m i t t i n g i n t h a t wavelength range. Assuming t h a t t h e e m i s s i o n r a t e changes by 100%, i n t h e f o l l o w -i n g s e c t i o n an e s t i m a t e w i l l be made of t h e s i z e of the f l u c t u a -t i o n c a u s i n g a g i v e n f r a c t i o n a l i n t e n s i t y change. I f the f l u c -t u a t i o n i s a r e g i o n of s i z e jf, c n l y t h a t p a r t of the f l u c t u a t i o n ( D 19 which i s moving at an appropriate velocity to e f f e c t the inten-s i t y in the wavelength i n t e r v a l contributes to the in t e n s i t y change. If the fluctuation i s moving at a uniform v e l o c i t y with an int e r n a l velocity dispersion less than the thermal speed, then the f r a c t i o n a l change in intensity i n one wavelength i n t e r v a l would be (2) where Xn i s the s i z e of the fluctuation moving with a ccamcn velo c i t y , in units of 1 0 1 1 cm, r ^ i s the size of the star i n units of 1 0 1 2 cm, A> the s p e c t r a l resolution i n units of 0.3A, which was the spectrograph resolution used. A t y p i c a l v e l o c i t y gradient i s found by taking a terminal velocity of 1000 Km s \u00E2\u0080\u0094 1 reached over a distance of 10 1 2 cm. For a fluctuation which has a velocity gradient which i s the same as the wind, the minimum volume emitting in a given wavelength i n t e r v a l would be just the velocity s h e l l thickness cubed. In th i s case the intensity fluctuation i s A more r e a l i s t i c s i t u a t i o n might be i f a fluctuation of size X has cnly a thi n slab of thickness As moving at the appro-priate velocity to be i n the desired wavelength i n t e r v a l . In this case A I > - * * * * - 8 * - \u00C2\u00B0 - A2 rn-*. Htr A S (4) 20 The X Cep time series r e s t r i c t s the magnitude of an inten-s i t y f l u ctuation to less than 2% over the six hour span. Equation 4 then l i m i t s the size cf the largest region to change in t h i s time to 5x10 1 1 cm. These observations have confirmed the v a r i a b i l i t y of the s t e l l a r Hoc l i n e p r o f i l e over times longer than one day, and con-c l u s i v e l y show that the variation i s due to the change i n the p r o f i l e , net changing t e l l u r i c l i n e s . The amplitude of the i n -tensity change i n any one p i x e l i s only s l i g h t l y greater than what might be due to noise, but considering that groups of more than 10 pixels show the same change gives considerable c o n f i -dence to the physical r e a l i t y of the change. The one time series of > Cep has no convincing evidence for any short term var i a t i o n , or evolution of the p r o f i l e . The signal to noise i n the time series spectra i s only about 50, which was a constraint imposed on the maximum integration time by the detector c e d i n g system. 21 CHAPTER 3. SUPERSONIC ACCRETION Optical observations of v a r i a b i l i t y are averages ever the entire volume of emission at that p a r t i c u l a r wavelength. If the fluctuations i n the wind contain components cn a small scale compared tc the scale of the wind, the detection of flu c t u a t i o n s by way of technigues in which the integrated l i g h t i s observed, are limited by the si g n a l to noise which can be acquired. The discovery of two X-ray binaries imbedded i n s t e l l a r winds, namely Cen X-3 and 30170G-37 (=HD153919) allows the p o s s i b i l i t y of using the X-ray source as a probe of the s t e l l a r wind. Since the X-ray luminosity i s d i r e c t l y related to the rate of accre-tion of a small f r a c t i o n of the s t e l l a r wind onto the neutron st a r , the int e n s i t y of the X-ray source can be used with the aid of a s u f f i c i e n t l y detailed understanding of the accretion pro-cess to derive estimates of the density and velocity i n the wind. This was the subject of the published paper which has been attached as Appendix 1. A summary of the p r i n c i p l e r e s u l t s of the paper which support the conclusions of t h i s thesis i s given below. A schematic drawing cf the supersonic accretion process i s shown i n Fiqure 5 and the regions referred to are numbered i n the figure. The incoming gas, region 1, i s moving at a speed V with respect to the neutron star of mass M. The streamlines are bent in by the neutron star\u00C2\u00BBs g r a v i t a t i o n a l f i e l d . The mass and the velocity define the accretion radius, RA=2GM/VZ, which gives (apart from an e f f i c i e n c y factor which i s close to one) the cross section for accretion of material. The incoming gas str i k e s a shock cone t r a i l i n g the neutron star, c a l l e d the F i q . 6: Supersonic a c c r e t i o n schematic. 23 sheath, region 2. The gas i s shock heated to a temperature of 3x10* (r/10 l 4cm)-* K, If the density i s s u f f i c i e n t l y high the gas cools. In the sheath the gas loses i t s component of v e l o c i -ty away from the neutron s t a r , joins the accretion column and starts f a l l i n g down, region 3. Wear the accreting neutron star the X-ray luminosity may be large enough to raise the tempera-ture by Compton heating, The column w i l l then expand cut tc an almost spherical inflow, region 4, Eventually the flow encoun-ters the magnetosphere of the neutron star, region 5, below which the dynamics of the flow are regulated by the magnetic f i e l d . The gas s t r i k e s another shock a short distance above the surface of the neutron star, region 6, where the k i n e t i c energy of i n f a l l i s converted into thermal energy which i s mostly ra-diated away as X-rays. The table below gives length and time scales c h a r a c t e r i s t i c of the d i f f e r e n t regions. TABLE 2: SCALES IH SUPEBSONIC ACCBETION region size scale time scale star 1 0 1 2 cm 1 day accretion column 10*0-11 Cm 500 seconds magnetosphere 1 0 8 \u00E2\u0080\u0094 9 cm 1 second neutron star 106 cm 1 millisecond An analysis cf t h i s model yields several quantities which are d i r e c t l y related to the major parameters of interest i n the s t e l l a r wind, the s t e l l a r wind density, n , and v e l o c i t y . The luminosity of the unobscured source i s L =4.7x103* n 1,'v f t-3(fl/fl0)3(H l (/1O\u00C2\u00AB ci)-\u00C2\u00BB|S erg s~\u00C2\u00BB, (5) where p i s a factor usually of order one giving the e f f i c i e n c y 2 4 of the accretion, n n = n 6 / ( 1 0 1 1 cnr 3) , and Vg=V/108 cm s~ l. The angle that the shock cone makes with the axis of the accretion cclunn i s 0=2.7\u00C2\u00B0 2.2 n,,-*/ l s V Q S Z / I S . ( 8 ) The electron scattering o p t i c a l depth up the sheath i s less than that up the column i f n uV 8- 2<3.5, which i s independent of the estimate of the column temperature. With these simple relations in hand and some X-ray data of an object that i s c l e a r l y fueled by a s t e l l a r wind i t i s possible to confirm the model of the accretion process outlined above. , More importantly the observa-tions can be used to derive the density and ve l o c i t y in the wind. Two f a i r l y good sets of published data exi s t for the source Cen X-3, which appears to be the clearest case of accretion from a spherical supersonic wind. The source 3U1700-37 <=HD153919) would appear to be a very strong s t e l l a r wind source from i t s op t i c a l spectrum (see Fahlman, Carlberg, and Walker 1977) a l -though there are s i g n i f i c a n t e f f ects i n the spectrum associated with the period of the neutron star orbiting the 06f primary. be These e f f e c t s may^due to a wake of disturbed gas t r a i l i n g the neutrcn star (see Appendix 1). Or they may represent a s i g n i f i -25 cant d i s t o r t i o n of the s t e l l a r wind i t s e l f . In any case, the X-ray data from 3U1700-37 has a lower count rate than Cen X-3, hence greater s t a t i s t i c a l errors. In Cen X-3 the observational situ a t i o n i s almost the exact reverse; the X-ray source i s one of the brighter sources in the sky, but i t s o p t i c a l companion i s a 14th magnitude OB star which has been poorly studied (Conti 1S78) . The X-ray data for Cen X-3 shows several scales of v a r i a b i -l i t y ; a 4.8 second pulsation period, ascribed to the rotation of the neutron star and i t s magnetic f i e l d , a 2.1 day o r b i t a l period sometimes superposed with \"anomalous dips\", and an aperiodic change in the mean i n t e n s i t y l e v e l with a time scale of order one month. A p a r t i c u l a r l y exciting observation was made by Pounds, et a l . (1S75), who observed regular dips occur-ring every o r b i t during a t r a n s i t i o n from X-ray low to high state. Jackson (1975) proposed that the two d i s t i n c t dips were due to the reduction of the received flux by scattering i n the two sides of the sheath of the accretion column. He deduced a velocity of the wind with respect to the neutron star of between 375 and 620 km s - 1 , and a column semi-angle of 20\u00C2\u00B0. From egua-tion (6) and the v e l o c i t i e s guoted by Jackson, the implied column temperature i s i n the range 3,5-9. 6x10s K. Schreier j\u00C2\u00A7t a l . (1S76) estimate the density i n the wind as 1-5x10** cm - 3. These two estimates are consistent with the l i m i t i n g temperature of 1.5x106 K from Eguation (7), Accepting Jackson's proposal that the double dips are due to scattering in the sheath, but using the theory developed in Appendix 1, more informaticr can be derived from the observations. Pounds, et a l . (1S15) note 26 that the r e l a t i v e depths of the dips decrease as the source turns on. Also, from inspection of t h e i r published data one can see that the dips appear to become single as the source turns on. A schematic tracing of the X-ray intensity i s shewn i n Figure 7. From these observations and the model cf the long term variations proposed by Schreier, et a l , A rough trajectory of the variation of the s t e l l a r wind density and v e l o c i t y car be plotted, which i s shown in Figure 8. At point A the source i s in the X-ray lew state and at B the high state. Adopting the estimate of Schreier gt a l . for the low state density as 5x10* 1 cm-3, and high state density of 10*1 cm~3, f i x e s the densities at point A and B, but not the v e l o c i t y . T h e data shows that as the wind density decreases allowing the source to become v i -s i b l e , the velocity must be such that the o p t i c a l depth up the sheath exceeds the column o p t i c a l depth, and be close to one i n order to provide the deep dips. This puts point A near the Tc=1 l i n e . As the wind density drops the dips have a decreasing f r a c t i o n a l depth, which means that the velocity must be dropping fast enough that the density and velocity are moving further below the ,^= 1 l i n e . Eventually the dips become single as the density and v e l o c i t y cross the TT5 > f c l i n e . The combined den-s i t y and velocity v a r i a t i o n i s such that the accretion rate, and hence the i n t r i n s i c luminosity only increase s l i g h t l y while going from low to high state. The source s e t t l e s down at the high state, point B, with a density of 10 1 1 cm - 3 and a wind ve-l o c i t y with respect to the neutron star of about 500 Km s~ l. The c p t i c a l depth up the column i s so small that dips are not 28 F i g . 8: Density V e l o c i t y V a r i a t i o n of Cen X-3 2 9 seen regularly. When the source starts to turn off the data suggests that a d i f f e r e n t density velocity trajectory i s f o l -lowed, such that the o p t i c a l depth up the column and sheath always remains small. As the source approachs the point A again i t i s obscured by the increasing density of the s t e l l a r wind. If the trajectory of the variation of the wind ve l o c i t y and density i s schematically correct i t i s possible to draw a con-clusion as to the driving force f o r the wind. The trajectory i n Figure 8 suggests a c o r r e l a t i o n between the wind velocity and the wind density which would imply the acceleration of the gas up to the location of the neutron star increases as the density in the wind increases. In the r a d i a t i v e l y driven wind of CAK the acceleration of the wind varies as n - , where i s a con-stant s l i g h t l y less than one., This implies that the r a d i a t i o n acceleration should drop as the density increases. On the other hand Hearn (1975) suggests that the wind i s i n i t i a l l y acce-lerated i n a hot corona. The corona would heated by shock waves which grow from i n s t a b i l i t i e s within the atmosphere. As w i l l be discussed i n Chapter 5, one of the dominant i n s t a b i l i t i e s pre-sent i s the thermal i n s t a b i l i t y which grows cn a time scale which varies as n - 2 . This i n s t a b i l i t y may provide the c o r r e l a -tion between wind ve l o c i t y and density. The one month scale cf the high low state v a r i a b i l i t y is an enigua. There appears to be no natural scale i n the wind to explain i t , so i t may be connected to the subatmosphere of the star (Cannon and Thomas 1977 and Thomas 1973). The observations of Schreier et a l . (1976) contain e v i -dence that there are small scale (cf order 10 1 1 cm) fluctuations 30 in the wind. The count rate c l e a r l y varies with an amplitude greater that the s t a t i s t i c a l error on time scales cf about one hour. This time scale, which has been set by the spacecraft e^*^ or b i t and pointing mode, i s much longer than the natural res-ponse time of the accretion process, which i s about ten minutes, from Table 2. I t would be extremely interesting to have data with a time resolution of a few minutes to see i f the f l u c t u a -tions i n the wind become time resolved., In summary the theory of supersonic accretion that was de-veloped and applied to a limited amount of data on Cen X-3 shows that there i s a positive c o r r e l a t i o n between the wind density and vind v e l o c i t y during a period of t r a n s i t i o n from X-ray low tc high state..,\u00E2\u0080\u00A2 31 CHAPTER 4. PHYSICAL DESCRIPTION OF THE GAS The s t a b i l i t y of the Hind w i l l be investigated with a li n e a r i z e d s t a b i l i t y analysis. The analysis requires that the prevailing physical conditions be specified. This chapter i s devoted tc the derivation of the reguired quantities. The d i f -f i c u l t physical quantities are those describing the i n t e r a c t i o n of the gas and the s t e l l a r radiation f i e l d , which are the rate of energy gain and l e s s , and the radiation acceleration. Since t h i s interaction i s probably the key to the s t e l l a r wind, an accurate physical description must be used. In order to derive the cooling rate, heating rate, and r a -diation force i t i s necessary to know the d i s t r i b u t i o n cf atcms over the various stages cf i o n i z a t i o n , and the rate of absorp-tion and emission of rad i a t i o n by the ions. These calculations reguire atomic constants tc describe the radiation processes, which then become functions of the l o c a l density, temperature, and radiation f i e l d . The radiation f i e l d i s influenced by the flow of the gas, so that a good approximation to the radiation f i e l d would re-quire knowing the d e t a i l s of the flow. As an approximation I have taken the unattenuated, but geometrically diluted s t e l l a r radiation f i e l d . . This offers the advantage of retaining a com-ple t e l y l o c a l analysis at the cost of oversimplifying the radia-t i v e transfer. The approximation of the unattenuated f i e l d w i l l have the e f f e c t of somewhat over estimating the radiation force because overlapping l i n e s are ignored. As discussed l a t e r the e f f e c t i s l i k e l y to be at most a factor of two. The gas i s assumed t c be i n ionization equilibrium, which 32 i s v a l i d for time scales longer than the recombination time scale, 30 (T/10* K) 1 / 2 n n - 1 seconds. Equilibrium implies that the rate of t r a n s i t i o n s out of an ionization state i s balanced by the rate i n . For element i the rate out of i o n i z a t i o n state j i s determined by the rate of ionization to the next higher ion and the recombination rate to the next lower ion. The rate into the ionization state i s determined by recombinations from above and ionizations from below. Algeb r a i c a l l y , n c,j-i ( ne c j - 1 + \"5 hJ-i > * n 1 n e \u00C2\u00B0 S j# where in^,T) i s the c o l l i s i c n a l i o n ization rate out of j ^ CJ i s the photoionization rate rt\u00C2\u00A3j{ng,T) i s the recombination rate from l e v e l j to j-1. These io n i z a t i o n balance equations were solved f o r as many atoms cf s i g n i f i c a n t s t e l l a r abundance for which gccd atomic data was available. The elements used are shown with their as-sumed abundances i n the accompanying table. , It would have been desirable to have included Nickel and Iron with their f a i r l y high cosmic abundance and great number of spectral l i n e s , but no r e l i a b l e and consistent set of data for a wide temperature range could be found. TABLE 3: ATOMIC ABUNDANCES ELEMENT Z ABUNDANCE Hydrogen 1 1.0 Helium 2 8.5x10-2 Carbon 6 3.3x10-* Nitrogen 7 9.1x10-5 3 3 Oxygen 8 6.6x10-* Neon 10 8. 3x 10-s Magnesium 12 2.6x10-s S i l i c o n 14 3.3x10-s Sulfur 16 1.6x10-5 These abundances were taken from Alien (1973). Standard rates were used for a l l the photoionization cross sections, recombination rates, and c o l l i s i o n a l ionzation rates. But since the gas has a f a i r l y high density (order 10*1 c a - 3 ) and i s in an intense radiation f i e l d i t i s necessary to make some corrections. The density e f f e c t s are allowed for by adding corrections to the recombination rate for three body recombina-tion s , and recombination to upper levels. A small correction for i o n i z a t i o n out of upper l e v e l s i s also included. The great-est d i f f i c u l t y i s allowing for the e f f e c t of both the radiation f i e l d and the density e f f e c t s on the di e l e c t r o n i c recombination rate. This process depends upcn captures to levels of large guantum number, and i t i s possible that these l e v e l s may be reionized before they can s t a b i l i z e by cascading down to lcwer l e v e l s . These e f f e c t s have been crudely allowed for by c a l c u l a -ting a m u l t i p l i c a t i v e correction factor, based on a f i t to the guantum mechanical ca l c u l a t i o n s cf Summers (1574), A l l these rates and corrections are discussed i n Appendix 2. The Ionization Balance The solution to the io n i z a t i o n balance equations i s very simple since the lowest l e v e l only interacts with the second l e v e l , and then the second l e v e l i s linked to the f i r s t and 34 t h i r d , and sc on., This gives the r a t i o of the population in a ion i z a t i o n state to the population i n the next lower l e v e l . A normalization completes the solution. , The equations are weakly nonlinear through t h e i r dependence on the electron density, but usually two or three i t e r a t i o n s s u f f i c e s for an accuracy of about 1 part in 10 6. The results are given i n terms of the ion i z a t i o n f r a c t i o n X ij for ion j of atom i , where xt*y summed over j i s unity. To get the number of atoms of type i , j we take the product X t j-A^n, where At' i s the abundance of atom 1. In Figure 9 the ionization balance for a gas of density 10 1 1 cm - 3, in the undiluted radiation f i e l d of the star, i s shown for a range of temperatures. It i s found that for the range of densi-t i e s of interest the reduction of the d i e l e c t r o n i c recombination rate by the density and radiation f i e l d e f f e c t s i s s i g n i f i c a n t and tends to s h i f t the ionization s l i g h t l y to higher stages of ioni z a t i o n . At very high densities the d i s t r i b u t i o n approaches to the .distribution expected f o r ITE. The heating and cooling r a t e for a gas of density 10-*1 cm - 3 in a undiluted radiation f i e l d are shown in Figure 10. , The plotted quantities are the c e d i n g and heating rates, - A and T , respectively. The plotted guantities are to be multiplied by the density squared to obtain the rates per cm - 3. The genera-l i z e d cooling rate i s taken as Jt* =n2 { A - r ). The quantity *\u00E2\u0082\u00AC/n2 i s plotted. The radiative equilibrium between the photoionization cea-tinq and radiative losses holds at temperatures of about 2x10* K for densities around 10*4 cm - 3. ; This i s shown i n Fiqure 10 for zero velocity of the qas with respect to the star. Of intere s t Sulfur S i l i c o n Magnesium He on Oxygen Nitrogen Carbon at Helium , \u00E2\u0080\u0094 i \u00E2\u0080\u0094 , \u00E2\u0080\u0094 \u00E2\u0080\u0094 ~ - j \u00E2\u0080\u0094 , .\u00E2\u0080\u0094\u00E2\u0080\u00A2 \u00C2\u00BB 1 5 6 7 Log T F i q . 9 : I o n i z a t i o n Balance 3 6 o Fig. 10: Heating and C e d i n g Bates for Solar Abundances 37 to the \"warm radiation acceleration\" model i s that near 2x10s K the loss rate i n an o p t i c a l l y thin medium i s at a maximum. Such a temperature would be very d i f f i c u l t to maintain in the gas, reguiring an immense input of energy from seme other heat source. Between 10* and 10 7 K the lo s s rate drops to a minimum where the ra d i a t i v e losses would be more e a s i l y balanced. The radiation losses from such a hot gas would consist l a r g e l y of X-rays, which would be suitable for producing the 0 VI ion, as has been suggested by C a s s i n e l l i and Olson (1978). The gas i s thermally unstable (see Fi e l d 1965) to both i s o -choric and i s o b a r i c disturbances when the temperature derivative of the generalized cooling rate at constant density i s negative. The temperature gradient i s not quite steep enough (logarithmic derivative of the cooling rate less than about -3) at any pcint to admit isentropic i n s t a b i l i t y , wherein ordinary sound waves gain energy i n the rarefactions and lose i t i n compressions. Even i f there i s a s l i g h t inaccuracy in the calc u l a t i o n s such that t h i s isentropic i n s t a b i l i t y condition could be met, i t would appear i n a very narrowly defined temperature i n t e r v a l . Calculations by Raymond et a l . (1S78) indicate that with the inclus i o n of the iron group elements the slope becomes even less steep, and the gas i s further away from isentropic i n s t a b i l i t y . The loss rate and i t s derivative turns out to be c r i t i c a l to the s t a b i l i t y cf an accelerating atmosphere, so i t has been plotted i t for the CNQ elements enhanced by a factor of 10 i n Figure 11. Obviously the abundance has a strong e f f e c t on the c o d i n g rate, since the CNO elements are responsible for the c o d i n g i n the range 10 s to 10* K. 39 The s t e l l a r wind i s usually o p t i c a l l y thin at o p t i c a l and longer wavelengths for continuum emission, but can become o p t i -c a l l y thick i n the resonance l i n e s , which provide the l i n e coo-l i n g as well as most of the radiation acceleration. Using the alt e r a t i o n to the l o s s rate of Hybicki and Hummer (1978) the reduced cooling rate i s shown i n Figure 12 for a velocity gra-dient of dv/dz=10 - 3. Note that the s p e c i f i c cooling rate (units of erg cm - 3 s~ l) w i l l s t i l l increase approximately l i n e a r l y with density, since the losses vary with the cooling rate i n erg cm+3 s - 1 times the density squared, over the o p t i c a l depth., This i s a very rough c a l c u l a t i o n , since no allowance has been made for the change of the l o c a l i n t e n s i t y due to the o p t i c a l l y thick l i n e s . Radiation Force The radiation force i s defined as ij C J (9) where ffi = J A-iiii , and m t- i s the atomic weight of the various ions. I f the unattenuated radiation f i e l d i s used i t provides an upper l i m i t to the radiation force. A more r e a l i s t i c e s t i -mate i s supplied by the method used by Castor, Abbott, and Klein (1975) , which i s based on an analysis of the radiative transfer in one spectral l i n e o r i g i n a l l y done by Lucy (1971). With the aid cf the Sobolev approximation the problem can be solved and i t i s found that the force due to li n e s i s ^ r a d - ^ rad \u00E2\u0080\u0094\u00E2\u0080\u0094-* (10) F i g . 12: C o o l i n g with o p t i c a l l y t h i c k l i n e s . 41 where T = Tte 2/(mc) f ; J (1) hL X ; J nc[ 1/2 ( I * / * 2 ) {dv/dz-v/r) +v/r ]-\u00C2\u00BB and Tfe2/(mc)=. 02654 9 \"rad i s acceleration in an o p t i c a l l y thin gas f/j(/) i s the o s c i l l a t o r strength for l i n e 1 of atom i i o n i z a t i o n state j , c i s the speed of l i g h t JJL i s cosine of the angle subtended by the s t e l l a r radius from the point i n the gas. In addition to the force on the l i n e s there i s the force on the electrons, - f t * * rr* ^c-q e - \u00E2\u0080\u0094 \u00E2\u0080\u0094 e \u00E2\u0080\u0094 J C Tl TYV <12) where F i s the fl u x integrated over a l l freguencies, and i s the Thomson cross section. The force on the electrons i n the undiluted radiation f i e l d i n a ecu t i e t e l y ionized gas i s 194.7 cm s - 2 . There also i s the force on the continuum, which i s usually guite small, with the undiluted radiation f i e l d at a density of 10** cm - 3 i t i s 63.48 cm s - 2 . The l i n e acceleration i s dominated by o p t i c a l l y thick l i n e s , and increases almost l i n e a r l y with the velocity gradient. A schematic of the acceleration as a function of dv/dz i s shown in Figure 11 below. In Figure 13, the acceleration i s a weak function of temperature i n the range 10* \u00C2\u00ABd= 9e 8 f t ) , <13> where t-cr er, ev t K (dv/dz)-*, ge i s the radiation force on the elec-trons, and v i s the thermal v e l o c i t y . I f i n d that K(t) = .067 t - Q - 9 * f o r n=10io, and M(t) = .022 t~\u00C2\u00B0 \u00E2\u0080\u00A2** for n=10 1 3 cm - 3, whereas CAK find .033 tr0,7, which i s good agreement. There are two reasons for the density dependence of the acce-l e r a t i o n . F i r s t , a few of the lines go from o p t i c a l l y thick to thin as the density gees down, and secondly, the ionization bal-ance i s density dependent i n t h i s c a l c u l a t i o n , through the a l l o -wance for c c l l i s i c n a l i o n i z a t i o n , and through the density depen-dence of the rate c o e f f i c i e n t s . The d e f i c i e n c i e s i n th i s c a l c u l a t i o n of the radiation force are due to a somewhat li m i t e d l i n e l i s t , mostly due to the lack of any iron group elements, and mere seriously a very simple treatment of radiation transfer. Within the approximation used these two d e f i c i e n c i e s cancel each other out to a certain ex-tent. The radiation force has been over estimated by not taking account cf overlapping l i n e s , which would involve formulating a model cf the atmosphere intervening between the point i n the gas and the st a r . The radiation force increases with the number of l i n e s present, but the flux available decreases as the number of lin e s goes up., Klein and Castor (1978) have reported on new 44 calculations made by Abbott of the radiation force. He finds that the o r i g i n a l CAK law i s bracketted by twc alternative tran-sfer schemes, and probably the CAK law represents a good approx-imation to the force. The calculations here are i n good agree-ment with the CAK law. The l i n e acceleration varies approximately as (n tJ/n) ((dv/dz) /n,-j )* where i s in the range 0.7 to 0.9. As the velo-c i t y gradient increases a l l l i n e s become o p t i c a l l y t h i n , and the force l e v e l s off at the maximum value. This means that the force depends on the abundances roughly to a power in the range of .1 to .3, which i s a very weak function. Therefore, the r a -diation force i s i n s e n s i t i v e to the assumed abundances for flows in r a d i a t i v e equilibrium because most of the l i n e s are o p t i c a l l y thick. One aspect of the radiation transfer which i s important to the analysis of the s t a b i l i t y of the flow i s the shape of the l i n e s , which can provide an immediate source of i n s t a b i l i t y , as has been reported by Nelson and Hearn (1978). The i n s t a b i l i t y they f i n d only acts i n subsonic flow. This has been l e f t out because i t i s dependent on the d e t a i l s of the radiation trans-fer. H2fSUium Balance The number of free zero order quantities can be reduced by reguirinq that the eguations of mass and momentum conservation be s a t i s f i e d . In the case of r a d i a t i v e equilibrium the tempera-ture i s determined by the balance of heatinq and coclinq, other-wise the temperature i s just a r b i t r a r i l y s p e c i f i e d . , 45 lo i l l u s t r a t e the solutions cf the mass and momentum equa-tions the one dimensional equation of mass conservation i s subs-ti t u t e d into the momentum balance equation (see CAK and i n Chapter 5, below) with zero temperature derivatives. where the form of the mass conservation equation for a spheri-c a l l y symmetric system has been used. Spherical qecmetry has been used partly because the density qradient remains negative even i f the velocity gradient acquires a small neqative value. The perturbations are in the form of plane naves, so a l l deriva-tives w i l l be made with respect to the heiqht z, instead of r. The independent variable i s chosen to be dv/dz. In Figure 14 the two sides of the momentum equation are shown as functions of dv/dz. For supersonic flow, v2>2RT, there i s always one dec-celerating solution, where the radiation force i s l e s s than the qra v i t a t i o n a l f i e l d . As can be seen from Fiqure 11, i f the ve-l o c i t y i s n * t too larqe there are two solutions in which the qas i s accelerated outwards. For t y p i c a l s t e l l a r wind conditions the two solutions have dv/dz approximately equal to 10-* and 1, The two solutions are acceptable l o c a l l y , but boundary and con-t i n u i t y conditions may rule out the hiqh qradient solution. By imposinq continuity of vel o c i t y from subsonic to supersonic flow CAK r e s t r i c t themselves tc the low qradient s o l u t i o n . The solu-tion with the larqe v e l o c i t y qradient i s acceleratinq so rapidly that the wind becomes o p t i c a l l y t h i n i n the resonance l i n e s . This means that i f the acceleration could be maintained over a distance of 0. 151 of the s t e l l a r radius, the qas would be acvinq a c c e l e r a t i o n g + g ,, low density - \u00E2\u0080\u0094 -rad g + 8 r a d > h i \u00C2\u00A7 h d e n s i t y dv/dz < 0 dv/dz > 0 dv/dz > 0, high gradient s o l u t i o n 2RT 2RT Fig. 14; The Momentum Equation Solution 47 at the terminal v e l o c i t y , although t h i s solution i s physically acceptable, observational evidence suggests that i t may not be realized. As can be seen from Figure 14, i f the v e l o c i t y becomes too large no accelerating solution can be found. The maximum velo-c i t y f c r which accelerating solutions exist varies with the gra-vity and the density of the gas, A table i s given below which indicates the maximum v e l o c i t y giving an outward acceleration. In Table 4 the Vmax column gives the maximum velocity at which a positive dv/dz can be found, the value of which i s given i n the next column. Two values for the gravity are used to show that the maximum velocity with dv/dz>0 i s mostly effected by the gas density. The gas was chosen to be in r a d i a t i v e equilibrium, which gives a temperature of 2x10* K. , TABLE 4: LIMITING VELOCITY FOR ACCELERATING SOLUTIONS g=10*cm s- 2 g=4000cm s~2 Vmax dv/dz Vmax dv/dz n=10\u00C2\u00BBo 4.1x108 .16x10-3 4.45x10\u00C2\u00AB .84x10-* n=10\u00C2\u00BB> 4.7x107 .16x10-2 5.2x10* .16x10-2 n=10\u00C2\u00BB2 5.7x10* .15x10-2 6.0x10* .64x10-2 Table 4 shows that the maximum velocity i s approximately inversely proportional to the density. The maximum velocity f o r acceleration decreases nearly to the sound speed at a density of I O 1 2 cm - 3. Below t h i s velocity the Sobolev approximation used for deriving the radiation acceleration i s i n v a l i d . I f a large portion of the flow, thicker than one Scbclev s h e l l , (the sound speed divided by the velocity gradient) ac-48 q u i r e s a v e l o c i t y which i s g r e a t e r than the maximum f o r a p o s i -t i v e v e l o c i t y g r a d i e n t i n momentum b a l a n c e , then the gas w i l l d e c c e l e r a t e . T h i s s i t u a t i o n c o u l d a r i s e i f the f l o w i s a chao-t i c medium i n which elements o f t h e f l u i d a r e p r o p e l l e d t o v e l o -c i t i e s i n excess c f t h e maximum f o r a c c e l e r a t i o n , o r have a den-s i t y i n c r e a s e which makes t h e v e l o c i t y g r e a t e r t h a n the maxiaum. The wind might c o n s i s t o f many, q u i t e l a r q e p a t c h e s , which a r e b e i n q a c c e l e r a t e d and d e c c e l e r a t e d w i t h r e s p e c t t o one a n o t h e r . Hhere these r e q i o n s c o l l i d e t h e i r s u p e r s o n i c v e l o c i t i e s v c u l d ensure shock h e a t i n q which would produce t e m p e r a t u r e s a p p r o p r i -a t e t o 0 VI and l i k e i o n s . T h i s shock heated qas would c n l y comprise a s m a l l p o r t i o n c f the t o t a l qas i n the f l o w , and a f t e r f o r m i n q would be blown away from t h e s t a r . , 49 CHAPTER 5. THE STAEILITY ANALYSIS The observations suggest that a s t e l l a r wind i s an extreme-ly variable, inhomogenous flow. On scales of a day to years there are large general variations, which may originate within the star. The X-ray observations suggest small scale f l u c t u a -tions cf order 10 1 1 cm. This observed v a r i a b i l i t y could have two causes: the flow may s t a r t out i n the lower atmosphere as smooth, and then enter a region of i n s t a b i l i t y where i t breaks up; or the existence of the flow may be depend upon some i n s t a -b i l i t y . In t h i s section the l o c a l s t a b i l i t y of the flow w i l l be investigated. This w i l l be done by considering the propogation of i n f i n i t e l y small disturbances, i . e . a l i n e a r i z e d analysis, with wavelengths short compared to the scale of variation within the s t e l l a r wind. This analysis i s directed towards finding i n s t a b i l i t i e s that are rapid amplifiers, i . e . the growth time scale i s shorter than the time to move one scale length i n the atmosphere; and absolute i n s t a b i l i t i e s , which can actually gen-erate o s c i l l a t i o n s cr lead to \"clumps\" within the wind, . One major l i m i t a t i o n of t h i s analysis i s that i t has only bee* done for cne dimensional wave propogation, that i s the waves can only have a velocity component which i s oriented along the d i r e c t i o n of propagation. Eor instance, t h i s immediately rules out the p o s s i b i l i t y cf the Rayleiqh Taylor instability,, (Krolik 1577, Nelson and Hearn 1978). Similar analyses, but with more approx-imations in the l i n e a r i z a t i o n , have been performed by Hearn (1972) and for guasars by Hestel et a l . , (1976). The basic equations that apply are the conservation cf mass 50 115) where -f i s the mass density and v i s the gas v e l o c i t y . The con-servation of momentum neglecting the v i s c o s i t y i s given by, T (16) where 9(which w i l l be sometimes abbreviated as g r) i s the acceleration due to r a d i a t i o n , P i s the gas pressure, and g i s the g r a v i t a t i o n a l acceleration. The conservation of energy i s expressed, 0* (17) where e and h are are respectively the s p e c i f i c i n t e r n a l energy and enthalpy. , ~\u00C2\u00A3 i s the generalized cooling rate i n the frame of the gas, defined as / = L-(1-v/c) G, where v i s the velocity of the gas r e l a t i v e to the star and L and G are the l o c a l s p e c i f i c c o d i n g and heating rates, respectively. , The thermodynamic re-lation s reguired are an eguation cf state P = kT(n+.ne) \u00C2\u00AB where n e = (j-1)n t*;. (18) The sum i , j i s over the ion i z a t i o n states and the elements, re-spectively. The number density of atoms and electrons are n and n The i n t e r n a l energy i s e = 3/2 kT(n+n ) + Zn \u00C2\u00A3; f:.t , (19) where \" X ; . - i s the i o n i z a t i o n energy of ion i , j with density n ty. J J The enthalpy i s defined as, h = e *P/^ . (20) For the conductivity \u00E2\u0080\u00A2% the standard value cf Spitzer (1962) has 51 been used. The above equations are l i n e a r i z e d i n order to obtain a dispersion r e l a t i o n , which i s a polynomial describing the prcpar gaticn of waves of i n f i n i t e s m a l amplitude. . The l i n e a r i z e d egua-tions are obtained by imposing a perturbation on the tempera-ture, density and ve l o c i t y cf the form Q(z,t) = Q D(z,t)+\u00C2\u00A3gt \u00C2\u00A9\u00C2\u00AB\u00C2\u00BB in an i n f i n i t e atmosphere. The c r i t e r i o n means that i n the neighbour-hood of the solution (oe,ko) to the two equations the root varies as u) = u0 + A ( k - k e ) 2 , where a i s a constant. This imp-l i e s that an absolute i n s t a b i l i t y i s a saddle point of the ima-qinary part of the frequency as a function of k. The imaqinary part cf the frequency w i l l be at a maximum with respect to rea l k at the solution, and t h i s freguency w i l l dominate the qicwth rate. These two nonlinear equations, D=0 and dD/dk=0, were solved simultaneously with the aid of a computer routine, using the l o c a l maximum of the imaqinary part of the freguency for re a l wavenumbers as a startinq point. An attempt was made to find common roots to the two equations by ccnstructinq the dis-criminant of the c o e f f i c i e n t s of the two equations. This was unsuccessful because of the im p o s s i b i l i t y of retaining s u f f i -cient numerical accuracy. In order to understand the dispersion r e l a t i o n and the physical o r i g i n of the roots, analytic expressions for the roots w i l l be derived for a number of simple l i m i t i n g cases. The roots i n a complex s i t u a t i o n can be understood as superposition of these several simple cases. These l i m i t i n g solutions have been derived with the aid cf numerical solutions, and unless noted the calculated roots plotted came from a dispersion r e l a -t i c n with c o e f f i c i e n t s calculated from a gas in an undiluted radiation f i e l d , with a density of 10 1 1 cm - 3, and a velocity of 100 Km s\u00E2\u0080\u0094*. The res u l t i n g equilibrium guantities are i n cgs u n i t s : T = 1.87x10* K; =0, d /'/dT = .46x10-*, d^/dn = 2.x10-i 3, dv/dz = .2x10~ 3, dn/dz = -2.1, g rad = 1.18x10*, dgr/dT = -0.3x10 -1 , and dgr/dn = -.15x10-*. I t i s found that the char-acter of the roots changes l i t t l e with a variation of the physi-c a l parameters around these values for t y p i c a l s t e l l a r wind con-dit i o n s . Case 1: Sound Haves In An Atmosphere The simplest case which has a non zero growth rate i s a wave propogating v e r t i c a l l y i n a s t a t i c , isothermal atmosphere, with no conduction or radiat i o n present. In t h i s case the d i s -persion r e l a t i o n as given i n the Appendix 3 reduces to Defining H=n/(dn/dz) , t h i s has solutions o>=0 and in the l i m i t of large and small k the n o n t r i v i a l roots become l a . - 5 * 00 to \u00E2\u0080\u0094\u00C2\u00BB C v (22) This i s e s s e n t i a l l y the well known solution of Lamb (1945) to the problem of wave propogation in an isothermal, exponential atmosphere. But note that the value of H, the scale height of F i g . 15: Pseudc Isothermal S t a t i c atmosphere Boots 56 the density gradient, used in the numerical calculations was not the isothermal scale height, but that the scale height was de-termined by the v e l o c i t y gradient through the mass conservation equation. In the short wavelength l i m i t (k-\u00C2\u00BBa>) the waves move at a phase and qroup v e l o c i t y equal to the ordinary sound velo-c i t y . Outward movinq waves are amplified and inward moving waves are damped at a rate such that the momentum carr i e d in the wave i s kept constant. These waves are not absolute i n s t a b i l i -t i e s . At lonq wavelengths (k\u00E2\u0080\u0094>0) the r e a l part cf the freguency goes to a f i n i t e l i m i t , c a l l e d the acoustic cutoff frequency, and the damping goes to zero. . This means that these waves have a phase velocity going to i n f i n i t y , but the group v e l o c i t y goes to zero and no energy i s prcpcgated. Physically t h i s cutoff results from the atmosphere as a whole moving with the wave mo-t i o n , rather than a wave propogating away from the source. The change over between the two l i m i t i n g solutions occurs for k of order H _ 1. The solution i s i l l u s t r a t e d in the accompanying Figure 15. In Figure 15, and a l l other graphs of the roots of the dispersion r e l a t i o n , the logarithm (base 10) of the r e a l and imaginary parts of the wave freguency are separately plotted against the logarithm of the wave number. , On the graph of the real part a symbol ^ o r O) on the l i n e means that the r e a l part i s negative. On the graph of the imaginary part the same symbols indicates that the wave i s unstable at that wave number, that i s , the imaginary part i s positive. Note that f r e -quently tfee two acoustic roots have an i d e n t i c a l magnitude, but opposite sign, so that in the plot the two l i n e s l i e on top of each other. 57 The plots are done for k ranging from 10- 1 5 to 1 cm - 1, which i s an u n r e a l i s t i c a l l y large range for the physical s i t u a -t i o n , but i s done to i l l u s t r a t e the asymptotic l i m i t s of the roots. The physically acceptable range of wave numbers i s for wave numbers les s than the a wavelength of a s t e l l a r radius, 10-n cm\u00E2\u0080\u0094i, to a wavecumber corresponding to one mean free path, about 10 - 2 cm\u00E2\u0080\u0094i. There i s a maximum freguency for which the solutions are v a l i d , set by the longer time scale, recombination or the elec-tron ion thermal equilibrium. The recombination freguency i s o>t&c = .188 n l k (T/10* K ) - i / 2 s - i , (23) and the electron ion equilibrium frequency i s UeL = 7x10* n \u00E2\u0080\u009E \u00E2\u0080\u00A2 (T/10* K ) ~ 3 / 2 s-\u00C2\u00BB... <2<4) The maximum frequency for which the calculations are v a l i d then i s the minimum of (Jr<.c and tJei. . The minimum freguency of i n -terest would be determined by the time for the complete replace-ment of the sta r ' s s t e l l a r wind envelope. This freguency i s about 6x10-5 s - i . Case 2: The Effect Of Conduction Allowing conduction a f f e c t s mostly the short wavelength roots. Taking the dominant terms i n the dispersion r e l a t i o n gives. 58 For t h i s case the dominant terms cf the roots f o r k are, f t v d (25) which i s a heavily damped non propagating disturbance. The sound waves are given as * * 1 / (26) which are isothermal sound waves, and always damped independent of d i r e c t i o n cf propagation. The numerical solution shows that the analytic solutions only apply for k>10 - 3, and that the slow root has a small r e a l part at short wavelengths. Case 3: Radiation E f f e c t s In the long wavelength l i m i t we expect radiation e f f e c t s to be dominant. The dominant terms of the dispersion r e l a t i o n becc me w * { . , f c . R ] t u>\u00C2\u00BB + u {c2 %4\u00C2\u00A3] + 0 \u00E2\u0080\u00A2 c * R. ct T * <29) which e s s e n t i a l l y i s the thermal s t a b i l i t y condition. A parcel of qas with d J^/dT <0 would probably tend to collapse. In a general case i f d J^dT were negative, part of the gas may cool and ccllapse, and other parts may rise i n temperature. The ex-istence of the hot, low density component depends on an appro-priate heat source to maintain a temperature of order 10* to 10 ? K, where the gas i s stable. If t h i s bistable mode i s possible within a s t e l l a r wind, i t may lead to a two component outflow with a cool (T around 2x10* K) and hot (T around 10 7 K) com-ponent. The hot component may be able to supply s u f f i c i e n t C VI atoms that there wculd be no need for a coronal region. The dominant terms cf the other two roots are sound waves, civ , l-Z R T + h rl-r - J U = ' ' ' PI ( 3 0 ) \"The roots are shown i n Figure 16, for a gas with a nonzero velo-c i t y and acceleration. From Equation 30 we make the discovery that deccelerating flows are unstable, and the numerical calcu-l a t i o n finds that i t i s an absolute i n s t a b i l i t y . An example of this i n s t a b i l i t y i s discussed l a t e r and i l l u s t r a t e d in Figure 20. Defining some basic time scales as t(dynamic) = V j d v / d z j , t (cool) = f c> R. / \ c| c|T I , P i g . 16: I s o t r o p i c B a d i a t i o n F i e l d 61 t (acoustic) = H/Cp where cx= Jk2BTi The condition for the i n s t a b i l i t y of deccelerating flews (Eg. 28) i s , t (dynamic) x t(cool) << (t (acoustic) ) 2 Note that the conductive damping dominates the roots for k > 10-*. Allowing a radiation acceleration, gives the roots as plotted in Figure 17. The asymptotic l i m i t s are not changed by the radiation acceleration, but the inward propagating acoustic wave i s unstable in the range cf wavenumber 10 _ 1 1< k < 10~ 7. The resul t i n g growth rate i s close to 1200 seconds, but the i n s t a b i -l i t y i s only amplifying. The cooling due tc c o l l i s i c n a l l y excited l i n e s may be d i -minished when the gas becomes o p t i c a l l y thick in the resonance l i n e s . The e f f e c t of t h i s has been approximated by turning the loss rate o f f , but leaving the heating on. The roots of the dispersicn r e l a t i o n i n th i s case are shown i n Figure 18. Be-sides the amplifying i n s t a b i l i t y from the radiation force there i s an additional range of i n s t a b i l i t y for both inward and cut-ward acoustic waves for 10~7< k <10 - s. This behaviour r e s u l t s from the term d jf/dn becoming si g n f i c a n t . Figure 19 shows the e f f e c t of a thermal i n s t a b i l i t y , d^/dt<0. The \"slow\" root has a rapid growth rate, which i s an absolute i n s t a b i l i t y . The acoustic roots are changed only s l i g h t l y , the amplification acting over a narrower range of wa-venumber, and not guite as r a p i d l y . Figure 20 shows the pressure dominated thermal i n s t a b i l i t y which i s present at 10 7 K. In t h i s case the flow speed i s sub-L O G- K F i g . 17: R a d i a t i o n Force Cn 6<4 sonic, and the radiation force i s less than 5% cf gravity. The i n s t a b i l i t y r e s u l t s when Eguation 28 i s v i o l a t e d . The acoustic waves have an absolute i n s t a b i l i t y at long wavelengths. The growth rate i s proportional to the element abundance through the cooling rate. An atmosphere of density 1 0 1 2 cm - 3 and a velocity of 100 Km s~ 1 exceeds the maximum ve l o c i t y for an accelerating solution. The roots when dv/dz i s negative are shown in Figure 21. Shis i s an absolute i n s t a b i l i t y at long wavelengths. In summary the roots of the dispersion r e l a t i o n can be un-derstood i n terms cf combinations of the roots which occur i n simple physical cases. For k> 10~* cm - 1, the conduction always provides a strong damping, es p e c i a l l y t c the slow wave. Thus i t can be concluded that i n the long wavelength l i u i t , k< 1 0 - 1 1 , the behaviour depends on the shortest time set by: the acoustic time scale, egual to the scale height divided by the sound speed; the dynamic time scale (dv/dz) - 1; $s. ? / \u00C2\u00AB ( T | . There i s always a \"thermal wave\", that i s , a slow moving wave, compared to the scund speed, with a growth rate set by the thermal time scale. The slow wave i s an absolute i n s t a b i l i t y i f the derivative d f/dT i s negative. I f t(acoustic) \u00C2\u00BB t (dynamic) then the acoustic waves have growth rates given by the dynamic time scale. These acoustic waves w i l l be absolutely unstable i f the velocity gradient i s nega-t i v e , A hot gas, with t(acoustic) < t (dynamic) w i l l have an absolute i n s t a b i l i t y a r i s i n g from the acoustic waves. At T=107 the growth rate i s about one hour, (one hour at 3x10* K where d//dt<0) which increases as n 2, u n t i l t (acoustic) exceeds L O G * F i g . 20: High Temperature I n s t a b i l i t y Fig. 21: Decceleraticn I n s t a b i l i t y at n=10 1 2 c u r 3 68 TABLE 5: THE DISPEBSICfr BELATION \u00C2\u00A3 PLOTTED Fig. n T V dv/dz d /dT g Bemarks 15 10*i 2x10* 0 0 0 0 pseudo isothermal 16 10*i 2x10* 100 2x10-* 5x10-s 0 i s o t r o p i c 17 10i\u00C2\u00BB 2x10* 100 2x10-* 5x10-s 1x10* standard case 18 10i* 2x10* 100 2x10-* 5x10-s 1x10* no cooling 19 10*\u00C2\u00BB 5x10s 100 7x10-* -3x10-* 1x10* thermal 20 10** 1x107 100 5x10-5 8x10-8 270 high temperature 21 10*2 2x10* 100 -6x10-* -3x10 -3 5x103 decceleration t(dynamic). The radiation acceleration only acts to provide an amplifying i n s t a b i l i t y for inward acoustic waves. This a m p l i f i -cation acts for an observaticnally interesting range of wave-lengths, from 5x10 7 to 5x10** cm, with growth times cf order 1200 seconds. on the basis of t h i s analysis the o r i g i n a l CAK \"cool\" atmo-sphere i s stable only i f no waves are sent into the accelerating wind from lower layers. Otherwise the radiation force acts tc provide an amplification of the inward moving (with respect to the gas, but outwards with respect to the star) acoustic wave. A corona with a temperature of several m i l l i o n degrees w i l l always have an absolute i n s t a b i l i t y , either the ordinary thermal i n s t a b i l i t y of the radiation losses, or the \"high temperature\" i n s t a b i l i t y outlined above, which has a growth time cf order an hour. The semi empirical model of the wind proposed by C a s s i t e l l i et a l . (1978) has the wind heated with an i n i t i a l acceleration in a hot corona, and then cooling i n the outer r a d i a t i v e l y acce-lerated 2one. The s t a b i l i t y analysis leaves no doubt that t h i s s i t u a t i o n would be expected to show f l u c t u a t i o n s . The hot corona i s subject to i n s t a b i l i t i e s which may be responsible for creating the shock waves to heat i t . Bemnants of these f l u c t u a -tions would be carried out into the wind where the length scales of 10 7 to 10** cm would be amplified. 69 CHAPTER 6. , CONCLUSIONS The program of o p t i c a l observations conclusively shows that s t e l l a r winds do vary on time scales of one day and more. , These observations were taken at s u f f i c i e n t l y high resolution that any variations cf the t e l l u r i c l i n e s could be separated from varia -tions of the s t e l l a r l i n e s . A star which has often been re-ported as varying. Lambda Cephei, was monitored with a time re-solution of one hour over a period of six hours but no s i g n i f i -cant variation was seen i n the H ^ l i n e . This n u l l observation puts an upper l i m i t of about 5x10 1 1 cm on the size of any \"blobs\" i n the wind. Day to day v a r i a b i l i t y was confirmed in X Cep and \u00E2\u0080\u00A2x'Cam, but not conclusively for & Ori. These variations may not be due to fluctuations within the wind i t s e l f since t h i s time scale i s long enough to allow complete replacement of the material i n the l i n e formation region. Causes of the long time scale v a r i a t i o n include rotation of the s t a r , i n t e r n a l o s c i l l a -tions of the star, or a variation of the the emergent flux and hence the driving force for the wind. The analysis of the X-ray observations cf Cen X-3 provides confirming evidence f o r the suggested mechanism causing the long term X-ray i n t e n s i t y variation reported by Schreier et a l . (1976)\u00E2\u0080\u00A2 That i s , the wind density varies s u f f i c i e n t l y that the source i s occasionally smothered by the opacity of the s t e l l a r wind. In addition I have found that as the density i n the wind changes, i t must be correlated with the wind velocity i n order to explain the changing character of the anomalous dips i n the i n t e n s i t y at non-eclipse phases. Besides these semi-regular dips the X-ray in t e n s i t y shows i n t e n s i t y fluctuations on a time 70 scalei. of less than one hour, which i s probably due to changes in the amount of mass being accreted by the neutron star, The natural source f o r the v a r i a t i o n in the accretion rate i s the variation of the density and velocity in the s t e l l a r wind with size scales of 1Q*o tc 10*1 cm. The t h e o r e t i c a l analysis of the s t a b i l i t y of a wind finds a number of sources cf i n s t a b i l i t y i n the flow. In the long wave-length l i m i t the highest growth rate, of order 10 seconds, usually holds for the thermal i n s t a b i l i t y which arises whenever the cooling rate minus the heating rate has a negative deriva-tiv e with respect to temperature. This s i t u a t i o n arises i n the temperature ranges of 3x10 s tc 107 K. The growth rate of t h i s i n s t a b i l i t y i s d i r e c t l y related to the cooling rate, which i s proportional to the abundances of the elements present. If this i n s t a b i l i t y operates, one would expect s i g n i f i c a n t differences between stars of d i f f e r e n t composition. That i s , stars with higher CMO or metal abundances would have a greater cooling rate for temperatures exceeding 10 s K i n an o p t i c a l l y thin gas., As a r e s u l t the thermal i n s t a b i l i t y would grow on a shorter time scale. This may have some bearing on wolf-Hayet s t a r s , which appear to have higher CNO abundances than OE stars, and d e f i n i -tely have higher tass loss rates. An amplifying i n s t a b i l i t y for acocstic waves which i s usually present i s the simple growth of wave amplitude due to the density gradient i n the atmosphere. In a moving atmosphere t h i s occurs on a time scale of the gas speed divided by the scale height, about 2x10 3 seconds. The decceleration i n s t a b i l i t y of acoustic waves i s an absolute i n s -t a b i l i t y , . The growth rate i s j d v / d z ) - 1 , usually cf order 1000 71 t o 10* seconds. I f the gas i s ve r y h o t , g r e a t e r than 10 s K, t h e r e i s an a b s o l u t e i n s t a b i l i t y which o p e r a t e s cn the time s c a l e c f an hour. The r a d i a t i o n f o r c e p r o v i d e s an a m p l i f i c a t i o n f o r wavelengths i n t h e range 10 8-10 1 1 cm on time s c a l e s of 1200 seconds. As a r e s u l t o f t h i s work, a number o f s u g g e s t i o n s can be made f o r f u r t h e r i n v e s t i g a t i o n . L i n e v a r i a b i l i t y s h o u l d be pur-sued i n o r d e r t o u n r a v e l the n a t u r e of the v a r i a t i o n . High r e -s o l u t i o n o b s e r v a t i o n s , w i t h good s i g n a l t o n o i s e must be ob-t a i n e d . These o b s e r v a t i o n s s h o u l d e i t h e r be made i n s p e c t r a l r e g i o n s f r e e o f t e l l u r i c l i n e s , o r at a s u f f i c i e n t l y h i g h r e s o -l u t i o n t o n i n i f f i i z e b l e n d i n g w i t h the s t e l l a r l i n e . A l o n g e r segment o f t h e X-ray d a t a s h o u l d be a n a l y z e d t o c o n f i r m t h e model g i v e n f o r t h e a c c r e t i o n p r o c e s s , and sore im-p o r t a n t l y t o e s t i m a t e t h e d e n s i t y and v e l o c i t y as a f u n c t i o n o f ti m e . The t h e o r e t i c a l a n a l y s i s o f i n s t a b i l i t i e s can be extended to a l l o w a v e r t i c a l and h o r i z o n t a l wave v e c t o r , and a l l o w the wave motion t c have a h o r i z o n t a l as w e l l as a v e r t i c a l component and t h e n t o more g e n e r a l waves, such as a l l o w i n g v c r t i c i t y . B e s i d e s d y n a m i c a l g e n e r a l i z a t i o n s , d i f f e r e n t s o u r c e s p e c t r a s h o u l d be a l l o w e d , p a r t i c u l a r l y X-rays. T h i s would a l l o w t h e a n a l y s i s t o be extended t o g u a s a r s and n u c l e i of g a l a x i e s . The purpose of t h i s whole s t u d y i s t o a c g u i r e i n f o r m a t i o n on t h e fundamental p h y s i c a l n a t u r e c f t h e mass l o s s i n t h e pre-sence of a s t r o n g r a d i a t i o n f i e l d . My t h e s i s i s t h a t t h e wind i s observed t o be v a r i a b l e , and t h a t t h e v a r i a b i l i t y on time s c a l e s c f a day or l e s s can be a t t r i b u t e d t c i n s t a b i l i t i e s which 72 exist i n the wind. It i s suggested that the presence of these i n s t a b i l i t i e s changes the fundamental dynamics of the solutions to the flow of the s t e l l a r wind. The length and time scales established in the analysis w i l l allow the nonlinear equations of the flow to be attacked with a knowledge of their l o c a l be-haviour. 73 BIELICGBAPHY 2. Akhiezer, A. I. ana Polovin, B. V., 1971, Sov A PJjs^ Ospekii , 14, 278. 3. Aldrovandi, M. V. and Pequignot, D., 1973, Astr. Ag., 25, 137. 4. Aldrovandi, M. V. and Pequignot, D., 1976, Astr . A j . , 47, 321. 5. Earlow, M. J . and Cohen, M. , 1977, A\u00C2\u00A3.J.., 213, 737. 6. Eeals, C. S., 1929, J.N., 90, 202.. 7. Eeals, C. S., 1951, Puh^DjjO, i x , 1. 8. Eers, A., 1975, Plasma Physics. Ed. D e i i t t , C. and Peyraud, J. . (Gordon and Ereach: New York), p. 121., 9. Eerthcmieu, G., Provost, J . and Bocca, A., 1S75, A s t r . J j 3 . , 47, 413. 10. Eethe, H. A. and Salpeter, E. E. 1957. The Quantum Mechanics of One And Two Electron Atoms, (Sprinqer-verlag: New York). 11. Brucato, B. J. , 1971, M.N., 153, 435. 12., Burgess, A. and Summers, H. P., 1976, M.J., 174, 345. 13. Burgess, A. and Summers, H.P., 1969, A\u00C2\u00A3.J., 157, 1007. 14. Cannon, A. J . and Pickering, E. C., 1918, Ann.t of Obs A Q\u00C2\u00A3 Harvard C o l l . M vols, 91-99. 15. Cannon, C, J . and Thomas, B. N\u00E2\u0080\u00A2, 1977, k\u00C2\u00A3.J., 211, 910. 16. C a s s i n e l l i , J . P. and Castor, J . I., 1973, A D . J . , 179, 189, 17. C a s s i n e l l i , J . P., Olson, G. I., and S t a l i o , B., 1978, AJE.J., 220, 573. 18. Castor, J . I., Abbott, D. C., and Klein, B. I. , 1975, J E . J . , 195, 157. 19. Castor, J. I., 1977, in Colloquium Internationaux Eu CJBS N O J S . , 250. Mouvements Dans Les Atmospheres S t e l l a i r e s . ed. Cayrel, B. and Steinberg, M\u00E2\u0080\u00A2, (CNBS: P a r i s ) . 20. Castor, J . I., 1978, at IAU Symposium No. 8 3, Mass Loss And 74 21. Chapman, E. E. and Henry, R. J . 8., 1971 Ap..J., 168, 169. 22. Chapman, E..D. and Henry, R. J . B., 1972 Ap.*?., 173, 243. 23. C o n t i , P. S. and F r o s t , S. A., 1974, Ap.J. L e t t . , 190, L137. 24. C o n t i , P. S., 1978, A s t r . Ap. , 63, 225. 25. Cox, 0. E. and Tucker, B. H. , 1969, Ap,J., 157, 1157. 26. Dupree, A. K., 1968, Astrophvs. L e t t . . 1, 125., 27. Dysthe, K. B., 1966, Nuclear Fus i o n. 6, 215. 28. Hearn, A. G., 1973, A s t r . .Ap., 23, 97. 29.. Fahlman, G. G., C a r l b e r g , fi. G., and B a l k e r , G. A. H., 1977, Ap.J. L e t t . , 217, L35. 30. F i e l d , G. B. and Steigman, G., 1S71, A j . J . , 166, 59. 31. F i e l d , G. B. , 1965, JLp.J., 1\u00C2\u00AB2, 531. 32. Flower, D.,E. 1968, i n IAD Symposium NOj. 34., P l a n e t a r y Nebulae, ed. Osterbrock, D. E. and O ' D e l l , C. R. ( R e i d e l : D ordrecht), p. 205. 33. Hearn, A. G., 1972, A s t r . Ap.., 19, 417., 34. Hearn, A. G., 1975, A s t r . Ap.. , 40, 355., 35. Henry, B. J . W. 1970, Ap.J., 161, 1153. 36. Hutchings, J . B. , 1976, Aja. J . , 203 , 438. 37. Jackson, J . C. , 1S75, JS.J. , 172, 483. 38. Johnson, L. C. , 1972, Ap* J . , 174, 227. 39. Jordan, C , 1969, M.J. , 142, 501. 40. Kato, T., 1S76, Ajp.. J . Supji. , 30, 397. 41. Kippenhahn, E., Mestel, L., And Perry, J . J . , 1975, A s t r . Ap., 44, 123. 42. K l e i n , R. I . and C a s t o r , J . I . , 1978, Ap.J., 220, 902. 43. K r o l i k , J . H., 1977, Phys. F l u i d s 20, 364. 44. l a c y , C. H., 1977, Ap.J., 212, 132. 75 45. Iamb, H, , 1945, Hydrody namies. (Dover: New York). 46. Earners, H. J, G..L. M. and Morton, D. C , 19*76, A\u00C2\u00A3.J. Su\u00C2\u00A3\u00C2\u00A3., 32, 715, 47. lamers, H.J. G, L. M., and Snow, 1 . P., J r . , 1978, A\u00C2\u00A3.J., 219, 504. 48. Lctz, W. , 1967, J f i . J . SupjS., 14, 207. 49. Lucy, L. B. And Solomon, P. M., 1970, j\u00C2\u00A3.J., 159, 879., 50. Lucy, L. B. , 1971, Ac.. J., 163, 95. 51. McWhirtler, E. 8. P., 1975, i n Atomic .And, Molecular \u00C2\u00A3Efi\u00C2\u00A3J5se,s In Astrophysics, ed. Hubbard, M. C. E, and Nussbaumer, H (Geneva Observatory: Sauverny, Switzerland), p. 205. 52. flestel, I., Moore, D. W., and Perry, J . J., 1976 Astr* lj\u00C2\u00BB, 52, 203. 53. Mewe, fi., 1972, Astr. Ap_. , 20,215. 54. Mihalas, D., 1972, Non-LTE Model Atmospheres f o r B and 0 Stars, NCAB-TN/STB-76, (National Center for Atmospheric Research: Bculder). 55., Milne, E. A., 1926, M. JJ., 85, 813. 56. Moore, C. E,, Minnaert, M. G. J., and Houtgast, J . , 1966, The Solar Spectrum. NBS Monograph 6 1. 57. Moore, C. E,, 1949, Atomic Energy l e v e l s , NBS Circular No. 467. 58. Morton, D . C and Smith, \u00C2\u00BB. H., 1973, Ajg.J. Su\u00C2\u00A3\u00C2\u00A3. , 26, 333. 59. Morton, D . C , 1967, Aj.J., 147, 1017. 60. Mushotzky, R. F., Solomon, P. M., And Strittmatter, P. A., 1972, Ap_. J. , 174, 7. 61. Nelson, G. D., and Hearn, A. G., 1978, Astr. An.. , 65, 223. 62. Pounds, K. A., Cooke, 6. A., Ricketts, M. J . , Turner, M. J. , and E l v i s , M., 1975, M.N., 172, 473. 63. Raymond, J . C , Cox, D. P., and Smith, B. H., 1976, J. \u00E2\u0080\u00A2 204, 290. 64. Bcsendahl, J..D., 1973, Afi.J., 182, 523., 76 65. R y b i c k i , G. B. and Hammer, D. G., 1978, Ap.J., 219, 654. 66. S c h r e i e r , E. J., S c h w a r t z , K., G i a c c o n i , R. , F a b b i a n c , G., and M o r i n , J . , 1976, A j . J . , 204, 539. 67. S e a t c n , M. J. 1958, E g V i JSod[A Phj^s., 30, 979. 68. S e a t c n , M. J., 1959, M.N., 119, 81. 69. S i l k , J. And Brown, B. L., 1971, A j . J . , 163, 495. 70. Snow, T. P. J r . and Morton, E. C , 1976, A p . J . Supja., 32, 429. 71. Snow, T. P. J r . , 1977, Ap . J . , 217, 760. 72. S o b o l e v , V, V., 1960. Moving En v e l o p e s of S t a r s * (Harvard U n i v e r s i t y P r e s s : Cambridge), 73. S p i t z e r , L, J r . , 1962, P h y s i c s oJ F u l l y . I o n i z e d Gases. 2nd ed., ( I n t e r s c i e n c e : New York) \u00E2\u0080\u00A2 74. Steigman, G., Werner, M. 8., and Ge l d o n , F.... B. , 1971, Ap.J. , 168, 37 3. 75. Summers, H. P., 1974, M.N., 169, 663. 76. Sunyaev, B. A. and V a i n s t e i n , L, A., 1968, A s t r o p h v s . l e t t ., 1, 193. 77. Thomas, B. N., 1973, A s t r . Ajg., 29, 297. 78. U n d e r b i l l , A. B. , 1960, i n JcJU o f SJbars A^d S ^ e l l a j c Systems, S t e l l a r Atmospheres. ed.,J. G r e e n s t e i n , (ti. .of C h i c a g o : C h i c a g o ) , p. 411. J r 79. Walker, G, A. H., B u c h o l z , V., Fahlman, G. G., G l a s p e y , J . , La n e - W r i g h t , D., flochnacki, S., and C o n d a l , A., 1S76, Pro c , IAU C o l l o q u i u m 40, ed. M. Duchesne, ( r e i d e l : D o r d r e c h t ) . 80. H i e s e , W. L., S m i t h , M. W., and Glennon, B. M., 1966, NSBDS-NBS4. 81. Wiese, W. L., S m i t h , M. \u00C2\u00AB., and M i l e s , B. M., 1969, NSBDS-NBS22. / 82. W i l s o n , B., 1962, J.O.S.B.T. 2. 477. 83. W r i g h t , A. E., and Barlow, M. J . , 1975, M.N., 170, 41. 84. Y c r k , D. G. , V i d a l - M a d j a r . A., L a u r e n t , C , Bonnet, fi. , 1S77, Ap.J. L e t t . , 213, 161.. 77 APEENDIX 1. SUPEBSGNIC ACCBETTON This appendix i s a reprint of the paper e n t i t l e d \"Badiative Effects i n Supersonic Accretion\", which appeared i n the A\u00C2\u00A7*rcjhxsical Journal Volume 220, p. 1041. The theory developed was used in Chapter III tc deduce the correlation between the wind velocity and wind density i n the observed i n t e n s i t y t r a n s i -tion cf Cen X-3. / T H E A S T R O P H Y S I C A L J O U R N A L , 220:1041-1050, 1978 March 15 \u00C2\u00A9 1978. The A m e r i c a n A s t r o n o m i c a l S o c i e t y . A l l r i g h t s r e s e r v e d . P r i n t e d i n U . S . A . R A D I A T I V E E F F E C T S I N S U P E R S O N I C A C C R E T I O N R. G. CARLBERG Department of Geophysics and Astronomy, University of British Columbia Received 1977 March 21; accepted 1977 September 22 A B S T R A C T Supersonic gas flow onto a neutron star is investigated. There are two regimes of accretion flow, differentiated by whether the gas can cool significantly before it falls to the magnetosphere. If radiative losses are negligible, the captured gas falls inward adiabatically in a wide accretion column. If the radiative energy-loss time scale is less than the fall time, the gas will cool to some equilibrium temperature which determines the width of the wake. A n accreting neutron star generates sufficient luminosity that radiation heating may determine the temperature of the accretion column, provided the accretion column is optically thin. Gas crossing the shock beyond the critical radius forms an extended turbulent wake which gradually merges into the surrounding medium. As a specific example, the flow for the range of parameters suggested for the stellar wind X-ray binaries is considered. Subject headings: shock waves \u00E2\u0080\u0094 stars: accretion \u00E2\u0080\u0094 X-rays: binaries I. INTRODUCTION Recent observations of X-ray binaries, at both optical (Conti and Cowley 1975; Dachs 1976) and X-ray (Jones et al. 1973; Pounds et al. 1975; Eadie et al. 1975) wavelengths, show phase-dependent absorption of radiation. It has been suggested that this is caused by a wake trailing the compact object which emits the X-rays. Models of the wake based on the X-ray observations were put forward by Jackson (1975) and Eadie et al. (1975). The general problem of a gravitating body moving through a gas at a velocity much greater than the sound speed was first discussed by Hoyle and Lyttleton (1939). More recently, wakes were discussed by Davidson and Ostriker (1973), Illarionov and Sunyaev (1975), and McCray and Hatchett (1975). These models are incomplete in that they lack a description of the gravitationally perturbed gas which is unbound, i.e., the far wake. Although most of these papers emphasize the importance of radiative effects, no clear analysis has been made of the variations in the flow of gas caused by radiative gains and losses. In this paper supersonic accretion onto a neutron star is considered. There are three basic physical parameters: the mass of the accreting body and the free stream velocity and density of the gas. The dynamics of the flow are essentially determined by the free stream velocity and the mass. The angular width of the accretion column depends on its temperature, which in turn is regulated by radiative cooling and heating and is sensitive to the gas density. The proposed description is worked out for linear motion, which is a good approximation for an accretion radius much smaller than the system dimen-sions. A schematic of the model is shown in Figure 1. The important regions are labeled: (1) the incoming supersonic gas; pressure forces can be neglected and streamlines are taken to be coincident with particle trajectories in a gravitational field; (2) the shock-heated sheath where the incoming gas impinges on the accretion column; the transverse component of the velocity is rapidly halted, providing pressure to contain the accretion column; (3) the accretion column, in which gas falls inward, toward the accreting body; (4) a region of spherically symmetric flow which may exist near the accreting body; (5) the base of the accretion column; beyond this the flow is regulated by the physics of the magnetosphere around the accreting object; (6) the accreting body, where the kinetic energy of the gas is liberated at a surface shock; and (7) the far wake, several hundred times the length of the accretion column. The density contrast between the far wake and the surrounding medium gradually goes to zero. One major qualitative aspect of this model is that there is no bow shock standing off from the front of the body which is distinct from the tail shock. U n -doubtedly there will be a preceding shock, but pressure waves generated there will not propagate very far in a transverse direction because the streamlines of the flow are bent in by gravitation. Consequently, the bow shock merges into the tail shock. In general, this will be the case for any body whose size is less than the \"accretion radius\" RA = 2GM/VQ2, where V0 is the free stream velocity. Calculations by Hunt (Eadie et al. 1975) indicate that part of the infalling column may \"miss\" the accreting body and force the leading shock forward. This occurs because small nonradial velocities increase toward the body, by conservation of angular momentum. This effect will be ignored. This model is to be applied to a neutron star orbiting a massive star with a strong stellar wind. For convenience, scaled variables will be used for the distance rv = rjRA: free stream density, w u = HQ /10 1 1 c m \" 3 ; free stream velocity, VB = V0j 103 cm s~ 1; and 1041 1042 CARLBERG Vol. 220 F I G . 1 . \u00E2\u0080\u0094 A schematic of supersonic accretion gas flow showing: (1) the incoming supersonic gas; (2) the shock-heated sheath; (3) the accretion column; (4) the possible spherically symmetric inflow at the bottom of the column; (5) the Alfven surface at the bottom of the column; (6) the accreting body; and (7) the far wake. mass of the accreting body, m \u00E2\u0080\u0094 Mx/M0. Similarly, in later calculations, the temperature of the column will be represented as T6 = T/IO6 K. This form of notation will be used throughout. II. D Y N A M I C S O F T H E G A S F L O W A gravitating body is placed in a uniform stream of gas moving at some velocity V0. To the point where the gas crosses the tail shock, we assume that the stream-lines of the flow can be found from particle dynamics, i.e., the flow is dominated by inertial forces. The velocity can be obtained from the equations of conservation of energy and angular momentum, Wr2 + V,2) GM Wo2 and V \u00E2\u0080\u0094 V0 \u00E2\u0080\u00A2 0) where Vr and V# are, respectively, the radial and tangential components of the gas velocity relative to the accreting object. The trajectories are given by (Ruderman and Spiegel 1971), I = 5* 0 + C0S & + ISin * ' (2) where s is the impact parameter and is the angle measured from the accretion axis. From these equa-tions, Danby and Camm (1957) obtain the density n = n(r, ) as n = 2 sin 0/2 2 + sin-+ sin2 - + sin 2 r A + sin2 (3) As a simplifying approximation, we take the case of (f> small and 2 \u00C2\u00AB RA/r, to obtain the relations, and V* \u00C2\u00AB V0(RJry*, 20 (4) d) The Sheath The gas crossing the shock has a discontinuity in its motion described by the equations for the conserva-tion of the total energy and of the normal components of mass flux and momentum. Assuming that a strong shock occurs (good for Mach numbers greater than ~3), the postshock density and temperature are (for a ratio of specific heats y = 5/3) P2 = 4Pl, and r a = ^ -^. (5) where R is the gas constant and 1 and 2 refer to the pre- and postshock conditions, respectively. Vn (x V^) is the component of velocity normal to the shock. An important point is that specific energy is conserved across a shock. If the gas is energetically unbound ahead of the shock, it will remain unbound behind the shock in the absence of cooling. On the other hand, if all the thermal energy is immediately lost, one finds that the gas is energetically bound for all radii less than RA. The sheath is bordered by the shock on the outside and the inward-flowing gas of the accretion column on the inside. The sheath is a dynamically defined region where the gas slows to a stop, changes direction, and joins the accretion column. 80 No. 3, 1978 RADIATIVE EFFECTS 1043 The semiangle to the shock cone will be approxi-mated by the semiangle of the accretion column. To this end we demonstrate that the thickness of the sheath is small for the case of a narrow shock cone. The radial-velocity component is approximately parallel to the shock and is continuous across the shock. In the limit of Vr = V0, one easily finds that gas entering the shock sheath at r 0 < RA will travel to a maximum distance r given by r0 (6) before the radial velocity is brought to zero. The gas would then join the accretion column. An estimate of the sheath thickness can be obtained by equating the mass influx between r0 and r, (r x r0, r \u00C2\u00AB RA) -nn0V0r02, to the mass flux through a cross section at distance r, H S K 02T7T S H', where vt' is the thick-ness of the sheath, ns the postshock density at r, and s the angle to the shock. This gives the semiangular width of the sheath as w r0 ar (7) Thus the maximum width of the sheath is only a function of distance from the accreting object, and for r \u00C2\u00AB RA the sheath width will be negligible. The above calculation assumed laminar flow and no premature mixing of sheath gas into the column, whereas it is quite likely that the sheath is turbulent. The Reynolds number in the sheath is 4 . 1 x 1 0 4 / j 1 1 r \u00E2\u0080\u009E 2 $ ~ 1 K 8 - 6 , indicating the possibility of tur-bulence. A turbulent sheath would come into equi-librium with the column more rapidly than laminar flow through mixing. As a consequence, a turbulent sheath would be even thinner than the limit set in equation ( 7 ) . b) The Accretion Column The mass flux in the accretion column is simply dM/dt = irp0V0sc2, where sc is the critical impact parameter, taken as the impact parameter of the streamline which would have a total energy of zero on the accretion axis. To allow for the thermal energy, a parameter (Z is introduced, such that the \" true accre-tion radius\" is equal to fiRA. In principle, jS is deter-mined once the physical parameters, the density, velocity, and accreting mass, are specified. The param-eter jS will be taken as the ratio of the specific kinetic energy (%V02) to the specific enthalpy (5RT0) if (\u00E2\u0080\u00A2JrKp2 < 5RT0), otherwise 0 = 1 , where T0 is the equilibrium temperature of the column at r\u00E2\u0080\u009E = 1 . The accretion rate is then dM/dt = 3 . 6 5 x 1 0 1 6 / J H / M 2 K 8 - 3 g s\"1 . ( 8 ) For a column in equilibrium, the transverse momen-tum of the incoming gas must be balanced by thermal pressure in the column, From this one obtains a relation between the central temperature of the column and the opening angle, V02 c 4(2Y'2Rp 2 . 1 3 x l O 7 ^ \" 1 K r a d \" 1 . ( 1 0 ) Pressure forces are unable to support the gas, and it falls inward toward the accreting object down the accretion column at a velocity v = (GM/r)112. Using equation ( 1 0 ) and mass conservation, we find that the equilibrium accretion column density is, for fiRA \u00C2\u00BB r, nc = 6 . 4 0 x 1 0 1 3 r 6 - a r t , - 3 ' a \u00C2\u00AB 1 1 K 8 * / 3 - 1 cm\"3 ( 1 1 ) The assumption that the width of the column is maintained by gas pressure is justified by the required default of any stronger forces, turbulence in particular. One can do a pressure confinement calculation similar to the one above by assuming a fully turbulent accre-tion column. The internal pressure in the column would be generated by the turbulent velocity, which can be taken to be a fraction/ of the velocity of fall. Requiring that the opening angle of the column be less than, say, 1 radian, we find that f is restricted by 2V2 r (12) Pi i (9) This implies that the turbulent velocity must become a smaller fraction of the fall velocity as it nears the neutron star; otherwise the turbulent pressure is im-possible to contain. But the Reynolds number of the gas increases inward (except for adiabatic infall), and we would expect the turbulence, if present, to increase and disrupt the column. Therefore, if the column exists, it must be in laminar flow. There are several reasons to think that laminar flow can obtain in the column. The turbulence would probably originate in the \"shear layer\" between the sheath and the column, but the Reynolds number in the sheath decreases down the column. In addition, the gas is being strongly accelerated only in the radial direction, which does not provide a driving force for turbulence. III. RADIATIVE EFFECTS If the gas is unable to cool before joining the accretion column, the column gas will fall adiabatically and will resemble the accretion scenario found by Hunt ( 1 9 7 1 ) , i.e., a very wide accretion column trailing the accreting body. Note that Hunt's solutions were obtained with essentially zero pressure at the boundary of the accreting body, and that the accretion rate would probably be diminished by the nonzero base pressures of a magnetospheric shock above the neutron star. On the other hand, if the gas cools much faster than any time scale for movement, a cold, narrow, high-density column will be formed. In order to determine which regime prevails, we compare the time scales for radiative energy-loss mechanisms with the time scale for infall of the gas, which is the basic and only uniquely identifiable dynamic time scale of the 1 0 4 4 CARLBERG Vol. 220 problem. The fall time from the accretion radius is approximately, '^ = 3 V ( G M ) = 2 5 0 R U 3 , 2 K R 3 / M S - \u00C2\u00B0 3 ) a) Cooling Time Scales In the absence of any heating, the temperature of the gas is entirely dependent upon whether or not a significant amount of cooling can take place in the gas before it reaches the surface of the accreting body. In this section an estimate is made of the cooling time scale, which divides the density-velocity parameter space into regions of cooling and no cooling. In the following, all radiative time scales will be defined as 3kT divided by the appropriate heating or cooling rate, where k is Boltzmann's constant. When the gas crosses the shock, the ions get most of the thermal energy, since they have a much shorter mean free path than the electrons. The electron-ion equilibrium time (Spitzer 1962) in the sheath is, with n = 4\u00C2\u00AB0, a maximum of r e q = 50.6^-^-^ s. (14) This time is compared with the fall time and is plotted in Figure 2. For the postshock gas, the equilibrium time decreases with density at the same rate as the cooling time and is always shorter than it. Therefore the postshock gas comes into collisional equilibrium and the electron and ion temperatures are assumed equal. The cooling from the postshock temperature can be taken from Figure 1 of Cox and Daltabuit (1971), omitting the cooling due to forbidden and semi-forbidden lines. The postshock cooling is assumed to be unaffected by any radiation present. Two assump-tions are made for the density of the postshock gas in the sheath. The rightmost line is drawn for the mini-mum possible postshock density, 4n0. This is an underestimate, since the density increases from its free stream value toward the accretion axis, by approxima-tion (4). The angle to the shock decreases with the temperature by equation (10), and choosing the minimum temperature in the column as 104 K results in the cooling line on the left. Gas flows with densities and velocities in between these two cooling lines may be subject to an instability from the cooling to noncooling state and vice versa. If hot, uncooled gas mixes into the accretion column and expands it such that the shock moves outward, it will decrease the density of the incoming gas, by approximation (4). If the density drops sufficiently, the incoming gas may no longer cool and the column will expand to its uncooled state. Consequently we take the rightmost cooling line (labeled 4n0) as the effective cooling line. A point of interest is that, for gas crossing the shock at a distance of less than 3 x 1010 cm, the postshock temperature is greater than 107 K, for which brems-strahlung is the dominant cooling mechanism, until the gas is close enough (see eq. [25]) to be Compton cooled. The time scale for bremsstrahlung losses varies with \u00C2\u00AB-lr1'2, which remains constant with distance in the sheath, whereas the dynamic time scale is decreasing as r 3' 2. Consequently, even though gas entering the column at large radii may cool, lower down the gas in the sheath may remain hot. As shown in the Appendix, the sum of the pressure force for the postshock gas FIG. 2 . \u00E2\u0080\u0094 T h e cool ing diagram constructed with the free stream density and velocity. Solid lines, regions of cool ing for the max imum (IO 4 K ) and m in imum (4n0) densities. Coo l i n g occurs to the right, i.e., higher densities, of these lines. Below the dashed lines (teq < '/) the electron and ion temperatures are equal. The hatched region is heated to the C o m p t o n equi l ibr ium temperature and is certainly subsonic, whereas below the dot-dashed line (/Wc < 1) the M a c h number of the incoming gas based on the C o m p t o n heating rate is less than 1. 82 No. 3, 1978 R A D I A T I V E E F F E C T S 1045 and the gravitational force is still directed downward, but eventually the excess energy must be lost if the gas is to be gravitationally bound to the neutron star. One way to do this would be through turbulent mixing at the boundary between the upward-flowing sheath gas and the downward-flowing accretion column. This would decrease the effective cooling time for the lower sheath by diluting the hot gas with the cooler, denser gas of the column. If a mixing process is required in order to capture gas entering the column at small radii, it would imply that, if gas at the top of the column is unable to cool, then the accretion may become very inefficient. b) Heating When the accreting body is a neutron star, the accretion luminosity may be sufficient to cause significant heating of the infalling gas. (An additional problem is that the incoming gas stream may be heated to a sufficiently high temperature that the assumption of supersonic flow is invalidated. This is considered in the Appendix.) Approximate rates for photoionization and Compton heating are derived and used to con-struct a heating diagram similar to the previous cooling diagram. Optical depths must also be considered, because radiative heating will be impossible if the optical depth up the column becomes too large. The accreting object is assumed to be a neutron star of radius 106cm. The entire kinetic energy of the infalling gas is converted into radiation at the surface shock. The resulting luminosity is L = 4.70 x 10 3 6/; 1 1K 8- 3w 3(^/10 6cm)- 1ergss- 1. (15) An upper limit to the gas temperature in the column due to photoionization heating is required. The calcu-lations of Hatchett, Buff, and McCray (1976). show that, for log \u00C2\u00A3 greater than 2, where \u00C2\u00A3 = L/nr2, the CNO elements are completely ionized. The heating will then be limited by the total recombination rate to the ground state, so the limiting photoionization heating rate is neafx/3, where ne will be taken to be the density in the column from equation (11), a is the recombina-tion rate to all levels for completely ionized oxygen, / (=10 - 3) is the fractional abundance of oxygen, increased slightly to allow for some carbon and nitrogen, and x/3 is the average energy deposited per ionization for a v~l spectrum. (The spectrum may not be but all spectra deposit an average energy of order x-) The recombination rate used is the expression given by Allen (1973) for the total recombination rate to the ground level, \u00C2\u00AB = 3 x 10-loZ2r-3'4. The resulting heating time is fph = 22.9T615'V'2\u00C2\u00ABir1>V4j81 x (//10-3)-HZ/8)-4s. (16) Comparing this with the fall time, we find that the temperature is limited by T6 < 1.89\u00C2\u00AB n 4 / 1 5V 1 5 (//10- 3 ) 4 ' 1 5 x (Z/8) I 6 / 1 5 . The expression used for photoionization heating applies only if heating in a static gas would be able to attain this temperature and the absorption and scattering optical depth up the column are less than 1. The static condition, log \u00C2\u00A3 > 2, is equivalent to VB < \.02Te2l3rv~116, which is indicated in Figure 3, and is always satisfied in the cooling region. The photo-ionization absorption can be estimated from the equations of ionization balance and of optical depth, Le-and w ( D = ^ p ^ G dr dr (18) where ne, nt, and na are the number densities of electrons, ions, and ground-state absorbers, respec-tively. The integrals in the exact formulae have been approximated by quantities integrated over frequency, and it is assumed that one ionic species is doing the major part of the absorbing at any given temperature. In the neighborhood of 106 K, the major absorber is O VIII. Setting / as the fraction of atoms that are absorbing X-rays, we find a solution to the above equations similar to Mestel's (1954), = i*4 Ufa(T) f ne (19) The bottom of the column, rb, has been assumed to be the magnetosphere at a distance near 108 cm from the center of the neutron star, and it is assumed that there is no significant opacity between the source of radiation and this lower boundary. This is consistent with the magnetospheric model of Arons and Lea (1976). Using the total recombination rate and 5 keV as an average photon energy, we find that the above integral becomes 1 0.9897V1 9'4 ln ( r / r ^ - 3 n u F 8 5 x (Rx/106 cm)(//10-3)(Z/8)4 . (20) For numerical estimates, the distance dependence of the optical depth will be ignored [In (r/r\u00E2\u0080\u009E) ~ 1]. This optical depth is meant to be useful only as an indica-tion of how radiative heating is attenuated. As it turns out, this estimate of the photoionization optical depth is always less than the electron scattering optical depth for the range of parameters plotted. A more precise calculation is required to estimate the trans-mitted spectrum as a function of frequency. The other major source of heating in an X-ray illuminated gas is Compton scattering. The Compton heating rate is given by Buff and McCray (1974) as - a-^kT L <7r me 2 the approximat ion used for the heating rate would be val id for a static gas. A s the flow parameters cross the T S > T C (dotted) line, the optical depth up the sheath exceeds that of the co lumn. and a describes the shape of the spectrum (a = 1.04 for a blackbody and i for an exponential spectrum). In principle, a and e are determined by the physical parameters n0, V0, Mx, and the radius of the accreting object. In view of the complexity of the calculation of the emitted spectrum, we chose to leave them as parameters. As typical values we chose a = \ and 0.4767'6rI)1'a)3-1m-2(i?x/106)-1 . (23) If the base of the column becomes optically thick, the radiation will be thermalized, so that Compton heat-ing becomes negligible as a result of the e parameter's being reduced. As discussed by Felten and Rees (1972), spectrum alteration begins when the optical depth T * = (3T F F T C S ) 1 / 2 exceeds 1, where r f f and T c s are the free-free and electron scattering optical depths. The calculation indicates that the optical depth at a photon energy of 5 keV is always much less than 1 as long as the flow is supersonic. As one gets closer to the source of radiation, the time scale for Compton heating decreases faster than the fall time. The bottom of the accretion column may be heated to Compton equilibrium, even though the regions higher up may be in photoionization equi-librium, or optically thick and cold. From equation (10) we see that the column is impossible to contain for VB < 0.9, and the column will widen and may even become spherical at the bottom. For electron scattering the new effective base of the column is at a distance where the Compton heating and the fall time scales are equal, r8 = 1.39/;112K8-6(\u00C2\u00AB/0.5)-2(\u00C2\u00A3/5 keV)~202. (24) The spherically infalling material below this has a negligible contribution to the optical depth, because the density is reduced by the much greater column angle. The electron scattering optical depth along the column is rc = 9.18\u00C2\u00AB11K87V2(/-6/108 cm)\"1'2. (25) The actual line plotted on Figure 3 for the electron scattering opacity assumes that the effective base of the column is at the Compton heated distance of equation (25) and that the column temperature is determined by the photoionization heating temperature of equation (17). The range of photoionization heating is thus extended by the reduction of column opacity. If one computes the Alfven radius from B2/8n = l/2py2, one finds that in some cases the Alfven radius exceeds the Compton heated radius and will determine the effective base of the flow. But these cases turn out to be in the region of the density-velocity diagram above the cooling line and hence of little interest to the heating calculation. The column is effectively optically thin to the side-ways loss of radiation because the sideways optical No. 3, 1978 RADIATIVE EFFECTS 1047 depth of the column is dominated by electron scatter-ing and is always less than 1 for the range of parameters plotted. The higher-energy X-rays will be attenuated by K shell absorption by elements with ionization potentials greater than the CNO elements. For a spectrum with a typical photon energy of 5 keV, the K shell cross sections of Daltabuit and Cox (1972) and the abun-dances of Allen (1973) suggest that the dominant K shell absorber will be silicon. The calculations of Hatchett et al. indicate that a typical silicon atom, for log \u00C2\u00A3 > 2, will have several electrons left, and therefore will have a cross section of order 10\"19 cm2, which is relatively independent of temperature and radiation flux. Combining this with a fractional abundance of 3 x 10~5, the effective cross section at the absorption edge is only 4.5 times the electron scattering cross section. Photons below the edge will be less affected, primarily interacting with only lower-abundance magnesium, and for those above, the cross section decreases approximately as E~3, until another edge, due to low-abundance sulfur, is encountered. In general we expect that K absorption will be of the same magnitude as the electron scattering. Similarly, the K shell photoionization heating rate is different from Compton heating only by a multiplicative factor of order 1, and will be ignored. c) Accretion Scenarios In Figure 3 a box has been drawn which encloses the suggested range of wind densities and velocities for stellar wind X-ray sources Cen X-3 and 3U 1700 \u00E2\u0080\u0094 37. Only a small part of this region is subsonic and beyond the description given here. For a given density and velocity, it is possible to qualitatively describe the flow. If the density and velocity parameters of the free stream lie above the cooling line, the accretion may be less efficient as pressure forces in the hot gas of the sheath become more important. This would be re-flected in a diminished luminosity. In general, flows with parameters above the cooling line would broadly resemble the scenario found by Hunt (1971). Captured gas falls inward with its temperature rising adia-batically. Below the cooling line, the gas temperature drops to some equilibrium value and falls down the accretion column. Although the K absorption edges will alter the spectrum somewhat, the line above which the electron scattering opacity exceeds 1 is almost co-incident with the cooling line, so that heating of the most distant parts of the accretion column will be possible below the cooling line. Typical maximum temperatures for Ve \u00E2\u0080\u0094 1 are Te = 1 at nn \u00E2\u0080\u0094 0.1 and Te \u00E2\u0080\u0094 3.5 at /in = 10. The column semiangles are 2\u00C2\u00B0 and 9\u00C2\u00B0, respectively. The gas will remain at the equilibrium temperature specified by the local radi-ation field, since all radiation time scales are more rapid than the fall time. Near the source Compton heating dominates, the temperature rises to the Compton equilibrium value, and the column expands so that the infall becomes almost spherical. It is of particular interest to compare the electron scattering opacity up the column with that up the edge of the sheath in the postshock gas. The opacity up the sheath is Evaluating this integral, and choosing (somewhat arbitrarily) the maximum extent of the column to be the distance at which the density has dropped to 4/3n0, i.e., rm = RAl(22), gives r s = 2 . 2 6 \u00C2\u00AB n K 8 2 r 6 - 2 , (27) whereas the electron scattering opacity up a column with a Compton heated base is 7.79V^TQ' 2 . A.S a result we find that the opacity up the sheath is greater than up the column if \u00C2\u00AB U K B \" S > 3.45, a result which is independent of the column temperature. This line has been included in Figure 3. If the stellar wind in which, the neutron star is em-bedded has velocity and density variations, this analysis predicts potentially observable effects. The most obvious is that, if the line-of-sight optical depth is constant, the X-ray luminosity responds to variations in n0V0~3 (eq. [16]) on times of variation longer than about two fall times, or 500VB~3 s. This variation reflects the local structure of the wind for regions of size greater than 2RA (5 x 10 1 0K 8~ 2 cm). Another expected effect is that, as n0V0~2 increases, the optical depth up the sheath will exceed that up the column. Thus, if the X-ray source is occulted by the accretion column, the X-ray absorption would change from a single dip (Tc > rs) to a double dip (TC < Ts). In addition, Jackson's calculations (1975) indicate that, if the gas fails to cool, the absorption up the sheath always dominates. IV. THE FAR WAKE The Reynolds number of the gas flow is extremely high (V0RAjv = lO 1 2 ^- 1 /?! 1 ^- 5 ' 2 / ; !! . ) , and the far wake is expected to be turbulent. Turbulence in super-sonic flows is not well understood, but experimental studies of supersonic wakes (Demetriades 1968) indicate that a phenomenological theory as outlined by Townsend (1976) provides a reasonable description of supersonic far wakes. Unfortunately, the dynamics of laboratory wakes are not dominated by a gravita-tional field, and therefore the applicability of the description to this case must be carefully considered. The subsonic theory is based on the observation that the wake remains self-similar with respect to a characteristic velocity and length scale. This is com-bined with the momentum equation from which all small terms have been dropped. The axisymmetric far wake is found to be self-similar with respect to the 1048 CARLBERG Vol. 220 half-width, /, and the turbulent velocity scale, u, denned by I \ (\u00C2\u00B1\113(_r_Y'3 where the so-called momentum radius RM (the radius such that the drag force is 1 /2pV02TrRM2) has been taken to be RA. RT is the turbulent Reynolds number, observed by Demetriades to be 12.8. The width scale of the wake implies that the small angle approximations for the density and transverse velocity apply for the exterior supersonic flow. Hence the effective exterior pressure will be the transverse momentum flux, which varies with distance along the wake. But the equations of momentum used to derive the length and velocity scale assumed that there was no pressure gradient in the free stream. For the gravitational wake, an order-of-magnitude estimate of the pressure gradient term in the small angle approxi-mation-finds that it is of order u2/l(RA/l), whereas the retained terms in the momentum equation are of order u2jl. Consequently, for RA \u00C2\u00AB I, the pressure gradient term can again be dropped. Using the half-width of equation (28), we find that RA/1 is of order (RJr)113. Thus the equations are consistent for r \u00C2\u00BB RA, but the crucial observation that a gravitational wake is self-similar is unavailable. Experiments also indicate that the flow may not be self-similar for distances of several tens of the momentum radius, but the deviation of the turbulent velocity scale from the self-similar value is a factor of 2 or less. For sufficiently low Reynolds numbers, part of the wake may be in laminar flow, ln this case the velocity defect on the axis is u/V0 = 3/2RA/r, and the half-width varies as (rRA)xlz (Lamb 1924). The Reynolds number grows with distance, and the flow will even-tually become turbulent. With the extreme Reynolds numbers present here, the wake is expected to become turbulent within the sheath of the accretion column. If the gas in the wake has no energy losses, i.e., 5RT + \v2 constant, the temperature on the axis is found to be T-Tn = ^ = 3.55 x }0eVe2l3m2l3rl2-2'3 K , (29) where TK is the temperature of the gas external to the wake. This temperature implies that the turbulent velocity scale is subsonic. The temperature becomes equal to 2Tm at R^Volc^, where RB is the Bondi radius, 10!4(104 K./T) cm, and cx is the sound speed at Tm. Experimentally it is observed (McCarthy and Kubota 1964) that the pressure is approximately constant across the wake. Equating pv'2 at the wake boundary to the gas pressure at the center gives the central density \u00C2\u00BB\u00E2\u0080\u009E = 1.06 x 1 0 1 1 n 1 1 r o 1 \" % - 1 ' 3 ' - i 2 ~ 1 ' 6 . (30) Note that this density is greater than that which would be found by using the external static pressure by a factor of ^f- = 376.07-4-1 J V ' 3 r i a - 5 / 6 \u00E2\u0080\u00A2 (31) The temperature estimate and density estimate of equations (30) and (31) assume that the wake is iso-energetic, but at these temperatures and densities, radiation cooling can be significant. Time scales of interest are the cooling time (where A is the cooling coefficient) 1 2 . 4 \u00C2\u00AB 1 X \" 1 K 8 5 / 3 / - 1 2 - 1 / 2 ( 1 0 - 2 2 ergs cm3 s\"7A) s, (32) turbulent dissipation of kinetic energy time scale, 3RT/(u3/l) = 3.04 x 10 4 K 8 1 / 3 m- 2 ' 3 r 1 2 5 ' 3 s, (33) and the turbulent time scale, 2.40 x 103 K 8 - 1 / - 1 2 m 2 ' 3 s . (34) These time scales imply that, for large distances, cooling removes most of the thermal energy from the wake. If the sound speed within the wake drops below the turbulent velocity scale, the turbulence would then become supersonic, leading to shocks which rapidly heat the gas, but the shocks would occur on the basic turbulence time scale, and would not be able to reheat the bulk of the gas. One might speculate that the temperature would decline to the minimum of 104 K, but with an extremely clumpy-distribution. If cooling is complete, the simple model used, which does not consider the energy budget, may break down com-pletely. Its value lies in the fact that, as the gas cools, the Reynolds number becomes even greater, and the dynamics of the gas flow in the far wake are almost certainly dominated by turbulence. One can combine the density and width to show a column density across the wake of sufficient size to produce optical absorption of radiation from the primary. That is n22l = 2 x 1033 c m - 5 for an optical depth in the wake of one at Ha, assuming the lower level is populated by recombinations at 104 K and depopulated by radiative transitions. This would be possible whenever the wake was silhouetted against the primary star. But these simple considerations fall well short of the ability to reproduce line profiles as seen by Conti and Cowley (1975). The wake will remain cold in the presence of an X-ray source for L/nr2 < 10, or distances from the X-ray source of r 1 2 > 3L371'2/;!!112. The absorption cross section for X-rays by a cold gas of cosmic abundances is about 1 0 - 2 2 ( \u00C2\u00A3 ' / k e V ) - 3 cm2. Again the column densities are adequate for X-ray absorption, but the absorption would be very sensitive to the inclination of the wake 86 No. 3, 1978 with respect to the X-ray star, since the far wake is very narrow. V. CONCLUSIONS The supersonic accretion of gas onto a neutron star has been described, working from the basic model as shown in Figure 1; the main features are the sheath and the accretion column. The angular width of the column, a measurable quantity in the X-ray light curve, is found to depend on the ratio of V2 to the column temperature, and therefore yields information about the local wind velocity provided the column tem-perature can be specified. An accurate estimate of the temperature would require a hydrodynamic calculation including radiation transfer, but upper limits to the temperature can be obtained by estimating the relevant heating and cooling rates. !The most important consideration in determining the thermal state of the gas is whether or not the gas can cool before it falls all the way down the accretion column. The cooling line of Figure 2 separates the flow into two main regimes. If the postshock gas in the sheath is unable to cool, it will fall inward adi-abatically in a wide accretion column, with the accretion efficiency (the /? factor) reduced by the thermal pressure. Below the cooling line of Figure 2, the gas will cool to an equilibrium value determined by the radiation field. In the region of the density and velocity parameters which apply to the stellar wind X-ray binaries, this means that the upper part of the accretion column will be photoionization heated to temperatures of order 10\u00C2\u00B0 K. The base of the column will be heated to the Compton equilibrium tempera-1049 ture, causing the pressure to rise sufficiently that the base of the column will spread to a broad inflow. It is predicted that the electron scattering will cause the X-ray light curve absorptions to change from single dips to double dips as n0V0~2 increases, if the param-eters are in the cooling region. The gas that is gravitationally perturbed but does not become bound to the neutron star forms the far wake. The high velocity and low viscosity indicate that the far wake is almost certainly turbulent. An exten-sion of the similarity description of supersonic wakes experimentally studied provides the temperature and density in the wake. But the cooling time of the gas in the wake is then found to be less than the basic turbulence time scale, which may mean that whole description is invalid. In spite of this, we suggest that the far wake is composed of a hot gas entering the wake and denser clumps of cold gas, a description which is marginally consistent with the \"wake\" observations of Conti and Cowley (1975). The model outlined here is intended to be useful for providing qualitative insight into the physics of super-sonic accretion. The numerical quantities employed are expected to be accurate to a factor of 3 or so, and should provide basic regimes which can be further explored with a numerical model. G. G. Fahlman provided invaluable advice and criticism, and read several rough drafts of this paper during the course of this research. General encourage-ment and useful discussions were provided by members of the UBC Institute of Astronomy and Space Science and the Dominion Astrophysical Observatory. RADIATIVE EFFECTS APPENDIX In the case of a heated gas, the thermal pressure forces may become large enough to destroy the assumption that the flow is dominated by inertial forces. In this Appendix the region of validity of the supersonic description of accretion is examined. The incoming free stream may be heated such that the Mach number, M = V0/{2yRT)1'2, becomes less than I. The maximum temperature which can be produced by Compton heating is 2.9 x 107(a/0.5) x (e/5 keV). This temperature can be attained for \u00C2\u00AB U F 8 \" Z > 13.8, which is shown as the crosshatched area in Figure 2. This gives only the area for which subsonic flow is guaranteed in the presence of Compton heating, but what is really wanted is a line on which the Mach number is equal to 1. If only Compton heating is considered, we find that (non-equilibrium) temperatures are produced such that the Mach number is less than 1 for ' ? n - / 8 ~ 4 > 17.2. This line (MC < 1) is shown in Figure 2. If photoionization heating is included, the subsonic region is increased very slightly at low velocities. We conclude that most of the box of Figure 3 is indeed in supersonic flow. The deviation of streamlines from particle trajec-tories will depend on the ratio of pressure forces to inertial forces. As a worst case we assume that the gas is fully Compton heated. The ratio of radial pressure force to gravitational force for gas outside the shock is 1^ dp /GM _ RT0 /GM p 8r SRT ( A l ) This implies that the net force is outward for r > 5RT, where the thermal radius is RT = GM/SRT0 = 3.19 x 10 1 07' 7 _ 1 cm. Similarly, in the transverse direction, the ratio of pressure to the momentum flux is r 5%. (A2) In the shock-heated sheath, the ratio of radial pressure force to gravitational force is 9/16, which will act only to reduce the effective mass of the gravitating object in the sheath. In the column the ratio of pressure forces to gravitation is, for a constant temperature, 3r/5RT. In general the pressure forces can be safely ignored, even in the presence of strong heating, provided that we remain in the area of validity of the supersonic flow assumption. 1050 CARLBERG R E F E R E N C E S A l l e n , C . W . 1973, Astrophysical Quantities (3d ed. ; L o n d o n : A th lone Press). A rons , J . , and Lea , S. M . 1976, Ap. J., 210, 792. Bond i , H. , and Hoy le , F. 1944, M.N.R.A.S., 104, 273. Buff, J . , and M c C r a y , R. 1974, Ap. J., 189, 147. Con t i , P., and Cowley, A . P. 1975. Ap. J., 200, 133. C o x , D. P., and Daltabuit, E. 197i, Ap. J., 167, 113. Dachs , J . 1976, Astr. Ap., 47, 19. Daltabuit , E., and Cox , D. P. 1972, Ap. J., Ill, 855. Danby , J . M . A . , and C a m m , G . L., 1957, M.N.R.A.S., 111, 50. Dav id son , K., and Ostriker, J . P. 1973, Ap. J., 179, 585. Demetriades, A . 1968, AIAA J., 6, 432. Eadie, G . , Peacock, A. , Pounds, K. A. , Watson, M . , Jackson, J . C , and Hunt , R. 1975, M.N.R.A.S. Short Comm., Ill, 35p. Felten, J . E., and Rees, M . J . 1972, Astr. Ap., 17, 226. Hatchett, S., Buff, J . , and M c C r a y , R. 1976, Ap. J., 206, 847. Hoy le , F., and Lytt leton, R. A . 1939, Proc. Cambridge Phil. Soc., 35, 405. j Hunt , R. 1971, M.N.R.A.S., 154, 141. I l larionov, A . F., and Sunyaev, R. A . 1975, Astr. Ap., 39, 185. Jackson, J . C . 1975, M.N.R.A.S., 112, 483. Jones, C , F o r m a n , W., Tananbaum, H., Schreier, E., Gursky , H., Ke l logg, E., and G iaccon i , R. 1973, Ap. J. (Letters), 181, L43. L a m b , H . 1924, Hydrodynamics (5th ed.; Cambr idge : C a m -bridge University Press). M c C a r t h y , J . F., and Kubo ta , T . 1964, AIAA J., 2, 629. M c C r a y , R., and Hatchett, S. 1975, Ap. J., 199, 196. Mestel, L. 1954, M.N.R.A.S., 114, 437. Pounds, K. A . , Cooke , B. A . , Ricketts, M . J . , Turner, M . J . , and Elvis, M . 1975, M.N.R.A.S., 172, 473. Ruderman, M . A. , and Spiegel, E. A . 1971, Ap. J., 165, 1. Spitzer, L. 1962, Physics of Fully Ionized Gases (New Y o r k : Interscience). Townsend, A . A . 1976, 77;e Structure of Turbulent Shear Flow (2d ed.; Cambr idge: Cambr idge University Press). R. G. CARLBERG: Department of Geophysics and Astronomy, University of British Columbia, Vancouver, B.C., Canada, V6T 1W5 88 APPENDIX 2. GAS PHYSICS P hot c i o ni z a t ion The simplest process to describe i s the photoionization rate which i s independent of the gas temperature and density, and i s simply given by where tf\"cj i s the photoionization cross-section of atom i , i o n i -zation l e v e l j . For these calculations the mean radiation f i e l d J was taken simply as the flux energing from a model atmosphere allowing for geometrical d i l u t i o n . The numerical computations used a model computed by Hihalas (1972), s p e c i f i c a l l y the Non-LTE 50,000 K, log g = 4 model. Mihalas gives the radiation f i e l d i n terms of the emergent f l u x , F , whereas the mean int e n s i t y 4 J i s needed for the i o n i z a t i o n and heating c a l c u l a t i o n s . To make t h i s change, the mean i n t e n s i t y was assumed to vary as the d i l u t i o n , H= 1/2 (1-(1-(r^/r) 2) l/ 2) , and the flux as ( r * / r ) 2 , where r i s the s t e l l a r radius. For the computations i l l u s t r a t e d here, a l l done at a ra d i a t i o n f i e l d corresponding to the surface of the s t a r , 4/ = 2 (fTy). Since the 0 VI ion i s of p a r t i c u l a r importance the pho-toi o n i z a t i o n rates of a l l ions up to comparable ionization po-t e n t i a l s (IP of 0 V i s 113.9eV) were included. The photoionization cross sections used in the c a l c u l a t i o n s are given below. The cross sections used are given in below. The Hydrogen cross section was taken from Bethe and Salpeter (1957). hi/ ( A l l . 1) 89 2 * r r l * AO* 77c) ^ e^p (- ^ <*> and l / j = 1^ w i t h a l l ^'s replaced by ji's. The constants are *=2. 182846 0 = 1.188914 Z r 2 Zb=1 k=( (hV-24. 587e\u00C2\u00A5)/13.59 8eV) V 2 The Helium II cross section i s the same as the p h c t c i c r i z a -t i c n cross section of Hydrogen but with Z=2 everywhere. 90 For the remaining ions the cross sections have been calcu-lated by various authors using the pri n c i p l e s of guantum me-chanics, and then making a f i t to a standard polynomial to rep-resent the data as a function of incident photon energy. Two forms for the polynomial are used here, one due to Seatcn (1958) ( A l l . 4) where i ^ , i s the threshold freguency and the other constants are f i t t i n g parameters. The other form i s a s l i g h t l y more general polynomial due to Chapman and Henry (1971) *\u00E2\u0080\u00A2<.\u00C2\u00BB)* 2 g 0(n) 1. 1330 1.0785 0. 9935 + . 25 28n~ , 1296n~^ g, (JJ) -0.4059 -0.2319 -n~\u00C2\u00BB {.6282-. 5598n~*\u00E2\u0080\u00A2 . 52 99n-*) g^(n) .07014 . 02947 n~* (.3887-1. 181n-*+1.470n-z) The values i n t h i s table were taken from Johnson (1972). In order to obtain the t o t a l recombination rate, recombina-tions to the l e v e l s n=1 to 9 summed together. Computations of the recombination rate f o r a l l of the other ions cf i n t e r e s t have been made by Aldrovandi and Peguignot 93 (1973), with errata i n Aldrovandi and Peguignot (1976). The data i s provided in the form of f i t s to simple functions. The radia-tive rate i s given by *r= *r\u00C2\u00BBel A**-* (n,n\u00C2\u00AB\u00C2\u00BB) + h ^ * A**-i (n,n 1\u00C2\u00BB\u00E2\u0080\u00A2) +h v t where i n general n* = n+1 and n'*>>n, n I f the gas becomes s u f f i c i e n t l y dense or i f the radiation f i e l d strong enough, the electron i n the high l y i n g guantum l e v e l n*1 can be either c c l -l i s i c n a l l y or r a d i a t i v e l y ionized out of the atom before i t has time to s t a b i l i z e by photoemission. A rough empirical correction factcr was devised to allow f o r t h i s decrease i n the di e l e c -tronic recombination rate. The p r i n c i p a l guantum number of the state at which half the captured electrons are s t a b i l i z e d by cascades to ground and half are reionized i s given by, 1(collisions)=<1.4x10 l s2 6T l/ 2/ne) l/ 7 Dupree (1968) 1 (radiative) -Z (3 flUn (1) / (HkT r < w j )) V 2 Sunyaev and Vainstein (1968) where W i s the geometrical d i l u t i o n factor of the radiation f i e l d approximated by a blackbcdy of temperature Trad* These numbers can be calibrated against the depression of the recom-bination rate calculated by Summers (1974). I t was found that data i s roughly f i t t e d by the multip l i c a t i v e factor f, such that ^=t *^(n=0,W=0), where f i s f = exp[-2. 303* (. C15*a2+,092*a) ] where a = 12.55-7*log10(1). 96 That i s , the adjusted recombination rate i s found by multiplying the value found from the f i t s given by Aldrovandi and Peguignot times the f factor given above. In addition to t h i s correction to the d i e l e c t r o n i c r a t e , the semicoronal approximation of Wilson (1962) has been used to add tc the radi a t i v e rate. This allows for some recombinaticn to upper l e v e l s , = 1.8X10-A\u00C2\u00BB \"&j(kT)-a/* Y c j ( A l l . 10) where Xcjl 1) ~ ^ / l 2 ( c o l l i s i o n s ) . In addition three body recom-bination makes s i g n i f i c a n t contributions at low temperatures, and i s simply approximated by (Burgess and Summers 1976) ^^=1,16xl0-\u00C2\u00ABJ3fV 2 f l e , (All.11) where J i s the charge of the ion. , C o l l i s i o n a l Ionization The rate of c o l l i s i o n a l i o nization for Hydrogen was also taker from Johnson (1972) as (All.12) where ^ Z f f t r \" t o i*3 V * vL = %I^/kT z' v = x 0 (I^/kT+r^) | (t) = E 0 (t)-2E, (t)*E^(t) , 97 and the Gaunt factors and x are as for the recombination rate in Hydrogen. Only i o n i z a t i o n from the ground state n=1 w i l l be considered, so r|=0.45 and b(= -0.603. A l l other ions have c c l l i s i o n a l ionization rates based upon an approximation investigated i n d e t a i l by l o t z (1967). A s l i g h t modification to the o r i g i n a l formula has been made by McWhirtier (1975) to allow for the decrease of the io n i z a t i o n rate at high temperatures. The rate i s given by where s goes from 1 to 2 in the calculations here, n (s) i s the number of electrons i n the subshell, and %\u00E2\u0080\u00A2\u00E2\u0080\u00A2{ 1) i s the normal J ionization potential as given by Allen (1973), yt'(2) i s \")(t''(1) plus the excitation energy of the lowest excited l e v e l i n the new ion with one of the inner s h e l l electrons removed. For i n -stance, the ionization of C II which has an e l e c t r o n i c c o n f i -guration of 1s 22s 22p can proceed with the addition cf 1S6659 cm - 1 of energy to C I I I 1s 22s 2 by removing the one outer s h e l l electron, or the i o n i z a t i o n can take 196659*52315 car - 1 and ionize to C I I I 1s 22s2p, by removing one of the two s s h e l l ele-ctrons. The energy 52315 cm-* i s simply the energy to go from C III 2 s 2 to C III 2s2p. The attached table gives ionization po-t e n t i a l s (in eV) and the number of subshell electrons. , The values were obtained from the tables of l o t z (1967) and Moore (1949) * TABLE 9: IONIZATION POTENTIAL AND SHELL ELECTBON POPULATIONS ATOM IP1 N1 IP2 N2 H I 13.598 1 HE I 2 4.587 2 z. (All.13) 98 HE II 54.416 1 C I 11.260 2 16.6 2 C II 24.383 1 30.9 2 c II I 47.887 2 323. 2 c IV 64.492 1 342. 2 c V 392.08 2 c VI 489.98 1 N I 14.5 34 3 20.3 2 N II 29.601 2 36.7 2 N II I 47.448 1 55.8 2 N IV 77.472 2 469. 2 N V S7.89 1 492. 2 N VI 552.06 2 N VII 667.03 1 0 I 13.618 4 28.5 2 0 II 35. 117 3 42. 6 2 0 III 54.934 2 63.8 2 0 IV 77.413 1 87.6 2 0 V 113.90 2 642. 2 0 VI 138.12 1 66 9. 2 0 VII 739.32 2 0 VIII 871.39 1 NE I 21. 564 6 48.5 2 NE II 40.962 5 66.4 2 NE III 63.45 4 86. 2 2 NE IV 97.11 3 108. 2 NE V 126.21 2 139. 2 NE VI 157.93 1 172. 2 NE VII 207.26 21072. 2 ME VIII 239.09 11106. 2 NE IX 1195.8 2 NE X 1362.2 1 MG I 7.646 2 60.420 6 MG II 15.035 1 67.809 6 MG I I I 80.143 6 118.768 2 MG IV 109.31 C -/ 14 4.42 2 MG V 141.27 4 17 2. 0 1 2 MG VI 186.51 3 201.22 2 MG VII 224. 95 2 241.14 2 MG VIII 265.92 1 283.38 2 MG IX 328.0 21680.4 2 MG X 367.5 11719.8 2 MG XI 176 1.8 2 MG XII 196 3. 1 SI I 8. 151 2 13.616 2 SI II 16.345 1 22.870 2 SI I I I 33.492 2 137,709 6 SI IV 45. 141 1 149. 358 6 SI V 166.77 6 217. 170 2 SI VI 20 5.08 5 250.48 2 SI VII 246.49 4 285.26 2 SI VIII 303. 16 3 321.76 2 SI IX 351. 1 2 371. 2 2 SI X 401. 4 1 422.4 2 SI XI 476. 1 22340.8 2 SI XII 523. 12388. 2 99 SI XIII 2438. 2 SI XIV 2673. 1 S I 10.360 4 20.204 2 s II 23. 33 3 33.747 2 s I I I 34.83 2 43.737 2 s IV 47.30 1 57.60 2 S V 72.68 2 24 3.31 6 s VI 68.05 1 258.68 6 s VII 280.01 6 342.45 2 s VIII 3 28.33 5 352.24 2 s IX 379. 1 4 402.8 2 s X 447. 1 3 46S.0 2 S XI 504.7 2 551. 2 2 s XII 565. 1 621. 2 s XIII 652. 2 s XIV 707. 1 s XV 3224. , 2 s XVI 3494. 1 Again, following Silson (1962) we make a small addition to the ionization rate allowing for high density e f f e c t s cf i o n i z a -tions out of upper l e v e l s , = 4.8x10-\u00C2\u00AB I-*/ 2exp {-Vy/kTJ/OCiyiz { c o l l i s i o n s ) ) . ( A l l . 14) Charge Exchange In order to increase the general usefulness of t h i s program the charge exchange rates of H+ \u00E2\u0080\u00A2 0 ^ 0+ + H H+ + N ^ N+ + H were included using expressions exactly as given by F i e l d and Steigman (1971) and Steigman, et a l . (1971). Since the tempera-tures here are usually i n excess of 10* K, the charge exchange rate i s at almost constant and at i t s maximum. , 100 The Heating Bate A l l heating i s due tc energy gain by photoionization, which is simply given by (All.15) The t o t a l gain i s \"J J ( A l l . 16) Cjgcl in.g Bates The emission of radiation i s calculated under the assump-tion the medium i s o p t i c a l l y t h i n . The cooling due to bremsstra-hlung i s (Cox and Tucker 1969). , - A 8 = 2.29x10-2* 1-1/2 n ^ / n 2 ( A l l . 17) where n\u00C2\u00BB i s the number density of Hydrogen. This loss mechanism dominates for temperatures i n excess of 10 7 K. The radiative recombination energy loss rate i s -A'Jj = * rCj ( V,;j., + kT) X t-j A c r- o . 0 7 / 3 + \u00C2\u00B1 J U V *\u00E2\u0080\u00A2 O.bH Lt~''* 1 ( A l l . 18) where 0=/X^/kT. The correction factor i n brackets was derived frcm the analysis of Seatcn (19 59) for the recombination process in Hydrogen. I t represents the correction to the radiative re-combination rate reguired to convert i t to the energy rate, ac-counting for the pref e r e n t i a l capture of slow electrons. The l o s s rate due to die l e c t r o n i c recombinations was e s t i -1 0 1 mated as, - A ^ =^7}<{^ yj H+ 6 E ^ ) X t J - At- . (All.19) The recombination radiation i s the dominant loss mechanism for temperatures of 2x10* K and l e s s . .The energy difference AErj i s taken as the lowest energy permitted t r a n s i t i o n tc the ground state. Between 2x10* and 10 7 K the dominant loss mechanism i s c c l -l i s i c n a l excitation of l i n e s . In princip l e a c a l c u l a t i o n of t h i s rate reguires the c o l l i s i o n a l cross section for excitation of a particular t r a n s i t i o n as a function of incident electron energy. With the aid of the Milne r e l a t i o n , which relates the c o l l i s i o n cross section tc the inverse process of photcabscrpticn, the loss rate can be approximated as (Mewe 1972) -A* = 1.7x10 - 3 T ~ V 2 f. \" I g U r f t] ( U e x p i-Altj (l)/kT)A c X(J' ( A l l . 20) where g i s a gaunt factor, ftJ (1) i s the f value for the t r a n s i -t i o n , and 4Et'j (1) i s the energy of the emitted photon. The g factor has been calculated by Mewe (1972) for many tr a n s i t i o n s and given a simple extension by Kate (1976) to cover a l l t r a n s i -tions. They both use the same f i t t i n g function for the Gaunt facto r , g = A* (By-Cy2+D) exp (y) E, 2 BY T C HIS NUMEEB C FE: GUESS AT ELECTRON DENSITY C WF: DILUTION FACTOR FOB B AEI AT ION FIELD C NLINE: NUMBER OF LINES IN CCCIING CALCUL C ATION C WLINE: WRITE INDIVIDUAL LINE CCCIING AND HEATING C TCL: TOLERANCE FCB CONVERGENCE OF BNOT,DENE,EQUILIBRIUM C TEMPERATURE C EGUIM: TRUE FOB FOBCING BALANCE OF HEATING AND COOLING C BATES C CHABGX TBUE FOB CHABGI EXCHANGE H-N, H-0 CALCULATIONS C NELMNT: # OF ELEMENTS STARTING WITH H IN IONIZATION CA C LCULATICN C USEFUL SOMETIMES IN EQUIM FOB PBELIMINAEY C ESTIMATE C DMAX: MAX FRACTIONAL CHANGE ALLOWED IN DELTA T, C PBEVENTS WILD OSCILLATIONS C DIELEC: FALSE TURNS ALL DIELECTBONIC RECOMBINATION OFF C MI-OOP: NUMEEB OF ITEBATIONS ALLOWED IN CONVEBGENCE TO C T IF EQOIM IS CN C TBAD: BACIATION TEMP EBATUBE OE PHOTON SOURCE C WFTBAD IS APPROXIMATELY A BRIGHTNESS TEMPEBATUBE C FUDGE: TBUE FOB BE DUCT ION OF DIELECTBONIC B ECO MB IN A HO N C WITH DENSITY C FUDGE FACTOR CALCULATION IS DESIBED C THEEEB: THREE BODY RECOMBINATION C DXENDT TRUE FOB COMPUTING EEBIVATIVES IN N AND T C SEBIES IS (T,N) , (T* { 1+-DERDEL) , N) ) . 0 AWAY FROM STAR, SIMILAB1Y Z INCREASES UPWARDS C V0=0. DV-0. NCNEQ=. TRUE. DTOL=.05 EYNEQ=. FALSE. , DYNIT-10 DYNTOL=1.E-3 SLIM=.5D0 TRAD=50000. CONDUC=.TRUE. DNON=1. E11 TNON=0. T= 1. E5 DRHO=0. DT=0. D2T=0. NIINE=874 CHF=1. SLAE=0. VERBOS=.TBUE. WLINE=. FALSE. SKIP=.FALSE. SPHEBE=. FALSE. STARMU=1. EST AR= 1. E12 158 BEHIND 1 BEHIND 2 C C BEAD IN LINE AND FLUX IAT A (USUALLY FILE L FOB C ED AT) C BEAD (2) ELINE,F,II,JJ,FLUX,FPT C C INCOMING FLUX HAS BEEN MULTIPLIED BY 4*PI/C C DG 49 18=1,96 DEL E (IN) =FPT (IN + 1) -FPT (IN) 49 CONTINUE BEHIND 3 BEAD(3) PHOT,PHEAT,SIGMA,FTOT1 BEHIND 4 BEAD (4,9774) (NDO (IJ) , NUF< IJ) , I J = 1, 76) 9774 FOBMAT (214) C C BEADS FILE CONTAINING INTEGBATTON FBEQUENCIES WHICH HAS C USED BY PHOTION C 998 BEAD (5,PABAM, END=999) IF(.NOT.VEBBOS) WLINE=.FALSE. IF {. NOT. NONEQ. AND. VO. EQ.O.) V0=1. D7 WBITE(6,1014) C C BEAD IN SET OF 5 OUTPUT QUANTITIES FROM HEAT COOL PBGGBAM C 1014 FOBMAT(*1\u00C2\u00BB) WBITE(6,PABAM) IF (SKIP) GO TO 111 DC 14 1=1,5 BEAD(1) D,T A,W FJ,FJ,GAIN,COOL,LBB EMS,LBRAD,LINLOS, $ EINT,ENTH,DE,ABUND,X,LCLX,CPHEAT,FTOT AD(I) = D AT (I) =TA AHFJ (I) = WFJ AGAIN(I)=GAIN AC0CL(I)=CO0L AE(I) =EINT AH (I) = ENTH ADE (I) =DE DO 15 IZ=1,9 IQ1=IZED (IZ) + 1 DO 15 J=1,IQ1 AX(J,IZ,I)=X(J,IZ) 15 CONTINUE 14 CONTINUE FA=ABUND(3)/3.0 35918E-04 WBITE (6,1000) FJ,FA,FTOT 1000 FOBMAT('OFJ,IA,FTOT*,2F11.6,E15.7,/* DEN,TEMP,WFJ,H* $ ,*EAT,COOL,\u00C2\u00AB, $ * EI NT, ENTHALPY, DEN E* ) DO 1 1=1,5 AGAIN (I) = AGAIN (I)/AD(I) 1 CONTINUE WBITE (6,1001) (SD (ID) , AT (ID) , AW FJ (ID) ,AGAIN (ID) , $ ACOOL (ID) ,AE(ID) , $ AH (I D) , ADE (IE) , ID = 1 , 5 ) 1001 FOBMAT(1X,8E15.1) C C SET UP CENTRAL VALOES C DEN=AD (1) DENE=ADE(1) T= AT (1) N0=DEN N10=DENE TG=T DEINT=AE(1) DENTH=AH (1) INTH= AH (1) ST4=SQRT(SNGL(T)*1.E-04) C C THERMAL VELOCITY ALREADY DIVIDED BY C C DO 11 1=1,9 VTHERM (I) = 4. 28 33E-5/SQRT (AM ASS (I) ) *ST4 11 CONTINUE DEN 1= ADE (1) C C CALCULATE UPPER AND LOWER DERIVATIVES DIFFERENCES C DNU=AD (4)-AE (1) DNL=AD(1)-AD(5) DTU= AT (2) - AT (1) DTL= AT (1) - AT (3) C C DERIVATIVES OF THE IONIZATION FRACTIONS C DO 16 IZ=1,9 IQ=IZED(IZ) DO 17 J=1,IQ GRIJ (J,IZ)=0. DXDT(J,IZ)=0. EXDN (J,IZ)=0. TF(AX(J,IZ, 1) .LE. 1.E-10) GO TO 17 DXDT ( J , IZ) = \u00E2\u0080\u00A2 5* ( (AX (J , IZ , 2) - AX (J , IZ, 1) )/DTU + $ (AX (J,IZ,1)-AX(\u00C2\u00AB3,IZ,3) )/DTL) IF ( E AES ( EX DT (J, IZ) *T/AX (J ,IZ,1) ). LE. DTOL) DXDT (J,IZ) $ =0.D0 DXDN(J,IZ) = . 5*( (AX (0,IZ,4) - AX (0 ,IZ, 1) )/DNU + $ (AX (J, IZ, 1) - AX (J , IZ , 5) )/DNL) IF (DABS(DXDN (J,IZ)*DEN/AX (J ,IZ,1)).LE.DTOL) DXDN (J, $ IZ)=0. DO 17 CONTINUE 16 CONTINUE IF ( .NOT.VERBOS) GO TO 20 DC 19 IZ=1,9 IQ=IZED(IZ) WRITE (6,1029) IQ 1029 FORMAT(* DXDT,DXDN FOB ATGM Z=*,I2) WBITE (6,1030) (DXDT(J,IZ) ,J=1,IQ) 1030 FOBMAT(1X,10E13.5) WBITE (6,1030) (0XDN(J,IZ) ,J=1,IQ) 19 CONTINUE 20 DO 18 IZ=1,9 IQ=IZED(IZ) DO 18 J=1,IQ X(J,IZ)=AX(J,IZ,1) 18 CONTINUE C C PHYSICAL DEBIV ATIV ES C DU=0. DL=0. CALL DERIV(ADE,DNEIN,DNEDT,DTOL) WBITE (6,1003) DNEDN,DNEDT 1003 FOBMAT {* DNEDN, DNEDT* , 2E1 5. 7) DN2=AD <1)*AD (1) CALX DEBIV(ACCOL,DLDN,DLDI,DTOL) DLDN= AD (1) *2.*ACOOL (1) +DL DN *DN2 DLDT=DLDT*DN2 BL=ACCCL (1) *DN2 WRITE(6,1005) DLDN,DLDT , EL 1005 FOBMAT ( * DLDN,DLDT,LOSS ES *,3E15.7) CALL DEBIV (AGAIN,DGDN,DGDT,DTOL) GAIN= AGAIN (1)*DN2 DGAINC=GAIN/2.9979E10 DGDN=AD (1) *2.*AGAIN (1) +DGDN*DN2 DGDT= DGDT* DN2 WRITE (6,1007) DGDN,DGDT,GAIN 1007 FORM AT ( * DGDN, DGDT , GAINS * , 4E15. 7) CALL DERIV(AE,DEDN,DEDT,DTGL) C= D EDT/KBOLTZ WRITE (6,1009) DEDN,DEDT,C 1009 FOB MAT (* DEDN, DEDT , CV \u00C2\u00BB , 4E15. 7) CALL DERIV (AH,DHDN,DHDT ,DTOI) C=DHDT/KBOLTZ C C CCNDUCTION CALCULATION FBGM SPITZEB C WRITE (6,1011) DHDN, DHDT, C 1011 FORMAT(* DHDN,DHDT,CP\u00C2\u00BB,4E15.7) IE (CCNDUC) GO TO 110 DKDT=0. DKDN=0. CONKAP=0. GO TO 113 110 CONTINUE COULOG = 9.00 + 3. 45*ALOG10(SNGI (T) )-1. 15*ALOG 10 ( DENE) CGNKAP=1.8E-5*T**2.5/COULOG DKDT=2.5*CONKAP/T+CCNKAP/ (COHLOG*T) *3.45 DKD N=CCNKAP/(COULOG*DEN\u00C2\u00A3) * (-1 . 15) C C THE MEAN MASS OF AN ATOM COEFFICIENT C 113 WRITE (6, 1013) CCNKAP, DKDT, DKDN 161 1013 FORMAT(* CONKAP,DKDT,DKDN\u00C2\u00BB,3E15.7) EHCCCN=0. BO 10 1=1,9 BHGCGN=RHOCGN*AMASS (I) *AEUND(I) 10 CONTINUE C C TEE FORCE DUE TO CONTINUUM RADIATION C RHOCV=RBOCON*1.660531E-2U 111 CONTINUE GRADC=0. DGCDN=0. DGCDT=0. DO 51 1=1,9 IZ=IZED(I) DO 52 J=1,IZ IJ=INDEX(J,I) XAC=ABUND(I) *X (J,I) IF(XAC.LT. 1.E-10) GO TO 52 KUB=NU0(IJ) NUE=NUF (IJ) IF(NUE. EQ\u00C2\u00AB NUB) GO TO 52 GCL=0. DO 50 IN=NUB,NUE GCL=GCL+.5* (FLUX(IN + 1) *SI GM A (IN \u00E2\u0080\u00A2 1, IJ) +FLUX (IN)* $ SIGMA ( I N , I J ) ) * $ DELE (IN) 50 CONTINUE GR ADC=GR ADC+GCL*X AC DGCDN=DGCDN*GCL*DXDNMTMCN* XNO*ECT+MTMKN*DN*ECT CRD (3,1) = $ VG**2*ECT+VG* (DVG*EKT-MTMCT*EKV-MTMKT*ECV+EKT*DV) $ +ECN*MTMKT* XN0 + EKN*H1MCT*N0+EKN*MBKT*DK + DVG* (-CONKAP) *DV-MTMKT $ *EKV*DV $ -BTMCN*DN* X(-CCNKAP) -MTMCH*N0*EKT-MTMKN*DN*EKT-HTMKN*NO*ECT CID (3, 1)=0. DO CRD (4, 1) =O.DO CID (4,1) = I VG**2*EKT+VG* (DVG* (-CONKAP) -MTMKT*EKV* (-CGNKAP) *DV $ ) + $ EKN*HTMKT*NO-XHTMCN*N0*(-CONKAP) -MTMKN*DN* (-CGNKAP) -MTMKH*N0* EKT CED (5,1) = $ -VG**2*(-CONKAP)+ MTMO*N0* (-CCNKAP) CID (5,1)=0.DO CRD (1,2)=0.D0 CID(1,2) = $ EWN *MTMCT*DN+D VG*ECT*D VG*E KT*DV-MTMCT*ECV-MTMCT* $ EWV*DV-\u00C2\u00ABIMCN* XDN*E8T*ECT*DV CRD (2,2) = $ VG* (-DVG*E8T+MTMCT*EWV-2*ECT-EHT*DV)-E8N*MTMCT*NO $ -EWN*MTMKT* XDN-DVG* EKT+MTMCT*EKV\u00E2\u0080\u00A2MTMKT*\u00C2\u00A3CV+MTMKT* EH V*D V+MTMCN* -$ N0*E8T + MTMO*DN X*EHT-EKT*DV CID (2,2)=0,D0 CRD (3,2) =0. DO CID (3,2) = $ -VG**2*EHT+VG* (MTMKT'*EWV-2*EKT) -EWN*MTMKT*NO-DVG* $ (-CGNKAP)+MTHKT X*EKV+MTMKN *NO*EBT-(-CONKAP)*DV CRD (4,2) = $ 2*VG*(-C0NKAP) CID (4,2) =0. DO CRD (5,2)=0.D0 CID <5,2)=0. DO CRD (1,3) = $ DVG*EwT-MTMCT*EwV4ECT+EHT*DV CID(1,3)=0.D0 CRD (2,3)=0. DO CID (2,3) = $ 2* VG* EHT- MT H KT * EH V \u00E2\u0080\u00A2 EKT CRD (3,3) = $ - (-COMKAP) CID (3,3)=0.D0 CRD(4,3) =0.D0 CID(4,3)=0.DO CRD(5,3)=0.D0 170 CID (5,3) = 0. DO CBD(1,4)=0.D0 CID <1,4) = $ -EST CBD (2,4) =0. DO CID(2,4)=0.D0 CBD (3,4)=0. DO CID (3,4) =0.D0 CBD <4,4)=0. DO CID(4,4)=0.D0 CBD (5,4)=0. DO CID(5,4)=0.D0 WBITE (6,1000) ((CBD (J.I) ,J=1,5) , (CID (J,I) , J= 1,5) ,1= $ 1,4) 1000 FOB H AT (' 0 * , 5D25.10 , /, 1 X, 5 D2 5.10) EETUBN 999 EETOBN1 END PROGRAM DISPER C C A EEOGBAM TO FIND THE BOOTS OF THE CUBIC DISPERSION C RELATION C FCUKD FOE THE CASE OF ONE DIMENSIONAL PLANE HAVES C SUBROUTINE COCALC USES THE OUTPUT FBOM COEF TO C ACTUAL FIND THE POLYNOMIAL C C C THE PARAMETERS CONTROLLING THE EEOGBAM ARE: C KLGG IF TBUE A LOGARITHMIC SERIES OF K ARE GENERATED C KMIN MINIMUM K VALUE (IF KLOG THEN THIS IS A LOG) C KM AX MAXIMUM K VALUE C PLREAL PLOT REAL FREQUENCIES C ELIMAG PLOT IMAGINARY IBEQUFNCIIS C KINC INCREMENT BETHEEN K VALUES C TGI IF TWO BOOTS L I E CLOSES THAN THIS THEY ABE C ASSUMED TO BE IDENTICAL C EBB ACCUBACY OF NDINVT SOL'N TO B=0, DD=0 C MAXIT # OF ITERATIONS IN NDINVT C EQTCL ACCURACY OF NEWTON ITERATION ^ BOOT IMPROVER C NFILE FILE LINENUMBEB WHERE COEFFICIENTS START C USUAL 1+A MULTIPLE OF 5 C LABEL TRUE TO LABEL PLOTS C SUEDIV # OF SUBDIV USED EY NEXT BOOT ESTIMATOB C PEMIN FOB DIFFEREN BEAL PLOT Y AXIS MINIMUM C PBINC SAME EUT INCREMENT FROM MIN FOR 10\" PLOT C PIMIN SAME AS PRMIN BUT FOR IMAGINARY PART C PIINC SAME AS PBINC EUT IMAG C NPRINT EVERY NPRINT'TH K VALUE AND ROOT OUTPUT TO PRINTER C ER C SEMIV IF TRUE OUTPUT W i l l ALLOW NPRINT TO TAKE EFFECT C REAL*8 DKR (210) REAL* 8 BD( 11), ID (11), BOOTS (4) ,ROOTI(4) ,PKR(5) ,PKI(5 $ ) EEAL*8 EDIS(12) ,IDIS (12) SEAL*8 DELTAK REAL*8 DIST,DISTOL ,RSAVE,DERR INTEGER*4 KEEG (3) , KEND ( 3) REAL*8 X(4) ,F(4),ACCEST (4) ,EEB EEAL*8 BBC (4) ,IDC(4) SEAL*8 DONE /1.D0/ EEAL*8 WIM (3,210) ,WB(3,210) BEAL*8 KMIN,KMAX,KINC,TOL BEAI*8 DSIGN, DEE AL, DIM AG, DLOG10, DABS, DMIN1 ,DMAX1 COMPLEX*16 DCMPLX,CECOT,CDLCG REAL*8 CDABS,BLOG LOGICAL FREEZE,NONEQ,KLOG,PLREAL,PLIMAG,SOLVEQ,KNEG $ ,MANY,RESTOR,FIRST,DUPLIC,LABEL LOGICAL*1 INSTAB(3,210),ALABEL(80) I NT EG EE* 4 SYM(3) /12,2,5/ REAL*8 GEEL,EQTOL,COMPAE,LDEI,DPS,DPI REAI*4 VPHASE(3, 210) ,VG(3,210) , AX (210) ,AY(210) C0MPLEX*16 CK,DC (4),DDC (4) ,8C,\u00C2\u00BBC2 ,WC3,NDC(3) ,DDIS, $ DDDK C0MPLEX*16 DwDK,PRED (3) ,SLOPE (3) LOGICAL V EE EOS, ZEE IMG, FANCY, SEMIV,LPRINT INTEGER*4 SUBDIV INTEGER*4 INSTBI (3, 40) , NMAXL (3) REAL*8 BMAXL (3) ,DAWIM CGHHCN /NEWT/ DDIS,MAXIT COMMON /CCCALC/ CDC,EC,NEC COMMON /CPOL/ RDC,IDC COMMON /CONTEO/ SOLVEQ COMMON /CNTB02/ MANY,BESTOR COMMON /DIS/ RDIS,IDIS NAMELIST /PARAM/ KLGG,KMIN,KMAX,PIREAL,PLIMAG,KlNC, $ TOL $ ,ERR,MAXIT,VERBOS,FANCY,PINC,PKMIN,MANY $ , RES TOR, EQTOL, NFILE,LABEL,S OBDIV $ ,PRMIN,PBINC,PIMIN,PIINC,NPRINT,SEMIV EXTERNAL FCN C C SET UP DEFAULT VALUES C SUEBIV=3 NFILE=1 TCL=1. D-6 EQTOL=1.D-14 FLIMAG=. TBUE. PLBEAL=.TBUE. , KLGG=. TRUE. , KMIN=-12.D0 KMAX=-2. DO KINC=.1D0 VEBBCS=. FALSE. , NPBINT=5 SEMIV=.TEUE. FANCY=.FALSE. ZEBIMG= .TBUE. , MANY=.FALSE. RESTOB=. FALSE. MAXIT-100 EEB=1,D-15 FIBST-.TBUE. LABFL=. FALSE. 99S9 IF(MANY.AND..NOT.FIBST) GO TO 9990 NFIIE=1 READ (5,PARAM,END=9S98) C C READ IN PARAMETER LIST OF WHAT TO DO C NFILE 1=NFILE- 1 IF(NFILEI.LE.O) GO TO 5 DO 6 ISKIP=1,NFILE1 READ (1) 6 CONTINUE 5 CCNTINUE IF(MANY.AND..NOT.FIRST) GO TO 9990 NK= (KMAX-KMIN)/KINC+1. 5 C C CALCULATE ABBAY OF K VALUES C IF (KLOG) GO TO 2 DO 1 I=1,NK DKB (I) = KHIN*(I-1) *KINC 1 CONTINUE GO TO 3 2 NKB= 1 SKE=NK DO U I=NKB,NKE DKB (I) = 10. DO** (KMIN \u00E2\u0080\u00A2 (I-NK B) *KINC) <4 CONTINUE 3 CONTINUE 99S0 WRITE (6,1000) 1000 FORMAT(* 1') IF (.NOT.LABEL) GO TO 9 BEAD(5,1066,END=9998) ALABEL WBITE (6,1067) ALABEL 1067 FORMAT(1X,80A1) 1066 FOBMAT(80 A1) 9 CONTINUE C C BLOT SCALING QUANTITIES C PINC=0. FKMIN=0, PRMIN=0. PBINC=0. PIMIN=0. EIINC=0. WRITE (6,PARAM) EIRSI=. FALSE. RRMN=9.E70 C C SEE IF ROOTS ARE PROPERLY SEQUENCED C RRMX=-9.E70 EIMN=9.E70 RIMX=-9.E70 DO 222 IN=1,3 VPHASE (IN, 1)=0. WMAXL (IN) = -9. D70 NHAXL(IN)=0 VG (IN, NK) =0. 222 CONTINUE CALL COCALC (69998) C C CALCULATE POLYNOMIAL COEFFICIENTS C GREL=0. DO 1=1 SOLVEQ=. FALSE. DDPLIC=.FALSE. C 17a C FIND FIRST BOOT OF POLYNOMIAL C CALL DISPCO(DKB(I)) CALL CPOLY1 (BBC, IDC,3,BOOTB,BOOTI,6999) 185 SOL VEQ=. TRUE. 99 DO 97 IN=1,3 CALL NEWTON (ROOTR (IN) ,ROOTI (IN) ,EQTOL) HR (IN,I) = BGGTR(IN) WIM(IN,I)=ROOTI (IN) 97 CONTINUE 181 DO 183 IN=1,2 IR=IN+1 184 IF(IR.GT. 3) GO TO 183 WC=DCMPLX (WR (IN, I) ,HIM (IN, I) ) HC2=DCMPLX(HR (IR,I) ,HIM (Ifi,I) ) C0MEAB=DMAX1 (CDABS (HC) , CDABS (HC2) ) IF(CDABS(HC-HC2)/COMPAR.LT.TCI) GO TO 170 IR=Ifi+1 GO TO 184 183 CONTINUE BUPLIC=.FALSE. 189 1=1+1 C C ESTIMATE THE VALUES OF THE NEXT SET OF ROOTS C IF (I.GT.NK) GO TO 98 DELTAK=DKR(I)-DKR(1-1) DO 96 IN=1,3 CALL ADVANC (ROOTR (IN) ,ROOTT (IN) , DELTAK , DHDK ,SUBDIV) VG (IN,I-1)=SNGL (DBEAL (DHDK) ) PEED (IN)=DCMPLX (ROOTR (IN) ,ROOTI (IN) ) 96 CONTINUE CALL DISPCO (DKR (I) ) C C CAICULATE COEFFICIENTS FOB NEXT K C GO TO 99 170 IF(DUPLIC) GO TO 188 I F (I. EQ. 1) GO TO 188 SOLVEQ=.FALSE. C A LI DISPCO (DKB (I) ) C C IF DUPLICATE ROOTS ABE FOUND C GO BACK TO POLYNOMIAL ROOT FINDER C TO SEE IF ANOTHER ROOT CAN BE FOUND C IF SO TRY TO PROPERLY ORDER ROOTS C HBITE(6,1099) I,IN,IB 1099 FOBMAT(\u00C2\u00BB SS6&S6SfiS8fiS6DUPLICATE ROOTS I,IN,IR\u00C2\u00AB,3I5 $ ) CALL CPOLY1 (RDC,IDC,3, ROOTR,ROOTI,\u00C2\u00A3999) NFIL=1 198 DISTOL=9.D70 DEBB=1.D70 DO 174 IN=1,3 IF (BOOTB (IN) . EQ.O. DO) GOTO 173 175 DEBB=DMIN1 (DABS (ROOTS (IN) ) , DEER) 173 I F (EOOTI (IN) . EQ.O. DO) GOTO 174 DERH=DMIN1 (DABS(BOOT! (IN)) ,DERR) 174 CONTINUE BEEE=DERR * 1.D-3 EPR=DLOG (CAES (DREAL (PRED(NFIL) ) ) +DERR) DPI=DLOG(DABS(DIMAG(PEED(NEIL))) + DERR) 194 DO 195 IN=NFIL,3 DIST-DABS ( (DREAL (PEED (NFIL) ) - BOOTB (IN) ) * $ (DPH-DLOG (DABS (ROOTB (IN) ) +DERR) ) ) $ +DABS ( (DIHAG (PEED (NFIL) )-ROOTI (IN) ) * $ (DPI-DLOG (DABS (ROOTI (IN ) ) +DER B) ) ) IF (DIST.GT,DISTCL) GO TO 1S5 DISTOL=DIST INFIL=IN 195 CONTINUE 196 WR (NFIL,I) = BOOTR (INFIL) WIM(NFIL,I)=ROOTI(INFIL) IF (INFI1.EQ.NFIL) GO TO 192 ROOTB(INFIL)=ROOTR (NFIL) BOOTR(NFIL)=WR(NFIL,I) ROOTI(INFIL)=ROOTI(NFIL) ROOTI(NFIL)=WIM(NFIL,I) 192 NFIL=NFIL+1 IF (NFIL. LT. 3) GO TO 198 GO TO 185 188 DUPLIC=. FALSE. GO TO 189 98 SOIVEQ=. FALSE. DO 100 1=1,NK LPRINT=. FALSE. IF (VERBOS) LP RINT=,TR U E. IF ( (\u00E2\u0080\u00A2 NOT. VERBOS. AN E. SEMI?) .AND. $ (MOD(I-1,NPRINT).EQ.0)) LPEINT=.TRUE. 166 DG 100 IN=1,3 C C CHECK FOB INSTABILITIES C IF(WR (IN,I) . EQ. 0. DO) GOTO 102 GR EL= WIM(IN,I)/WR (IB , I) IF (WIM (IN,I) . LT.O. DO) GOTO 102 INSTAB(IN,I)=.TBUE. IF (I. EQ. 1) GO TO 105 DAWIM=WIM (IN,I) I F (DAWIM.LT.O. DO) WMAXL (IN) = 0. DO IF(WMAXL(IN) .GT.DAWIM) GO TO 101 IF (INSTBI (IN, NMAXL (IN) ) .EQ.I-1) GO TO 104 105 NMAXL (IN) =NMAXL (IN) +1 104 INSTBI (IN, NMAXL (IN) )=I WM AXL (IN) =DAWIM GO TC 101 102 INSTAB(IN,I)=.FALSE. 101 CCNTINUE IF(DKB(I).EQ.O.DO) GO TO 288 C C FIND PHASE AND GBOUP VELOCITIES 176 C CHECK FOB ICCAL MAXIMA IF ONSIAELE C VPHSSE (IN,I)=SNGL (WE (IN , I)/IKE (I) ) 288 IF (LP.BINT) WBITE (6,1033) I , DKR (I) , WR (I N,I) ,WIM (IN, $ I) $ ,VPHASE (IN,I) ,GREI,VG(IN,I) 1033 FORMAT(1X,13,3D26.14,3D 15.5) C C SAVE PLOT SCALING MAX AND MIN C IF(KLOG) GO TO 131 BIMX=AMAX1 (BIMX,SNGL(WIM(IN,I)) ) EIMN=AMIN1(RIMN,SNGI(WIM(IN,1) ) ) EBMX=AMAX1 (RRMX ,SNGI (WB (IN, I) ) ) RBMN-=AMIN1 (RRMN,SNGL (WR (IN,I) ) ) GO TO 100 131 R=SNGL(DABS (WR(IN,I))) IF (B.EQ.O. ) GO TO 132 RRMX=AMAX1(RRMX,R) BRMN= AMIN1 (FRMN, R) 132 R= SNGL (DABS (HIM (IN ,1) ) ) IF (B.EQ.O. ) GO TO 100 RIMX=AMAX1(RIMX,R) RIMN= AMIN1 (BIMN, B) 100 CONTINUE SOLVEQ=. FALSE. IF(.NOT.PLREAL) GO TO 300 WBITE(6,1011) RRMX,EFMN,RIMX,RIMN 1011 FORMAT(' RBMX,RRMN,RIMX,RIMN',4E15.7) IF (. NOT. KLOG) GO TO 201 RIMX=ALOG10 (RIMX) RIMN=ALOG10 (RIMN) \u00C2\u00A3BMX=ALOG10 ( B R M X ) BRMN=ALOG10 (BBM 8) 201 IE (. NOT. FANCY) GO TO 260 RMM=AINT (ALOG10 (ABS (BRMN) ) ) IM=RRMN/10.**RMM HMIN=IM*10.**RMM BMM= AlNT(ALGG10 (BRMX-RBMN)) IM= (RRMX-RRMN) /10. **RMM WDX=IM*10.3?*RMM C C PLOTTING C DO SCALING C PLOT AXES C IAE EL C PLOT ROOT LINES C APEIY SPECIAL SYMBOLS IE C REAL(8)<0, OF IMAG(W)>0. C GO 10 261 260 WMIN= ERMN WDX=(RRMX-REMN)/10. 261 IF (PINC. EQ.O.) PINC= (KMAX-KMIN)/1 0. IF(PKMIN.EQ.O.) PKMIN=K MIN INC=(NK-1)/20 177 IF(PRMIN.NE.O.) WMIN=Ffi MIN IF(PBINC.NE. 0. ) WDX=PBINC CALL AXIS(0.,0.,'HAVE NUMBEE* ,-11,10.,0.,PKMIN,PINC $ ) CALL AXIS{0.,0.,\u00C2\u00ABBEAL ANG FBEQ\u00E2\u0080\u00A2,13,10.,90.,WMIN,WDX $ ) IF (.NOT. KLOG) GO TO 262 CALL SYMBOL (\u00E2\u0080\u00940.3,9.0,.14,* LOG-LOG * ,90.,7) 262 IF (.MOT. LABEL) GO TO 263 CALL SYMBOL(i.,9.75,.14,ALABEI,0. ,80) 263 CONTINUE S=10./(NK-1.) DO 206 1=1,NK AX(I) =(1-1.) *S AY (I)=0. 206 CONTINUE IF (KLOG) GO TO 250 DO 202 IN=1,3 DO 203 1=1,KK AY (I) = (SNGL (WB (IN,I) ) -WMI N) /WDX 203 CONTINUE CALL LINE(AX,AY,NK,1) 202 CONTINUE CALL PLOT(12.,0. ,-3) GO TO 300 250 DO 212 IN=1,3 NIN=0 DO 213 1=1,NK IF (SB (IN,I) . EQ. 0. DO) GOTO 214 AY (I) =(SNGL(DLOG10 (DABS (WE (IN,I) ) ) ) -WMIN) /WDX IF(HB(IN,I).GT.0.DO) GO TO 213 IF (I. EQ. 1) GO TO 240 IF(WR (IN,1-1) .LT.O.DO) GO TG 241 NIN=NIN+1 KBEG(NIN)=1 KEND (NIN) = I GO TO 213 241 KEND (NIN) \u00E2\u0080\u0094 I GO TO 213 240 NIN=1 KBEG (NIN)=1 KEND (NIN) = I GO TO 213 214 AY(I)=0. 213 CONTINUE CAIL LINE (AX, AY, NK, 1) IF(NIN.EQ.O) GO TO 212 DG 248 K=1,NIN KB=KBEG(K) KE= KE ND (K) DO 249 I=KB,KE,INC CALL SYMBOL (AX (I) , AY (I) ,. 14,SYM (IN) ,0.0,-1) 249 CONTINUE 248 CONTINUE 212 CONTINUE CALL PLOT (12. ,0. ,-3) 178 300 IF(.NOT.PLIMAG) GO TC 399 DO 369 1=1,NK AY (I) =0. 369 CONTINUE IF{ .NOT. FANCY) GO TO 360 E M M= AINT(AL CG10 (ABS (BIMN) ) ) IM=EIMN/1Q.**8MM WMIN=IM*10. **.BJ1M BMM=AINT (ALOG10 (BIMX-EIMN) ) IM=(BIHX-BIMN)/10.**BMM WDX=IM*10.**BMM GO TO 361 360 WMIN= BIMN WDX=(BIMX-RIMN)/10. IF (1DX. EO. 0.) GO TO 399 361 CONTINUE IF (PIMIN. NE. 0. ) WMIN=PIMIN IF(PIINC.NE.O.) WDX=PIINC CALL AXIS(0.,0.,\u00C2\u00BBWAVE NUMEES\u00E2\u0080\u00A2,-11,10.,0.,PKMIN,PINC $ ) CALL AXIS(0.,0.,\u00E2\u0080\u00A2IMAG ANG FEEQ\u00C2\u00BB,13,10.,90.,WMIN,SDX $ ) IF(KICG) GO TO 350 DO 252 IN=1,3 DO 251 1=1,NK AY (I)= (SNGL "Thesis/Dissertation"@en . "10.14288/1.0085756"@en . "eng"@en . "Astronomy"@en . "Vancouver : University of British Columbia Library"@en . "University of British Columbia"@en . "For non-commercial purposes only, such as research, private study and education. Additional conditions apply, see Terms of Use https://open.library.ubc.ca/terms_of_use."@en . "Graduate"@en . "Radiation driven instabilities in stellar winds"@en . "Text"@en . "http://hdl.handle.net/2429/21496"@en .