Документ взят из кэша поисковой машины. Адрес оригинального документа : http://www.sao.ru/hq/lizm/conferences/pdf/1991/1992_1_p48-p58.pdf Дата изменения: Wed Feb 24 18:49:35 2010 Дата индексирования: Sat Sep 11 19:42:31 2010 Кодировка: Поисковые слова: m 17

"Stellar magnetism"

International

Conference Proceedings,

1992 ,

92-102 .

SYNT H - A CODE FO R FAST SPECTRA L SYNTHESI S N.E. Piskunov
Helsink i University , Finlan d

ABSTRACT We intend to describe here in full details ou r procedur e for calculating the synthetical stellar spectrum . It features high accuracy an d high computin g speed. In combinatio n wit h an interactive progra m for comparin g synthetical spectru m wit h the observations this procedur e prove d to be a usefull an d reliable tool for abundanc e analysis, spectral lines identification an d other practical purposes. Th e compute r implementatio n o f spectral synthesis - progra m SYNT H - i s written i n F O R T R A N 7 7 an d ca n b e installed o n several types o f computer s fro m 38 6 base d PC/A T machin e t o mainframes . ROTAT E -the progra m for compariso n the results o f SYNT H wit h the observations, is mad e independen t fro m a specific graphics interface whic h make s it portable to a different graphical environments . INTRODUCTION Th e ultimate goal of every spectral synthesis progra m is to extract som e astrophysical information fro m compariso n of the synthesized spectra to the observation. Therefore the calculations mus t carefully approximat e the process of radiative transfer throug h the atmospher e of the star an d also account for transformation of the spectru m inside the telescope-spectrograph-detector syste m use d in observations. An d so we are going to follow the pat h of the light, propagating fro m the internals of the star "toward s a computer" , whic h produce s a plot of residual intensities as a function of wavelength . At this point the compariso n of the observed an d compute d spectra becom e meaningfu l an d ca n privide astrophysical information. In our synthesis procedur e we consider only the absroption lines. Furthermor e we treat line formation in the local thermodynamica l equilibrium (LTE ) approximation , whic h basically mean s that w e assum e the temperatur e o f the radiation t o b e equal to the temperatur e of stellar matte r in every point in the atmosphere . Th e whol e procedur e ca n be extende d for the stars where LT E approximatio n is not valid, but it falls outside of the scope of this pape r (for non-LT E treatment see, for example , Aue r 1976, an d the bibliography therein). We will also restrict ourselves wit h the case of plane-parallel an d static atmosphere . Thi s mean s that the calculations will be less accurate for giants an d super giants an d they are probabl y not applicable to the fast expandin g envelopes. Anothe r restriction arised fro m ou r disk integration
92


procedure whic h assume s no rotational broadenin g of spectral lines an d the spherical shape of the star. Th e rotational broadenin g is importan t thoug h for compariso n with observations, so we included it in the cod e that does the compariso n using the approximation, suggested b y Gra y (1976) 1. RADIATIVE TRANSFER We are not going to repeat here the textbook s on the radiative transfer (see, for example , Sobole v 1960, Mihala s 1978, an d the references therein). Just a fe w definitions: · specific intensity I stands for the amoun t of energy, emmite d throug h a unit area aroun d a given point on the surface of the star at a given wavelengt h in a given direction. Thu s I is a function of wavelengt h , an d the direction of light propagation. I n the case o f plane-parallel atmospher e w e hav e a symmetr y relative to the radius pointing to the surface elemen t unde r consideration that is wh y only on e angle is sufficient to describe the direction. We expect the atmospher e o f the star t o b e homogeneous , an d the sam e mode l atmospher e (see below ) an d chemica l composition are valid for every surface element . So finally we have : , where : is the wavelengt h an d is a cosine of the angle betwee n the radius to a given point an d the direction of propagation of I. As we will see soon, is a very convenient paramete r for characterizing the direction. · monochromati c flux is the total amoun t of energy radiated by the star toward s the observer at a given wavelength . So ca n be calculated by integrating specific intensity over the visible hemisphere : (1) Ove r visible hemispher e is the projected elementary area on the stellar surface, which is getting smaller close to the limb. is a portion of the energy whic h is absorbe d · Continuou s opacity coefficient fro m a bea m of radiation wit h a specific intensity I in the elementar y volum e at a given wavelengt h in a bound-fre e an d free-free transitions processes. is only a function of the physical conditions inside In LT E approximatio n the elementar y volum e (chemical composition, gas an d electron pressures, an d temperature ) an d the wavelengt h . So if the radiation wit h specific intensity / will cross an elementar y volum e wit h opacity , it will lose · I of its energy.
93


