Determination of the mass loss rate and the terminal velocity of stellar winds. I Genetic algorithm for automatic line profile fitting

0.35 MB·English

by  L. Georgiev

To appear in “Conference Title (2001)” RevMexAA(SC)
DETERMINATION OF THE MASS LOSS RATE AND THE TERMINAL
VELOCITY OF STELLAR WINDS. I GENETIC ALGORITHM FOR
AUTOMATIC LINE PROFILE FITTING
Leonid Georgiev,1 Xavier Hernandez,1
RESUMEN
Se presenta un nuevo m´etodo de ajuste automa´tico de perfiles P Cyg observados en espectros UV de vientos
estelares. La funcio´n fuente de las l´ıneas se calcula usando la aproximaci´on de Sobolev y el flujo emitido se
5
0 obtiene atraves de la soluci´on exacta de la equaci´on de transporte (similar al m´etodo SEI descrito por Lamer
0 et al. 1987). La calidaddel ajuste se evalua usando el estimador Likelyhood. La maximisaci´ondel Likelyhood
2 se hace atraves de un algoritmo gen´etico. Las ventajas de nuestro m´etodo comparado con otros m´etodos

n similares son su estabilidad y su poca sencibilidad de las condiciones iniciales. Ademas, el m´etodo garantiza
a la localizaci´on del maximo global de la superficie del Likelyhood, cual no se puede hacer con los algoritmos
J
cl´asicos. Presentamos la implementaci´on del metodo, las pruebas de su funcionalidad usando datos sint´eticos

8 y reales y la estimacio´n de los rangos de confianza de los resultados.
2