Lin e opacity coefficient is th e sam e as continuous opacity coefficient except, that th e energ y is absorbe d in a bound-boun d transition processes. Optical dept h is a convenient substitute for geometrical distance scale in. stellar atmosphere . Optical dept h betwee n tw o points A an d is defined as:

(2)

wher e d x i s a n elementar y geometrica l distanc e alon g th e straight line betwee n A an d B . Th e su m represent s th e total absorptio n coefficient. Fo r furthe r convenienc e th e direction of th e positive optical dept h is set so , tha t if w e mov e fro m th e surfac e toward s th e cente r o f th e star, will increase . O n th e surfac e o f th e star = 0 .

Mode l atmospher e is a table of basic physical parameters , usually temperatur e an d electron an d gas density, compute d along som e standar d dept h scale i n th e stellar atmosphere . Mode l atmosphere s are suppose d t o give a n adequat e description of physical conditions in the atmospher e (unfortunately, it is not alway s true). A s a dept h scale on e usually us e a monochromati c optical dept h at a fixed wavelengt h (for example , ) or a colum n density (the total mas s i n a radial colum n wit h 1 cm 2 cross-section measure d fro m th e surface t o a given point i n the atmosphere) . Th e state-of-the-art model s hav e bee n calculated b y Kuruc z (see Kuruc z 1971 , Kurucz , Avret t an d Peytrema n 1974 , Kuruc z 1991a , Kuruc z 1991b) .
N o w w e ca n writ e th e basi c equatio n o f th e radiativ e transfer for specific intensity I :

(3) wher e - is a source function whic h in LT E approximatio n is equal to th e Plan k function an d thus depend s only o n an d temperature . Th e solution o f (3) for an y point M on th e surface of th e star ( =0 ) ca n be expressed in th e integral form : (4) an d substituting (4) to integral (1) we obtain th e expression for monochromati c flux radiated fro m th e atmospher e of the star at wavelengt h (for derivation see Mihala s 1978) :
(5 )

94


where

- is a second exponential integral.

W e solve equation 5 numerically using a n 8-th order Gaussia n quadratur e formul a (for description of numerica l integration wit h Gaussia n quadratur e see for example : Press et al. 1986) : (6) where an d , show n in Tabl e 1, are node s an d weights calculated for corresponding orthogonal polynomials. Th e only complication here is that are given in the monochromati c optical dept h scale , bu t , as a function of temperature , can be obtained fro m the mode l atmospher e only in the dept h scale of the model . Table 1

We ca n solve this proble m an d obtain the values of equation that links the monochromati c optical dept h scale
:

using a differential to the "standard" scale

if the dept h scale of the mode l atmo sphere is the standard optical depth , if the dept h scale of the mode l atmo sphere is the colum n density, wher e all the functions in the right part (K , an d ), ca n be foun d fro m mode l atmospher e an d therefore are know n as a functions of the mode l atmospher e dept h scale. Th e initial condition for equation 7 is set on the surface where : = =0 . Finally we will list the steps of synthetical spectru m calculations for a single wavelength:
std

95


· Start fro m the surface, that is wher e the monochromati c optical dept h · Integrate equation 7 fro m current a corresponding value of ; · Fin d the value of the Plan k function · Ad d t o the sum m 6 ; to the next ;

= 0;

fro m Tabl e 1 an d obtain

· Incremen t i an d repeat the procedur e for the next nod e until the deepest fro m Tabl e 1 is reached. Th e continuous opacity an d the line opacity are not given in the mode l atmosphere . The y should be calculated using those physical parameter s whic h are available fro m the model . We are going to discuss the opacities calculations in the next section. 2. OPACITIES CALCULATION

Thre e processes contribute to the continuous opacity: free-free an d bound-fre e atomi c transitions an d scattering on mos t abundan t particles in the atmosphere . In our calculations we include the following agents of continuous opacity: neutral hydrogen , helium , silicon, aluminium , magnesium , carbon, nitrogen an d oxygen ; positive ions H , He , Si , Mg , Ca ; negative ions H an d He ; the scattering o n free electrons, neutral hydroge n an d heliu m atoms . W e adopte d formulas an d interpolation tables fro m Kuruc z (1971) with som e corrections an d extensions. Line opacity is du e to the bound-boun d transitions in a particular ion. Several ions of different atom s ca n contribute to the line opacity at a given wavelength , so the total line opacity is a summ :
+ + + + + 2

(8) Lines wher e is the central opacity of line l an d H(a v) is the Voigt function whic h describes the wavelengt h dependenc e of the bound-boun d absorption:
y

(9) wher e the parameter s an d are defined as:

96


where : is the line dampin g paramete r an d is the Dopple r width . We use the approximatio n of Voigt function suggested by Dobriche v (1984), whic h combine s high accuracy an d computin g speed:

,(10) wher e Voigt widt h is defined as:
( )

includes therma l an d microturbulent broadening: (12) Th e value of the dampin g paramete r in mos t of the cases is determine d by three processes: radiative, quadratic Stark an d va n der Waal s damping : (13) Th e corresponding constants for a numbe r of lines are available fro m experiment s or quantu m mechanica l calculations. For other lines mor e simple approximation s of the impac t theory ca n be use d (see Gra y 1976) . Th e central opacity is proportional to the population of the lower energ y level by the atoms , responsible for the line formation, , to the oscillator strength f and to the stimulated emission factor: (14) wher e m an d e are the mas s an d the charge of electron. Th e population of the lower energy level is obtained fro m Sah a equation whic h gives the ratio of the numbe r of atom s in tw o consecutive ionization stages N lon an d Njon+1 as a function of the numbe r density of electrons N e an d the temperatur e T : (15) wher e Uion an d Uion+1 are the partition functions for tw o stages of ionization an d is the energy difference betwee n those tw o stages.

9 7