ABSTRACT
1
v A new method for automatic fitting of P Cyg line profiles in UV spectra of stellar winds is presented. The
9 line source function is calculated using Sobolev approximation and the emergent flux is obtained by exact
3
integration of the equation of the radiation transport (similar to the SEI method described by Lamers et al.
6
(1987)). ThequalityofthefitisevaluatedusingtheLikelyhoodestimator. ThemaximizationoftheLikelyhood
1
0 is done by a genetic algorithm. The advantagesof our method with respect to other similar approachesare its
5 robustnessanditsinsensibilitytotheinitialguess. Inaddition,thealgorithmguaranteesthelocalizationofthe
0 global maximum of the Likelyhood hypersurface, which is not the case for classical minimization algorithms.
/
h Here we present an implementation of the genetic algorithm for line profile fitting, its tests on both synthetic
p and real data and and estimation of the confidence limits of the results.
-
o Key Words: — — —
r
t
s 1. INTRODUCTION the fitting of the spectrum requires high quality
a data which covers a wide range of wavelengths and
: Many hot stars lose mass through a supersonic
v with high dispersion, this again limits the method
wind. The terminal velocity of these winds, and the
i to bright or nearby stars. An alternative has been
X details of the mass loss rate play an important role
proposed by Lamers et al. (1987) and Groenewe-
r instellarevolution. Stellarwindscruciallydetermine
a gen & Lamers (1999). This method is based on a
theinteractionofthestarwithitssurroundinginter-
code which approximates the opacity of the wind
stellarmedium,thisbeingoneofthecentralfeedback
as a simple function of the distance from the stel-
mechanisms, which in turn modulate star formation
lar nucleus. The code calculates the source func-
at galactic levels. Estimating stellar mass loss rates
tion in a given line within the Sobolev approxima-
is therefore a highly relevant problem.
tion (Sobolev, 1957) and then calculates accurately
The best way of determining the mass loss rate
the emitted profile. The fit of the observed profile
of a star is by studying the radio emission from
givestheterminalvelocityofthewindtogetherwith
the wind itself (Lamers et al. 1999) however, with
the total optical depth along the line of sight. The
current observational techniques this is limited to
mass loss rate is estimated from the later (see next
nearbystarswithhighmasslossrates. Theotherre-
section). Eventhoughtheobtainedvalueissensitive
liable method is basedona detailed modeling ofthe
tocertainassumptionsregardingtheionizationfrac-
stellar spectrum. The recent development of realis-
tion of the element whose line is fitted, the method
ticcodeswhichincludemanyatomiclinesmakesthe
has the advantage of being applicable to faint stars
latter approach the preferred one. Unfortunately,
forwhichonly lowresolutionandlowsignalto noise
spectra are available.
1InstitutodeAstronom´ıa,UNAM,M´exico.
1
2 GEORGIEV & HERNANDEZ
The main disadvantage of this approach is the opacity of the line, ν is the laboratory frequency of
0
large number of free parameters which have to be the transition and c is the speed of light. T, α and
1
adjusted. Groenewegen & Lamers (1989) proposed α are free parameters of the model, which describe
2
a method for determination of the parameters by a the radial dependence of the optical depth. The ve-
”classical” minimization of the χ2 = (1 F /F ) locity field of the wind is adopted to be a standard
o c
−
where F is the observed flux and F is the calcu- β-law
o c
P
lated one. But the multi-dimensional fitting con-
β
vergestothecorrectsolutiononlyifthe initialguess R0
is close to the χ2 minimum. The complex topology V(r)=V∞ 1− r (2)
(cid:18) (cid:19)
associated with the high dimensional hypersurface
where V∞ is the terminal velocity. The winds
being treated leads to the appearance of numerous
of the hot stars are not smooth but have random
local minima. The “classical” minimization meth-
motions. The shape of the blue wing of the P Cyg
ods usually get trapped in the local minimum close
absorptioncomponentsuggestthatthe intrinsic line
totheinitialguessmissingtheglobalone. Thesome-
profile has a width on the order of hundreds of kilo-
what lengthy numerical procedure involved, again
meters per second. The physical nature of these
together with the large number of dimension of the
random motions is not clear but they act as an ad-
problem, makes a dense sampling of the parame-
ditional line broadening mechanism similar to tur-
ter space completely impractical. Recent years have
bulence and are frequently called ”turbulence” al-
seenthe developmentofnon-standardmaximization
thoughtheymayhavelittletodowithit. Wemodel
techniques such as simulated annealing and genetic
this chaotic motions assuming the intrinsic line pro-
algorithms, particularly in connection with astro-
file to be Gaussian defined as
physical applications (e.g. Teriaca et. al 1999 in
fitting radiativetransfermodels to solaratmosphere 1 − v 2
data; Sevenster et al. 1999 in fitting parameters of φ(v)= √piV exp Vturb (3)
turb
Galactic structure models to observations or Met- (cid:0) (cid:1)
The turbulent velocity, V , can be variable
turb
calfe 2003; for the fitting of White dwarf astroseis-
throughout the wind, but because its physical na-
mologyparmeters). Thesemethodsprovideaframe-
ture is not clear there is no model of its possible
workwhichina reliablemanneridentifies the global
variability. To keep the model as simple as possible,
maxima (or minima, as the case might be) of com-
we assume V to be constant and a parameter of
turb
plex multi-dimensional surfaces, with a close to op-
the model.
timal use of computing resources.
Finally, we set the innermost point of the grid,
In this paper we present an application of a ge-
R to the location where V(R ) is equal to the
min min
netic algorithm used to obtain the optimal param-
sound speed V . This is done by setting the pa-
sound
eters of a stellar wind, within the Sobolev approxi-
rameter R in (2) as
0
mation, for the fitting of P Cyg line profiles of UV
resonance lines. V β1
sound
Section (2) describes the radiative transfer R0 =Rmin 1.0 (4)
model. In section (3) a description of the genetic " −(cid:18) V∞ (cid:19) #
algorithmisgiven,withuseofitinanumberofsyn- Equations (1), (2) and (4) determine the distri-
thetic test cases and real data examples appearing butionofthecentralopacityofthelineχ asafunc-
0
in section (4). Section (5) presents our conclusions. tionofthe radiusr. Within the Sobolevapproxima-
2. THE ALGORITHM tion, one can now calculate the line source function
2.1. Radiative transfer code S and then can apply a formal solution of the equa-
Thecoreofthemethodisacodewhichsolvesthe tion of radiative transportthroughout the wind and
radiativetransportproblemundertheassumptionof calculatetheemergentflux(seeGeorgiev&Koenigs-
a spherically symmetric stellar wind. Following the berger,2004for more details ofthe code). The code
SEI code (Lamers et al. 1987) we approximate the is designed to workin 3D geometry,but for the cur-
radial Sobolev optical depth τ(v) as: rent test purposes, the solution is restricted to the
caseofsphericalsymmetry. Finally,wecalculatethe
χ c
τ(v)= 0 =T vα1(1 v)α2 (1) mass loss rate as follows (Groenewegen & Lamers,
ν0ddvr − 1989). The line opacity χ0 is
where v = VV(∞r) is the velocity of the wind, nor- πe2 ν0dV(r)
malizedtotheterminalvelocityV∞,χ0isthecentral χ0 = mcflunl =τ(v)∗ c dr , (5)
GENETIC ALGORITHM ... 3
where n is the population of the lower level, f point is taken from the signal-to-noise ratio of the
l lu
is the oscillator strength and we do not take into observations. Setting a large σ to certain points
i
account the correction for stimulated emission. The forces them to participate less in the final fit, which
level population n can be written as allowsustoexcludeeffectivelyspectralregionswhere
l
overlap with other lines occur. Hence, the final re-
nl =Nionqex =Natomqionqex =NHAEqionqex, (6) sult depends on the assigned σi. In practice, the
bandsofthe observedspectrumaffectedbyspurious
where q and q are the ionization and excitation
ion ex details are easily detectable and the assignment of
fractions and A is the chemical composition of the
E the weights is straightforward.
element relative to Hydrogen by number. Then the
Thesyntheticprofileiscalculatedoveragridad-
level population n is related to the mass loss rate
l justed for better sampling the velocity field in the
by the equation of continuity
wind. The calculated profile is then interpolated on
A M˙ q q the wavelengths of the observed profile using mono-
E ion ex
n = , (7) tonic cubic polynomials (Steffen, 1990). Finally,
l µ 4πr2V(r)
equation (9) is calculated and used as a measure of
where µ is the average molecular weight. Substitut- the quality of the fit.
ing(7)into(5),onecanobtaintheproductM˙ q q
ion ex
at each point r as 2.2. The genetic algorithm
Once we have a method for calculating the line
dV(r) µ 4mc ν
M˙ q q =τ(v) r2V(r) 0. (8) profile which corresponds to a particular set of six
ion ex dr A e2f c
E lu parameters,andan estimator ofthe goodness of the
Then, if the parameters V∞,Vturb,β,T,α1 and fit,wemustnowdeterminewhattheoptimalparam-
α are given,equation (1) and (2) together with the eter vector is; i.e. we must find the combination of
2
atomic data andthe chemicalcompositionallowthe sixparameterswhichmaximizesourLikelihoodfunc-
estimation of the product M˙ q q . The method tion. This will yield the physical parameters of the
ion ex
described in this paper suggests a better way to ob- stellar wind associated to the observed line profile,
tainthesixparametersthatdescribethevelocitylaw from which the mass loss rates and wind terminal
andtheopticaldepth. Itmaynotimprovethepreci- velocities can be estimated.
sion of the calculatedmass loss rate if the dominant A dense exploration of our six-dimensional pa-
uncertainty is the ionization fraction. rameter space is unfeasible. Classical maximization
Applying the above procedure, we are able to techniques arealsounreliablein connectionto prob-
calculate any line profile, once the six parameters lems of such high dimensionality, and liable to get
of the method V∞,Vturb,β,T,α1 and α2 have been trapped by local maxima (Press et al . 1989). In
specified. The next stepis the selectionofastatisti- fact wehaveencounterednumeroustest caseswhere
cal estimator of the goodness of fit to be associated local maxima exist in the Likelyhood surface of our
to a particular set of parameters and a given ob- problem. Webypasstheseobstaclesusinggenetical-
servedprofile. From a Bayesianperspective, the op- gorithms in searching for the maximum of the Like-
timal suchindicatoris the Likelyhood, calculatedas lyhood surface.
the probability that a certain data set be observed, Theideaistosimulateapopulationoforganisms
given a model. We hence compare a modeled pro- (wecallthem”bugs”)whichbreedandevolvefollow-
fileagainstagivenobservationthroughaLikelihood ingprescriptionsbasedonthatbiologicalsystemsare
estimator defined as: thought to follow, the result being a progressive in-
creaseinthefitnessofthepopulation,withthefittest
individualineachgenerationeventuallyreachingthe
logL(V∞,Vturb,β,T,α1,α2)= (9) absolute maximum of the fitness surface. For the
2 case at hand, a ”bug” is a set of wind parameters
1 F F
i,o i,c
log exp − described above.
i "(2π)1/2σi −(cid:18) 2σi (cid:19) # Once an observed line profile has been picked,
X
. the firststepistoselectN randompointsinourpa-
Where F and F are the observed and calcu- rameter space, each a parameter vector for the stel-
i,o i,c
lated flux respectively at a certain wavelength, λi lar wind, (V∞,Vturb,β,T,α1,α2). The ranges over
andthefactorsσ aretheobservationalerrorsatthe which these parameters are to be chosen have to be
i
same λ . The weight σ assigned to each observed determinedbyinspectionoftheobservedlineprofile.
i
4 GEORGIEV & HERNANDEZ
Onlyamplemarginsenclosingtheexpectedvalueare “inferiors”mate, only one child ensues. The process
needed. We turn to this issue in the examples sec- isrepeateduntilNoffspringhavebeenobtained,the
tion. new generation has been constructed, and it com-
Thesixcoordinatesofeachpointarethenturned pletely replaces the former.
into binary numbers, with a pre-determined resolu-
The parameter S hence represents a selection
tion and then combined in a string. We have used 8 P
pressure. If chosen close to 1, most individuals will
bitsperparameterforourimplementation,resulting
be deemed “superior”, and hence the second gener-
in a string of N =6 parameters 8 bits =48 bits.
G × ation will be distributed in parameter space much
This selection gives to our search algorithm a res-
like the first one was. However, as S takes larger
olution of Range /28 in each of the six dimensions P
i values,the secondgenerationwillbe made upofthe
where Range1−6 are the ranges for each parameter sons of (mostly) the “superior” individuals of the
chosen above.
former, yielding a gradual climb up the Likelihood
Eachofthe N randomlyselectedpointsbecomes
surface of the problem. Next the genotypes of the
a “bug” of the first generation. The string of N 1s
G new generation are turned back into coordinates in
and 0s whichcorrespondsto its 6 coordinatesin pa-
our parameter space, and then into phenotypes -the
rameterspacebecomesits“genotype”,andtheLike-
associated Likelihood. The ranking is repeated, and
lyhoodassociatedto that pointbecomes its “pheno-
the cycle iterated.
type”. The N members of the first generation are
thenrankedaccordingtothe valueoftheLikelihood In this way, we have introduced mating, selec-
associated to each, with those ranking above N/S tion, and chance mutations, the three main ingredi-
P
being deemed “superior”, and those ranking below ents of biological evolution, albeit in a highly sim-
this threshold “inferior”. SP is a parameter of the plified manner. Selection forces subsequent gener-
simulation which has to be a number greater than ations up the Likelihood surface, with the random
1.0, typically around 1.5. mutationsensuringthatafractionoftheindividuals
This completes the first generation, which are are constantly sampling new regions, which in turn
now “mated” to produce the second one. This pro- guaranteestheabsolutemaximumwilleventuallybe
ceeds bythe randomselectionofpairsofindividuals found.
from the pool of N members of the first generation.
In many classical maximization algorithms one
Once a pair is selected, “offspring” are calculated.
is called to evaluate the gradient of the surface at
The first C genes of one of the parents are taken
P
anygivenpoint. Whenthe surfaceisnotananalyti-
for first C genes of the offspring, and the remain-
P
cal one, but the result of a lengthy numericalproce-
der are takenfrom the other parent. C is a “cross-
P
dure, as in our case, the gradients could quite easily
ingpoint”,chosenatrandominthe interval(0-N ),
G
be dominated by roundoff errors. The method de-
each time an offspring is constructed.
scribed calls for no gradients of the surface, only its
In this way, the offspring will have “genes” in
valueatthepointsbeingtested. Withtheresolution
some aspects resembling those of one parent, and in
described here, a dense sampling of our Likelihood
some cases the other, the corresponding position in
surface would have required (28)6 =2.8 1014 eval-
parameter space will hence share some coordinates ×
uations of our wind model. We have used N=100
with one parent, and some with the other, with one
members per generation, and both in synthetic ex-
of them being a mixture. Also a fraction M of
F ampleswheretheanswerwasknowninadvance,and
randomly selected genes in the new generation are
in tests with real data, have found convergence in
“mutated” i.e. they are flipped from 1 to 0, or 0 to
200 generations, i.e. only 2 104 costly line profile
1, as the case might be. MF is a second and last calculations. ×
parameter of the method, which together with S
P
must be chosen and tuned through some testing of The method described here differs only slightly
the algorithm in the context of the particular prob- from the general purpose genetic maximization al-
lem being treated. We have used MF close to 0.01 gorithmPIKAIA,describedinCharbonneau(1995),
which means 1% of the genes are mutated andhavingbeenusedsuccessfullyinanumberofas-
The number of offspring to be produced by each trophysical applications to date (e.g. Teriaca et. al
pair depends on the fitness of the parents, if both 1999, Metcalfe 2003). The differences between the
are“superior”individuals,3offspringarecalculated. twoarelimitedtotheway“selection”istreated,the
Thepairingofa“superior”andan“inferior”individ- implementation described here giving results more
ualyieldsonly2offspring,andwhentwounfortunate suitable for the particular problem we are treating.
GENETIC ALGORITHM ... 5
-7.2 -7.1 -7 -6.9 -7.2 -7.1 -7 -6.9 -7.2 -7.1 -7 -6.9 -7.2 -7.1 -7 -6.9
2060 2060
2040 2040
2020 2020
2000 2000
1980 1980
1960 1960
1.1 1.1
1 1
0.9 0.9
0.8 0.8
0.7 0.7
-7.2 -7.1 -7 -6.9 -7.2 -7.1 -7 -6.9 -7.2 -7.1 -7 -6.9 -7.2 -7.1 -7 -6.9
Fig.2. Theellipsesshow1−σconfidenceintervalsonthe Fig.3. Theellipsesshow1−σconfidenceintervalsonthe
recoveredparameters(circles)forthefittothesynthetic recoveredparameters(circles)forthefittothesynthetic
profile with S/N=100. The input values are shown by profilewithS/N=50. Theinputvaluesareshownbythe
thesolid dots. solid dusts.
3. TEST CASES -7.2 -7.1 -7 -6.9 -7.2 -7.1 -7 -6.9
2060
The first test of the algorithm is made on syn-
2040
thetic data. We calculate the line profile of the
CIV1548/1550˚A doubletwithapredefinedvelocity 2020
law and wind density distribution and then fit this 2000
profile using the genetic algorithm. Fig. 1 showsthe
1980
original profile and the resulting fit for three differ-
ent signal-to-noise ratios. The parameters used to 1960
generate the profile and the result of the fitting are 1.1
summarized in Tab. 1.
The genetic algorithms do not include an easy 1
way of estimating the errors of the results. The
methodonlyservestolocatetheglobalmaximumin 0.9
an efficient manner. We estimate the confidence in-
tervals of the recovered parameters by Monte-Carlo 0.8
simulationsof synthetic line profiles. Inthis way,all
uncertainties inherent to the method and the mod- 0.7
-7.2 -7.1 -7 -6.9 -7.2 -7.1 -7 -6.9
elled data are taken into account. We select three
signal-to-noise ratios: 100, 50 and 10. For each of
themwe generate20differentline profiles,usingthe
Fig.4. Theellipsesshow1−σconfidenceintervalsonthe
same synthetic data but adding different and inde-
recoveredparameters(circles)forthefittothesynthetic
pendent noise. The obtainedline profileswerefitted
profilewithS/N=10. Theinputvaluesareshownbythe
by the genetic algorithmusing differentfirst genera- solid dusts.
tions of ”bugs”. Figs. 2,3 and 4 show the confidence
intervalsobtainedthroughthismethodfortherecov- variance analysis is performed, yielding the param-
ered parameters for the simulated profiles of Fig. 1 eter’s standard deviations, which together with the
Once the Monte-Carlo simulations are all done, relevant covariancescan be used to construct statis-
we find the means for all the recovered parameters, tically meaningful confidence ellipses. Fig. 2 shows
andofthe recoveredmasslossrates. Next, afull co- the 1 σ ellipses for the caseof S/N =100,with the
−
6 GEORGIEV & HERNANDEZ
Fig. 1. Test of thefittingalgorithm. A syntheticprofileofCIV1548/1550˚A (squares) isshown together with thebest
fit (solid line) for threeadopted S/N ratios.
soliddotsindicatingtheinputvalues,andthecircles
showing the ones reached by the fitting algorithm.
The mass loss rate is recoveredalmost exactly, with
a 1 σ interval of 30%. A negative correlation is
seen−between the inferred M˙ and the recovered ter-
minalwindvelocity,whilemasslossrateandβ show
aweakpositivecorrelation. Figs.3and4areequiva-
lenttoFig.2,butforthecaseswithS/Nvaluesof50
and 10, respectively. The same qualitative features
can be seen, with the ellipses getting only slightly
larger, to reach 40% for the S/N =10 case. It is re-
assuringthaterrorsontheinferredmasslossrateare
not highly sensitive to the value of the S/N of the
simulatedprofile,makingthemethodusefullevenin
cases where observationsare not of the best quality.
Also, we note that no systematics appear on the re- Fig. 5. Fit of the C IV 1548/1550˚A line in ζ Pup
coveredM˙ , we have constructed unbiased estimator (=HD66811). The thin line is the observed spectrum.
oftheimportantphysicalparameterwesetouttoin- The dotted line is thefit
fer. Wehavetostressthattheerroronthemassloss
rate are the formal errors of the fit. The errors re- the profile of Si IV 1398/1402˚A doublet in the IUE
lated to the determination of the ionization fraction spectrum of HD 30614. The comparison between
will be treated in the next paper of the series. the results of our fit and the parameters obtained
by Groenewegen & Lamers (1989) (Tab. 2) shows a
As can be seen from Tab. 1 in all cases the al- good agreement. Our solution shows a slower veloc-
gorithmfinds thecorrectmaximumofthelikelihood ity law (large β) but our final fit has higher Likely-
estimator, the mass loss rate is recovered success- hood than the fit obtained with the published pa-
fully, always well within 1 σ of the input value, rametes. For both stars we obtain a larger V∞ and
−
even at low S/N values. As an additional check we lowerV . AsshowninFigs.6,thetwoparameters
turb
triedtofitthelineprofilesofrealstars. Fig.5shows are correlated and one could expect this behavior.
the profile of C IV 1548/1550˚A doublet observed in Butinbothcasestheautomaticfitgivesvaluesvery
IUEspectrumofζ Pup(HD668110)andFig.7shows close to the more human controlled solution.
GENETIC ALGORITHM ... 7
TABLE 1
PARAMETERS OF THE SYNTHETIC PROFILE AND THE RESULTS OF THE FITS.
Parameter Exact value S/N=100 S/N=50 S/N=10
V∞ (km s−1) 2000 1997 19 2007 23 2014 30
± ± ±
V (km s−1) 100 114 12 107 13m 101 9
turb
± ± ±
β 1.0 0.99 0.06 0.99 0.05 0.87 0.08
± ± ±
T 2.0 2.01 0.06 1.96 0.08 2.015 0.11
± ± ±
logM˙ qionqex (M⊙/yr) -7.0 -7.00 0.13 -7.00 0.14 -7.02 0.15
± ± ±
TABLE 2
PARAMETERS OF THE WIND OF HD30814 AND HD66811 DETERMINED IN THIS WORK
TOGETHER WITH THE SAME DATA OBTAINED BY GROENEWEGEN & LAMERS (1989).
Parameter HD66811 C IV 1548/1550˚A HD30814 Si IV 1398/1402˚A
This work Published This work Published
V∞ (km s−1) 1550 16 1550 50 2137 21 2200 60
± ± ± ±
V (km s−1) 240 24 190 70 287 30 290 70
turb
± ± ± ±
β 1.13 0.1 0.7 0.1 1.18 0.1 0.7 0.1
± ± ± ±
logM˙ qionqex (M⊙/yr) -6.99 0.15 > -7.14 -7.27 0.15 > -8.4
± ±
120
100
80
Fig. 7. Fit of the Si IV 1392/1402˚A line in HD 30614.
The dotted line is thefit
1960 1980 2000 2020 2040 2060
robustness of the method even in case of low S/N
data. The method is especially useful for objects
Fig. 6. Correlation between V∞ and Vturb calculated for for which only low resolution and low S/N data are
syntheticprofile with S/N=50.
available. Their spectra show only few spectral fea-
tures and the application of full NLTE codes is not
4. CONCLUSIONS AND DISCUSSION
practical. Our method was successfully applied by
The determination of the basic stellar parame- Arrieta&Stangellini(2004)forcentralstarsofLMC
ters like wind velocity and mass loss rate is a very planetary nebulae. Also, the use of Monte-Carlo
time consuming process. In this paper we present a simulations on the recovered parameters allows for
method which allows an objective determination of the secure and objective determination of statisti-
these parameters based on a completely automatic cally meaningfull confidence intervals around recov-
procedure. The test cases shown above support the ered parameters, indeed, the full covariance matrix
8 GEORGIEV & HERNANDEZ
quentlyhaveoverlappedabsorptionlinesofotherel-
ements, interstellar lines or absorption components
ofthesamestarandthesameelementbutwhichare
not formed in the wind. The weights σ should be
i
chosenlowerintheaffectedbandssotheLikelyhood
is calculated using only the pixels which belongs to
the line. Fig 8 show the fit of the P V 1117/1128˚A
line in the star HD5980 as an example of the per-
formance of the code on a highly contaminated line.
TheFUSEspectraarecharacterizedbyahugenum-
ber of interstellar absoption line. In the above case,
a small weight (large errors) were assigned to the
bands withhighcontamination. As seenfromFig8,
Fig. 8. Fit of the P V 1117/1129 line in the LMC star the code finds a satisfactory fit which reproduce the
HD5980. The thin line is the observed spectrum. The
P Cyg line profile, avoiding the contamination.
dashed line is the fit. Small weight were applied to the
Another problem is that the parameters can be
interstellarabsorptionline,whichcontaminatethePCyg
obtainedonlyfromagrid. Asaresult,thecodenever
profile.
convergestotheexactsolutionbuttothenodeofthe
parameter’sgridwhichisclosest,butnotexactlyat,
isavailable. But,wehavetostress,thattheparame-
the maximum ofthe likelyhoodsurface. Thatmight
terobtainedbythismethodisnotthemasslossrate
M˙ itself, but rather the product M˙ q q . There cause the systematic error seen in Fig. 4. One pos-
ion ex
sible solution to this problem is a second run of the
is no easy way to separate these quantities on the
code with smaller ranges of the parameters around
base of low dispersion and low S/N data on which
the solution found. Another approach is to com-
only few lines are measurable. Our method helps to
bine a genetic algorithm with a classical minimiza-
obtainamorereliableestimatesoftheproductinan
tionmethod. Thegeneticalgorithmguarantees,that
easier way, but the final decision on how to disen-
tangleM˙ andtheionizationandexcitationfractions the solution is near the global maximum of the like-
lyhood surface and then the classical minimization
depends on the studied object and on the problem
algorithmstarts fromthat point andfinds the exact
which needs to be solved.
solution. Until now we have explored only the first
An important advantage of our method, com-
option. Thecombinationofthealgorithmsandtheir
pared to other fitting algorithms is a minimal effect
applicationtoawiderangeofobjects,togetherwith
of the initial guess on the final solution. This is be-
theproblemofseparatingthemasslossratefromthe
causethemethoddoesnotrelyonaseriesexpansion.
ionization fraction will be subject of the next paper
Thegeneticalgorithmrequiresonlyacrudemargins
of this series.
on the parameters to find the maximum of the like-
lyhood surface. In case the maximum is not within
We thank Gloria Koenigsberger for many use-
the proposedranges,the population migrates to the
full discutions and for critical reading of the
parameter’s limit which is closer to the maximum.
manuscript. This work was supported by CONA-
Asubsequentrunwithnewrangesthatincreasethis
CyT grants 40864 and 42809 and UNAM/DGAPA
limit usually leads to a successfully solution. The
grant IN107202.
price to pay is only an increased number of pro-
file evaluations. But the fast computers make this
REFERENCES
price affordable. The calculation of one generation
of”bugs”takeslessthan2minutesonaPentiumIV Arrieta,A., Stangellini,L., 2004, ESO AstrophysicsSym-
class computer. Usually fifty to a hundred genera- posia, ”Planetary nebulae beyong the Milkey Way”
tionsareneededtofitonelineprofile. Each”bug”in ed.Stangellini,L,J.R.Walsh&N.Douglas(inpress).
Charbonneau, P., 1995, ApJS,101, 309
a generation is independent of the rest of the popu-
Georgiev,L.N., Koenigsberger,G. 2004, A&A,423, 267.
lation,soasimpleparallelizationincreasesthespeed
Groenewegen, M. A. T. & Lamers, H. J. G. L. M. 1989,
ofthecalculationalmostproportionaltothenumber
A&AS,79, 359.
of available processors.
Lamers, H. J. G. L. M., Cerruti-Sola, M., & Perinotto,
The method has however few disadvantages. M. 1987, ApJ, 314, 726.
First of all, the choice of the weights σi in equa- Lamers, H. J. G. L. M., Haser, S., de Koter, A., & Lei-
tion (9) is not automatic. The observed profiles fre- therer, C. 1999, ApJ, 516, 872.
GENETIC ALGORITHM ... 9
Metcalfe, T.S. 2003, ApJ,587L, 43.
Rybicki,G.B., Hummer,D.G., 1978, ApJ, 219, 654.
Sevenster,M.,Saha,P.,Valls-Gabaud,D.,Fux,R.1999,
MNRAS,307, 584.
Sobolev,V.V. 1958, Theoretical Astrophysics, V.A. Am-
bartsumian, Ed. (Pergamon Press, New York).
Steffen,1990, A&A,233, 493.
Teriaca, L., Banerjee, D., Doyle, J.G. 1999, A&A, 349,
636.
LeonidGeorgiev:InstitutodeAstronom´ıa,UniversidadNacionalAut´onomadeM´exico,ApartadoPostal70264,
CD Universitaria,CP 04510 , M´exico DF, M´exico DF ([email protected]).
Xavier Hernandez: Instituto de Astronom´ıa, Universidad Nacional Aut´onoma de M´exico, Apartado Postal
70264, CD Universitaria, CP 04510 , M´exico DF, M´exico DF ([email protected]).