3. NUMERICAL CALCULATION OF THE MONOCHROMATIC FLUX
Th e equatio n 7 i s solve d startin g fro m th e surfac e wher e boundar y conditio n i s set . Th e integratio n i s carrie d usin g a 6-t h orde r Rung e - Kutt a procedur e described , fo r example , b y Pres s e t al . (1986 ) althoug h w e woul d no t recommen d t o us e th e cod e fro m thi s book . M u c h bette r F O R T R A N implementatio n ca n b e foun d i n N A G an d I M S L scientific softwar e libraries . W e ar e usin g a slightl y modifie d versio n o f a 6-t h orde r Rung e - Kutt a routin e fro m I M S L library . B y usin g 6-t h orde r algorith m w e achiev e tw o goals : th e accurac y ca n b e checke d withou t extr a integratio n wit h a smalle r stepsize , an d i f th e accurac y criteriu m i s no t satisfied onl y thre e mor e point s shoul d b e adde d t o obtai n th e integra l fo r hal f stepsiz e interval . I n eac h poin t o f th e integratio n w e hav e t o comput e th e righ t par t o f th e equa tio n 7 . Th e continuou s opacit y usuall y behave s quit e smoothl y through,th e atmosphere , s o w e calculat e i t i n advanc e fo r eac h poin t o f th e mode l an d us e th e paraboli c interpolatio n durin g th e integratio n o f "equatio n 7 . I n ou r cod e w e ar e us in g th e continuou s opacitie s subroutine s fro m Kurucz' s A T L A S - 5 progra m (Kuruc z 1971) . Lin e absorptio n m a y hav e ver y stee p changes . T o sav e som e computin g tim e w e tabulat e th e numbe r densit y o f all neutra l atom s an d ions , whic h produc e lin e absorptio n i n a give n spectra l interval , plu s th e numbe r densit y fo r hydroge n an d helium , neede d i n va n de r Waal s broadenin g calculations . Th e stepsiz e o f mode l atmospher e i s usuall y too bi g fo r thi s sor t o f tables , s o w e ad d extr a point s (sometime s u p t o 10 0 point s fo r on e ste p o f th e mode l dept h scale ) i n suc h a wa y tha t linea r interpolatio n provide s th e requeste d accuracy . Th e electro n numbe r densit y and , th e temperatur e ar e interpolate d fro m th e origina l model . Th e partitio n functio n hav e bee n pretabulate d b y Kuruc z an d w e ar e jus t usin g hi s subroutin e P F S A H A (Kuruc z 1971 ) fo r solvin g th e Sah a equation . Th e dampin g constant s ar e calculate d usin g Lindhol m approximatio n o f th e impac t theor y (Mihala s 1978 ) i f mor e accurat e value s ar e no t supplie d wit h th e othe r lin e parameters . Durin g th e integratio n o f equatio n 7 w e comput e th e centra l opacitie s (equatio n 14 ) an d th e Voig t functio n (equatio n 10 ) value s fo r eac h lin e an d evaluat e th e tota l lin e absorptio n coefficient usin g equatio n 8 . Afte r th e last nod e i n Tabl e 1 ha s bee n reached , w e ca n fin d th e valu e o f monochromati c flu x usin g quadratur e formul a 6 .

4.

WAVELENGTH GRID SELECTION

Th e avoi d redundanc y o f th e calculation s w e us e a procedur e fo r dynami c con structio n o f th e wavelengt h gri d insid e th e spectra l interval . Thi s procedur e mini mize s th e numbe r o f mes h point s b y choosin g th e variabl e stepsiz e whic h provides ] a fixe d accurac y o f linea r interpolatio n o f th e synthetica l spectru m fo r an y inside


our interval. Th e procedur e starts wit h computin g monochromati c fluxe s i n bot h ends of spectral interval. Nex t it add s a poin t in th e cente r of ever y spectral line. If tw o lines are closer tha n certain limit, onl y on e central wavelengt h (the shorte r one ) i s added . The n extr a point s ar e inserted i n th e middl e o f eac h pair o f neighbourin g lines. A t th e last stag e th e procedur e check s ever y poin t for th e accurac y o f linear interpolation an d adjust s stepsize, i f needed , th e kee p th e interpolation error unde r a requeste d limit. Let's see ho w t o estimat e th e stepsize whic h provide s a certain accuracy o f linear interpolation. I n a give n wavelengt h , whic h falls betwee n an d , th e monochromati c flu x i s approximate d as: (16) Tw o errors contribut e t o th e valu e o f Error : th e error o f linear interpolation itself wit h a finite differences an d th e error of approximatio n of th e first derivative of formula. Bot h o f the m ar e proportiona l t o . No w i f w e ad d a third an d let h b e th e stepsize betwee n an d , th e point i n th e middl e o f total error ca n b e estimate d as: (17) an d for a give n valu e o f Erro r th e stepsize ca n b e foun d fro m (17) . Not e tha t th e procedur e basically add s mor e point s i n place s wit h highe r curva ture (for example , line cores) . 5. CONTINUUM CALCULATION In order to be able to compar e the synthetical spectru m to the observations we need to normalaiz e it to the continuu m level. Th e integration of the continuous flux is not very different fro m the monochromati c flux integration. We take advantag e of a small an d smoot h variation of the continuu m level inside the synthesized interval (typically less tha n 10 0 ) an d calculate continuous flux only in the end s of ou r spectral interval. Fo r all other points continuu m level is obtained by linear interpolation. Anothe r difference fro m monochromati c flux calculation is the equation 7 which is transforme d to: if the dept h scale of the mode l atmo sphere is the standard optical depth , (18) if the dept h scale of the mode l atmo sphere is the column,density, because no line opacity should be included.
99


Th e procedur e of wavelengt h grid construction, described in the previous section is applied after the normalization to the continuu m level is done .
6 . A C C U R A C Y

Three parts of the synthetical spectrum calculations are crucial for the accuracy: construction of the line optical dept h scale (the integration of equations 7 an d 18), calculation of the monochromati c flux (equation 6) an d selection of the wavelength grid, described in section 4. Rung e - Kutt a algorithm, used for integration of equations 7 an d 18, ca n be tune d to an y required accuracy. We usually select 3 10-5, whic h is consistent with the accuracies achieved in the other parts. Flu x integration accuracy is mor e difficult to check. To do it we compare d our results wit h a direct disk integration of equation 1 an d we foun d that the 8-th order Gaus s quadratur e formul a provides the accuracy better the n 10 - 4 for a wid e range o f mode l atmosphere s wit h effective temperature s fro m 350 0 u p t o 2000 0 K . For the wavelengt h selection procedur e we usually pu t the accuracy criterium t o 10 -4 . Th e resulting synthetical spectru m ca n b e interpolated t o an y wavelength within a given spectral interval wit h the accuracy better the n 5 10 - 4 whic h corre spond s to signal-to-noise ratio of 2 000.
7 . L I M B D A R K E N I N G C A L C U L A T I O N

W e also comput e a continuu m lim b darkening coefficient, whic h ca n b e used later for convolution of the synthetical spectru m wit h the rotational profile. To do this we integrate equation 4 the sam e wa y as we did wit h equation 6 an d find the values of specific intensity for continuu m at 3 points on the disk at 0, 0.6 an d 0.95 stellar radius fro m the center. Th e lim b darkening coefficient is the n foun d from least square approximatio n of the conventional formula:
8 . C O M P A R I S O N W I T H T H E O B S E R V A T I O N S

Th e compariso n wit h the observations is don e wit h a separate progra m whic h includes an interactive graphics interface. It show s the observed an d the calculated spectra o n the screen an d wit h the help o f men u on e ca n convolute the synthetical spectru m wit h the rotational profile, the macroturbulenc e profile (using the approx imatio n formul a fro m Gra y 1976 ) an d the instrumental profile (Gaussian) . Fo r the rotational profile w e followed the approximatio n suggested b y Gra y (1976): (19)

10 0


wher e th e kerne l o f th e convolutio n i s a functio n o f lim b darkening :

(20 )

Th e synthetica l spectru m ca n als o b e shifte d i n wavelengt h an d zoome d fo r fin e tunin g o f th e shifts an d broadenings . Th e equivalen t widt h calcutatio n simplifie s th e estimat e o f chemica l abundances .

9. COMPUTER IMPLEMENTATION
Bot h synthetica l spectru m calculatio n cod e ( S Y N T H ) an d th e progra m fo r com parin g th e synthetica l spectru m t o th e observation s ( R O T A T E ) ar e writte n i n con ventiona l F O R T R A N 7 7 an d w e hav e bee n usin g the m o n AT/38 6 typ e persona l computer . Compile d wit h N D P F O R T R A N compiler , S Y N T H take s approximatel y 2 second s pe r on e wavelengt h point , whil e computin g 2 5 spectra l line s o n a 2 0 M H z machine . Th e onl y critical requiremen t wa s th e presenc e o f expende d memory , bu t wit h 4 Mbyte s o f tota l m e m o r y siz e w e wer e abl e t o calculat e 3 0 A interva l wit h 8 0 spectra l lines . Bot h program s ca n produc e colou r graphic s usin g N D P L O T library . T o mak e th e cod e less dependen t o n graphic s packag e w e m a d e a simpl e inter mediat e set o f interfac e subroutines , whic h includ e all scaling s an d calls t o graphic s library functions . S o wit h simpl e substitutin g o f thes e calls w e ha d n o problem s runnin g S Y N T H a s batc h jo b o n th e V A X (wit h n o graphics ) an d R O T A T E a s a n XWindow s client o n Su n workstatio n wit h all bell s an d whistle s availabl e there . S Y N T H require s tw o sort s o f inpu t data : line s list an d mode l atmosphere . Line s list include s th e spectra l interva l t o b e synthesize d an d th e numbe r o f line s followe d b y th e parameter s fo r eac h line . Line' s parameter s includ e th e n a m e an d ionizatio n stag e o f chemica l element , centra l wavelength , excitatio n energ y (e V o r c m ) , mi croturbulen t velocit y (km/s) , lo g g f an d thre e dampin g constant s (radiative , Star k an d va n de r Waals ) i n exponentia l o r logarithmi c for m (se e exampl e below) . I f som e o f th e dampin g constant s ar e no t known , the y shoul d b e se t t o 0 an d S Y N T H wil l us e th e approximatio n formul a fo r them . Her e i s a n exampl e o f on e spectra l lin e description :
-1

'F E 1' , 6165.36 , 33412.71 , 0.5 , -1.670 , 7.943 , -6.156 , -7.83 3 Th e filenam e o f th e mode l atmospher e dat a i s include d wit h th e line s list. Usuall y mode l atmospher e include s th e chemica l compositio n data , use d durin g th e mode l calculation . W e rel y o n thi s abundance s value s a s a defaul t fo r spectra l synthesis , bu t the y ca n b e change d fro m th e line s list. Fo r example , t o chang e th e defaul t abundanc e o f iro n an d nike l on e shoul d ad d th e followin g lin e t o th e line s list: 'FE:-4.4','NI:-5.8','END ' Mode l atmospher e dat a beside s th e defaul t abundance s shoul d includ e th e optica l

101


dept h scale, temperature , gas an d electron numbe r densities an d the density of stellar matter . A s the result SYNT H produce s a file , whic h contains the wavelengt h grid, the corresponding residual intensities an d the lim b darkening coefficient in a forma t suitable for ROTATE . ROTAT E accepts the dat a produce d b y SYNT H an d the fil e wit h observations (if any ) in the sam e forma t as the synthetical spectrum . Th e results of interactive fitting (rotational, instrumental an d macroturbulen t broadenin g an d wavelengt h shift) wit h ROTAT E are stored i n the outpu t file , whic h ca n by use d for plotting or further processing. REFERENCES Auer , L.: 1976, J.Q.S.R.T, 16 , 931 Dobrichev , V.: 1984, Comptes rendus de L'Academic bulgare des Sciences, 37(8) , 991 ) Gray , D.F.: 1976, The Observation and Analysis of Stellar Photospheres. Ne w York , Joh n Wile y an d Son s Kurucz , R.L.: 1971, S.A.0. Special Report No. 309. Cambridge , Mass . Smithsonia n Astrophysical Observator y Kurucz , R.L.: 1991a, Harvard-Smithsonian Center for Astrophysics Preprint Series No. 3181. Kurucz , R.L.: 1991b , in Stellar Atmospheres: Beyond Classical Models, eds. L. Crivellari, I . Hubeny , D.G . Hummer , N A T O AS I Series, Kluwer , Dordrecht , i n press Kurucz , R.L., Peytremann , E . an d Avrett, E.: 1974 , Blanketed Model Atmospheres for Early-Type Stars. Washington , D.C. , Smithsonia n Institution Mihalas , D. : 1978 , Stellar atmospheres. Sa n Francisco, Freema n an d Co . Press, W.H. , Flannery, B.P., Teukolsky, S.A. an d Vetterling, W.T. : 1986, Numerical Recipes. Cambridge , Cambridg e University Press Sobolev, V.: 1960, Moving Envelopes of Stars. Cambridge , Mass. : Harvar d Univer sity Press.

10 2