NLTE models of line-driven stellar winds I. Method of calculation and first results for O stars

0.33 MB·English

by  Jiri Krticka

Astronomy&Astrophysicsmanuscriptno.nlte1 February2,2008
(DOI:willbeinsertedbyhandlater)
NLTE models of line-driven stellar winds
I. Method of calculation and first results for O stars
J.Krticˇka1,2,3 andJ.Kuba´t3
1 U´stavteoreticke´fyzikyaastrofyzikyPˇrFMU,CZ-61137Brno,CzechRepublic
4
0 2 DepartmentofPhysicsandAstronomy,UniversityofGlasgow,GlasgowG128QQ,UK
0 3 Astronomicky´ u´stav,AkademieveˇdCˇeske´republiky,CZ-25165Ondˇrejov,CzechRepublic
2
Received1July2003/Accepted24November2003
n
a
J Abstract.New numerical models of line-driven stellar winds of late O stars are presented. Statistical equilibrium (NLTE)

2 equations of themost abundant elements aresolved. Properly obtained occupation numbers areused tocalculate consistent
1 radiativeforceandradiativeheatingterms.Winddensity,velocityandtemperaturearecalculatedasasolutionofmodelhydro-
dynamicalequations.Contrarytootherpublishedmodelsweaccountforamulticomponent windnatureanddonotsimplify

1 thecalculationoftheradiativeforce(e.g.usingforcemultipliers).Wediscusstheconvergencebehaviourofourmodels.The
v abilityof our models topredict correct values of mass-lossratesand terminal velocitiesof selectedlateO stars(mainly gi-
4
antsandsupergiants)isdemonstrated.Thesystematicdifferencebetweenpredictedandobservedterminalvelocitiesreported
9
intheliteraturehasbeenremoved.Moreover,wefoundgoodagreementbetweenthetheoreticalwindmomentum-luminosity
1
relationshipandtheobservedoneforCygOB2supergiants.
1
0
4 Keywords.stars:winds,outflows–stars:mass-loss–stars:early-type–hydrodynamics
0
/
h
1. Introduction librium equations. Despite the abovementionedenormous ef-
p
- fort in wind model calculations, many unsolved problems re-
o NLTEline-drivenstellarwindmodelshavebeencalculatedby
main which may be connectedeither with the simplifyingas-
r variousauthors(cf.deKoteretal.1993,Hillier&Miller1998,
t sumptions mentionedabove or with some others, e.g. neglect
s Vink etal. 2000, Pauldrachet al. 2001, Gra¨feneret al. 2002).
a ofwindnon-stationarity(cf.Runacres&Owocki2002),wind
ThesesophisticatedmodelscoverawiderangeofOBandWR
: 3Dstructure(Dessart&Owocki2003),orwindmagneticfields
v starparameters,anditmayseemthatthereisnoneedforany
i (ud-Doula&Owocki2002).Amongtheseunsolvedproblems
X additionalNLTEwindmodels.However,webelievethisisnot
is the long-standingmomentumproblemof WR stars (forthe
thecase.Firstly,thesemodelsemploydifferentsimplifyingas-
r latestprogressseeGra¨feneretal.2002)orproblemswithfitting
a sumptions which cannot be necessarily fulfilled in all model
ofX-raylineprofiles(Krameretal.2003).
stars. For example,windvelocityis oftencalculatedusingei-
Many of the simplifying assumptions employed in mod-
therasimplifiedβ-velocitylaw(Vinketal.2001,Gra¨feneret
ern wind codes are probably acceptable for high density
al. 2002) or using simplified line-forcemultipliers(Pauldrach
radiatively-drivenstellarwindsofOstars,buttheymaybreak
et al. 2001, Kudritzki 2002) of Castor et al. (1975) and Puls
down for low density stellar winds. Whereas the bulk wind
et al. (2000). Some of these models (e.g., Vink et al. 2001)
material consists of hydrogen and helium, line-driven stellar
donotallowindependentcalculationofobservedterminalve-
windsareacceleratedmainlybytheabsorptionofradiationin
locities.Moreover,temperaturestratificationisoftenneglected
the resonancelines of trace elementslike C, N, O, or Fe. For
(cf. Kudritzki 2002) assuming that the radiative force is not
low-density stellar winds this multicomponentnature may be
significantly influenced by temperature (see Pauldrach 1987)
important for the wind structure itself, influencing both tem-
or calculated using a simplified expression (Pauldrach et al.
peraturestratification(Springmann&Pauldrach1992,Krticˇka
1994, Vink et al. 1999). Although the calculation of correct
& Kuba´t 2001, hereafter KKII) and wind dynamics (Babel
temperaturestratificationisincludedinmodernwindcodes(cf.
1995, Porter & Skouza 1999, KKII, Owocki & Puls 2002,
Pauldrachet al. 2001), the latest publishedgrid of NLTE cal-
Krticˇka & Kuba´t2002). Clearly, such a complicatedsituation
culationsofwindtemperaturedatesbacktothepaperofDrew
cannotbeadequatelydescribedbymodelswithatemperature-
(1989),whousedarelativelyconstrainedsetofstatisticalequi-
independentradiativeforce,orwitharadiativeforceartificially
Sendoffprintrequeststo:J.Krticˇka, dependingonwindtemperatureviathethermalvelocityashas
e-mail:[email protected] been done, e.g., by KKII. Note that the artificial dependence
2 J.KrticˇkaandJ.Kuba´t:NLTEmodelsofline-drivenstellarwindsI.
on the thermal velocity may be removed using an alternative hydrogen and helium if Population I/II stars are consid-
descriptionoftheradiativeforcebyGayley(1995). ered).
Thesituationisespeciallyappealingnow,whenwindmod- – Leveloccupationnumbersofmodelionsareobtainedfrom
elssuitablydescribingtheatmospheresoffirstgenerationstars NLTErateequationswithaninclusionofasuperlevelcon-
havebeencalculated(Kudritzki2002),becausesomeofthese cept(Anderson1989).
modelsarealsoinfluencedbymulticomponenteffects(Krticˇka – Radiation transfer is split into two parts. Line radiative
et al. 2003, see also Kudritzki 2002). Moreover, predicted transfer is solved in a Sobolev approximation (Rybicki
mass-loss rates of low-metallicity winds of some SMC stars & Hummer 1978) neglecting continuum opacity sources,
are more than one magnitude higher than the observed ones lineoverlaps,andmultiplescattering.Continuumradiative
(Bouret et al. 2003). These discrepancies may be connected transferisformallysolvedusingtheFeautriermethodwith
withproblemswithsourcefunctiongradients(Owocki&Puls an inclusion of all free-free and bound-free transitions of
1999). modelions, howeverneglecting line transitions. Thus, we
Althoughsatisfactorydeterminationofmass-lossratesand neglectline-blocking.
terminalvelocitiesofOstarscanbeperformedfromobserva- – For the solution of NLTE rate equationstogetherwith the
tionsindifferentwavelengthregions,thesituationintheBstar radiative transfer equation we apply the Newton-Raphson
domainisnotsoclear.EvenformassiveBstarsthereisasig- method for the radiative transfer in lines and lambda iter-
nificant discrepancy between observations in different wave- ationsin continuum.The simplificationofa lambdaitera-
lengthregions(Vinketal.2000).Evenworse,formostmain- tion may, of course, fail for optically thick continua, like
sequenceBstarsthereisnotrustworthydeterminationofeither the HeII continuum (however, see Sects. 2.3 and 3 for a
mass-loss rates or terminalvelocities, and, moreover,it is not discussionofthisissue).
clearwhichBstarspossesastellarwind(seeBabel(1996)for – Derivedoccupationnumbersareusedtocalculateradiative
thediscussionofthisproblem).Thisisachallengeforsophis- forceinthehydrodynamicequations.
ticatedwindmodelstohelpsupplementmissingobservational – The radiative cooling/heating term is obtained by adding
data. There exist severalsubclasses of B stars whose peculiar thecontributionsfromallfree-free,bound-freeandbound-
behaviour may be connected with their stellar wind. This is bound transitions of model ions. For the calculation of
epecially the case of Be stars, because the plausible explana- this term we use the thermal balance of electrons method
tion of the Be phenomenon is still missing (excluding some (Kuba´tetal.1999).
binaries, for which a binary hypothesis of Kˇr´ızˇ & Harmanec – Thedescriptionofatransitionregionbetweenaquasi-static
(1975)offersatrustworthyexplanation). stellar atmosphere and supersonic stellar wind is simpli-
Last but not least, there is an ambitious idea of the fied.Thisregionisimportantforthecontinuumformation
windmomentum-luminosityrelationship(seeKudritzki&Puls (Pauldrach 1987). This subsonic region was not included
(2000)andreferencestherein)whichmaybeusedasaninde- in our previous models since we did not solve the radia-
pendentreliablemethodfordeterminationofstellardistances. tivetransferequation(seeKKII).Toobtaincorrectionizing
Thus, several types of detailed wind models are necessary to fluxweincludedthissubsonicregioninourmodels,i.e.the
verifyalloftheassumptionsofline-drivenwindmodels. boundaryvelocity at the stellar surface is set to some low
Motivatedbytheseconsiderationsweofferanewindepen- value well below the sonic speed. The particular value of
dent type of NLTE line-driven wind model to obtain correct the base velocity is selected in such a way that lowering
values of radiative force, mass-loss rates and terminalveloci- of this value by a factor of 10 does notchange the model
tiesinthecaseoflow-densitystellarwinds. structure(i.e.mass-lossrateandterminalvelocity)bymore
than10%.Formostofthewindmodelsthismeansthatthe
boundaryvelocityisofthe orderof10−3 ofthe hydrogen
2. Modelassumptions
thermalspeed.
– The calculation of models below and above the critical
Hydrodynamicmodelequationssolvedinthispaperareessen-
pointissplitintwoparts,thusweneglecttheinfluenceof
tiallythesameasinKKII.Aninterestedreaderwillfindmore
the radiationemitted in the regionabovethe criticalpoint
details aboutthe solutionof multicomponenthydrodynamical
on the regionbelow the critical point. However,the influ-
equationsin KKIIandKrticˇka(2003). Modelspresentedhere
ence of the subcritical parts on the supercritical ones is
differ in the method of calculation of the radiative force and
takenintoaccountcorrectly.
the radiative heating/cooling term, and these new models as-
– We neglect Gayley-Owocki heating (Gayley & Owocki
sume slightly different boundary velocities. In the following
1994) which is not important for stars considered in this
sectionswesummarizeourmodelequations.Thebasicmodel
paper(seeKKII).
assumptionsarethefollowing:
– We also neglect the X-ray radiation which may influence
– Weassumeastationaryspherically-symmetricflow. theionizationbalance.
– Wesolvethecontinuityequation,themomentumequation
and the energy equation for each component of the flow,
namely for absorbing ions, nonabsorbing ions (hydrogen
andhelium),andelectrons.Nonabsorbingionsareionsfor
whichtheeffectofaradiativeforcemaybeneglected(e.g., TheseassumptionsarediscussedinmoredetailintheSect.5.
J.KrticˇkaandJ.Kuba´t:NLTEmodelsofline-drivenstellarwindsI. 3
2.1. Equationsofstatistical equilibrium Table 1. Atoms and their ionization stages included in the
NLTE calculations. In this table, level means either an indi-
Ourtreatmentofstatisticalequilibriumequationsisverysimi-
viduallevelorasetoflevelsmergedintoasuperlevel.
lartothatofPauldrach(1987).We solverateequationsinthe
formof(Mihalas1978)
Ion Levels Ion Levels Ion Levels
N P −N P =0, (1) HI 9 NeIV 12 SIII 10
j ji i ij HII 1 NeV 17 SIV 18
j6=i j6=i
X X HeI 14 NeVI 1 SV 14
where N is the number density of ions in state i relative to HeII 14 NaII 13 SVI 1
i
the totalnumberdensityofgivenatom,andP ,P are rates HeIII 1 NaIII 14 ArIII 25
ij ji
of all processes that transfer an atom from a given state i CII 14 NaIV 18 ArIV 19
CIII 23 NaV 16 ArV 16
to state j and back, respectively. They are sums of radiative
CIV 25 NaVI 1 ArVI 1
andcollisionalrates(bothexcitation/de-excitationandioniza-
CV 1 MgIII 14 CaII 16
tion/recombination).Thesystemofrateequations(1)isclosed
NII 14 MgIV 14 CaIII 14
by the equation of normalisation of relative number densities
NIII 32 MgV 13 CaIV 20
foreachmodelatom(i.e.abundanceequation) NIV 23 MgVI 1 CaV 22
NV 13 AlII 16 CaVI 1
Ni =1. (2) NVI 1 AlIII 14 FeIII 29
Xi OII 50 AlIV 14 FeIV 32
This equation can be used instead of any equation from the OIII 29 AlV 16 FeV 30
OIV 39 AlVI 1 FeVI 27
set (1). However, it is computationally convenient to replace
OV 14 SiII 12 FeVII 1
therateequationforthelevelwiththehighestrelativenumber
OVI 20 SiIII 12 NiIII 36
density.
OVII 1 SiIV 13 NiIV 38
WeusetherelativenumberdensityNi insteadofthecom- NeII 15 SiV 1 NiV 48
monly used number density to be consistent with the hydro- NeIII 14 SII 14 NiVI 1
dynamicalpartofthecode,becauserelativequantitiesaread-
vantageousfromthenumericalpointofview.Notethatrelative
whereα isthecorrespondingphotoionizationcross-section,
numberdensitieswereusedforthesolutionoftheequationsof i,ν
J is the meancontinuumintensity (line absorptionspectrum
statisticalequilibriumbyLucy(2001). ν
isneglectedincalculationsofbound-freerates),ν istheion-
i
izationfrequencyofleveli,T istheelectrontemperature,and
e
2.1.1. Modelatomsandrates
anasteriskdenotesLTEvalues.
Upward (C ) and downward (C ) collisional rates (for
RadiativeexcitationratesR aregivenby ij ji
ij
bothexcitationandionization)arecalculatedusingtheexpres-
N R =N B J¯ , (3) sions
i ij i ij ij
whereJ¯ = ∞φ (ν)J dν isprofile-weighedlineintensity N C =N n Ω (T ), (8a)
ij 0 ij ν i ij i e ij e
and ∗
R 4π2e2 N C =N n Ni Ω (T ), (8b)
Bij = hνijmecfij, (4) j ji j e(cid:18)Nj(cid:19) ij e
where νij is line frequency, fij is corresponding oscillator where Ωij(Te) is the integratedcollisional cross-section (col-
strength,eiselementarychargeandme iselectronmass. lisionstrength)andne istheelectronnumberdensity.Finally,
Similarly,theradiativede-excitationratesRji are ratecoefficientsPij forbothexcitationandionizationaresums
ofradiativeandcollisionalterms
N R =N (A +B J¯ ), (5)
j ji j ji ji ij
P =R +C . (9)
ij ij ij
wheretheremainingEinsteincoefficientsaregivenby
WesolveNLTErateequations(1)foreachoftheelements
A = 2hνi3jB , B = giB , (6) listed in Table 1 separately. This allows us (in addition to a
ji c2 ji ji g ij fixed temperature) to fix the number of free electrons during
j
thesolutionoftheequationsofstatisticalequilibrium.Thisstep
andg ,g arestatisticalweightsofinvolvedlevels.
i j greatly reduces the computationalneeds, since instead of one
Radiativeionizationandrecombinationratesaregivenby
large rate matrix we have to deal with a moderate number of
∞ α smallermatricesfortheelementsconsidered.
N R =4πN i,νJ dν, (7a)
i ik i hν ν ThemodelatomsarepredominantlytakenfromTLUSTY
Zνi models(Hubeny1988,Hubeny&Lanz1992,Hubeny&Lanz
N ∗ ∞ α 2hν3
NkRki =4πNk Ni hiν,ν c2 +Jν e−khTνe dν, 1995,Lanz&Hubeny2003).Thesetofmodelatomsisslightly
(cid:18) k(cid:19) Zνi (cid:20) (cid:21) modifiedforthecaseofhotstarwindmodelling.Forthispur-
(7b)
poseweuseddatafromtheOpacityProject(Seaton1987,Luo
4 J.KrticˇkaandJ.Kuba´t:NLTEmodelsofline-drivenstellarwindsI.
&Pradhan1989,Sawey&Berrington1992,Seatonetal.1992, 2.2.1. Lineradiative transfer
Butler et al. 1993, Nahar & Pradhan 1993) and from the Iron
TheaveragedlineintensityiscalculatedusingtheSobolevap-
Project(Hummeretal.1993,Bautista1996,Nahar&Pradhan
proximation(Rybicki&Hummer1978)
1996,Zhang1996,Bautista&Pradhan1997,Zhang&Pradhan
1997,Chen&Pradhan1999). J¯ =(1−β)S +β I , (14)
ij ij c c
We use a similar strategy for the preparation of model
atoms that enter the calculations as was already applied in whereescapeprobabilityfunctionsβ andβ aregivenbyinte-
c
TLUSTY. With the exception of Fe and Ni, levels with low grals
energies (relative to the ground level of a given ionization
stage) are treated individually (neglecting the hyperfine split- β = 1 1dµ1−e−τµ, (15)
c
ting),whereashigherlevelswithsimilarpropertiesaremerged 2Zµ∗ τµ
into superlevels. For Fe and Ni all levels (except the ground 1 1 1−e−τµ
ones) are merged into superlevels. More details on the con- β = dµ , (16)
2 τ
struction of model atoms and atomic data used are given in Z−1 µ
AppendixA. µ = 1−R2/r2 1/2, and I = 4H . The Sobolev optical
∗ ∗ c c
Procedures for the calculation of Ω (T ) are taken from
ij e depthτ isgivenby(Castor1974,Rybicki&Hummer1978)
µ
the TLUSTYcode.Thepartitionfunctionapproximationsare (cid:0) (cid:1)
mostly from Smith & Dworetsky (1988) and Smith (1992). πe2 n n r
τ = i − j g f , (17)
Moreover,thepartitionfunctionapproximationstakenfromthe µ m ν g g i ijv (1+σµ2)
e ij (cid:18) i j(cid:19) ri
TLUSTYcodeareappliedforC, N, O,FeVI,andFeVII.
wherev istheradialvelocityofabsorbingions,n ,n arethe
ri i j
numberdensitiesofindividualstates,
2.1.2. Superlevels
n =N z n (18)
i i atom H
ComparedtothepioneeringworkofAnderson(1989),weusea
simplifiedmergingofindividuallevelsintosuperlevels.Given (similarly for n ), where z = n /n is the number
j atom atom H
thesublevelslwithenergiesE ,relativenumberdensitiesN , densityofagivenatomrelativetothehydrogennumberdensity
il il
andstatisticalweightsgil,thetotalpopulationofsuperlevelis nH,
ρ
givenby n = p , (19)
H
m +z m
H He He
N = N , (10)
i il
whereρ isthedensityofapassivecomponent(hydrogenand
l p
X helium).ThevariableσwasintroducedbyCastor(1974),
andtheenergyofsuperlevelis
dlnv
σ = ri −1. (20)
E = lgilEil. (11) dlnr
i
g
P l il A suitable method for the calculation of escape probability
functions β and β was described by KKII. Finally, the line
Thus, we effectively neglecPt the exp(−E /kT ) term. c
il e sourcefunctionis
Similarly,the photoionizationcrosssection of superlevelα
i,ν
iscalculatedby S = NjAji . (21)
α = lgilαil,ν, (12) ij NiBij −NjBji
i,ν
g
P l il
2.2.2. Continuumradiative transfer
andtherelativenumberdensityPofasublevelofagivensuper-
levelisapproximatedas For the calculation of the continuum intensity we use a
procedure for the solution of the transfer equation based
N = gilNi . (13) on the Feautrier method in the spherical coordinates (for a
il
lgil description, see Mihalas & Hummer 1974 or Kuba´t 1993).
However, for the solution of the obtained tridiagonal system
NumericaltestsshowedthatthePneglectoftheexp(−Eil/kTe) of equations we use the numerical package LAPACK
termdoesnothaveasignificanteffectonthecalculatedradia- (http://www.cs.colorado.edu/˜lapack,
tiveforce. Anderson et al. 1999) instead of the classical Gaussian
scheme. All bound-free transitions of explicit ions, free-free
transitionsofhydrogenandheliumandlightscatteringdueto
2.2. Radiativetransfer
free electrons contribute to the source function.The radiative
Radiative flux at the base of the stellar wind H is taken transferequationissolvedonlyabovethehydrogenionization
c
from H-He spherically symmetric NLTE model atmospheres edge.Forfrequencieslowerthanthehydrogenionizationlimit
of Kuba´t (2003). The line transfer equation and continuum we use an optically thin approximation J (r) = 4W(r)H ,
ν c
transferequationaresolvedseparately. whereW(r)= 1(1−µ )isthedilutionfactor.
2 ∗
J.KrticˇkaandJ.Kuba´t:NLTEmodelsofline-drivenstellarwindsI. 5
2.3. SolutionofNLTEequations of vector h), since the level populations depend on the solu-
tionofradiativetransferequationincontinuum.However,itis
For the solution of the NLTE rate equations (1), (2) together
possible to obtainapproximatevaluesof these derivativesne-
with the radiative transfer equation we use a combination of
glectingtheexplicitdependenceofratesP onthecontinuum
lambdaiterationsandofacompletelinearization.Forthecon- ij
intensity.DerivingEq.(23)withrespecttohandneglectingthe
tinuumtransferweusesimplelambdaiterations.Althoughthe
dependenceonJ (N,h)weobtainasystemofequationsfor
stellar wind is optically thick near to the star below the HeII ν
∂N /∂h
ionization limit, these iterations converge relatively well (see i
also Sect. 3). In lines, we use complete linearization of the
∂P ∂P ∂N ∂P ∂N
radiativetransferequationina simpleSobolevapproximation jiN + ji iN + jiN +P j −
j j j ji
∂h ∂N ∂h ∂N ∂h
(neglectingcontinuumradiation,see Eq.(14)).Thelineinten- j6=i(cid:20) i (cid:18) j (cid:19) (cid:21)
X
sityJ¯ijdependsonlyonthelevelpopulationsofthegivenatom ∂Pij ∂Pij ∂Nj ∂Pij ∂Ni
N + N + N +P
forthiscase.Thus,thestatisticalequilibriumequationsfordif- i i i ij
∂h ∂N ∂h ∂N ∂h
ferentatomsareindependent(forafixedhydrodynamicalstruc- Xj6=i(cid:20) j (cid:18) i (cid:19) (cid:21)
ture).Hence,thefollowingschememaybeobtainedfromsta- =0. (24)
tisticalequilibriumequations(1),(2)forthecalculationofthe
All calculated derivatives ∂N /∂h are inserted into the
correctionδN of relativenumberdensitiesin thegivenitera- i
i
Newton-Raphson matrix of hydrodynamic equations (corre-
tionstep
spondingto theradiativeforce,radiativeheating/coolingterm
∂P ∂P andionizationcharge,seeKKIIandKrticˇka2001).
N jiδN + P +N ji δN −
j i ji j j Note that the system of equationsfor derivatives (24) has
∂N ∂N
Xj6=i(cid:20) i (cid:18) j(cid:19) (cid:21) the same form as the system of NLTE rate equations (22),
− P +N ∂Pij δN +N ∂PijδN =0, (22) however with a different right-hand side. Thus, because the
ij i i i j
∂N ∂N LAPACK packageusesanLU decompositionforthesolution
j6=i(cid:20)(cid:18) i(cid:19) j (cid:21)
X of the system of linear equations, the values of ∂N /∂h can
i
where the terms ∂Pij/∂Nj are calculated using Eqs. (14) – be calculated using data from the previous solution of NLTE
(21). For details see Appendix B. NLTE equations (22) are equations(22).
solved for each element separately. Note that Eq. (22) is an This method is especially useful in the case when the in-
application of a Newton-Raphson method to statistical equi- tensityofradiationinthecontinuumdoesnotdependonlocal
libriumequations(1).ForthesolutionofthesystemofNLTE properties(i.e.,foranopticallythincontinuum,whenitisgiven
equations(22) we use the numericalpackageLAPACK. Note by the intensity emerging from the underlying photosphere).
that we can again take full advantage of LAPACK subrou- From the mathematical point of view, this method splits the
tines designed for general band matrices because our NLTE Jacobi matrix into several parts and it is equivalent to a lin-
rateequationmatriceshaveabandstructure. earizationofthewholesystemofmodelequations(i.e.hydro-
dynamics, NLTE equations and radiative transfer). However,
from the numerical point of view, it saves a huge amount of
2.4. Derivatives ofoccupationnumbers
computertimebecauseitisnotnecessarytoinvertlargematri-
The hydrodynamic equations are solved using a Newton- ces.Asimilarmethodhasbeenused,e.g.,byAnderson(1985),
Raphson method. To obtain a well convergedmodel we have and Kuba´t (1994, 1997) referred to it as to the implicit lin-
toknowtheexplicitvaluesofderivativesoftheradiativeforce earizationoftheb-factors.
and of the radiativeheating/coolingterm with respectto fluid
variables(i.e.densities,temperatures,etc.).However,because
2.5. Calculationofthe radiative force
the radiative force and the radiative heating/cooling term de-
pendon levelpopulations,we have to calculate derivativesof TheradiativeforceiscalculatedusingtheSobolevapproxima-
levelpopulationswith respectto fluidvariables.Thesecanbe tiondirectlyasthesumofcontributionsofindividuallinesafter
calculatedinthe followingway.Forclarity,we rewriteNLTE Castor(1974),
rateequations(1)inasymbolicform
8π v 1
grad = r νH µ 1+σµ2 1−e−τµ , (25)
Pji(N,Jν(N,h),h)Nj− i ρic2 r lines cZµc
j6=i X (cid:0) (cid:1)(cid:0) (cid:1)
X
whereτ isgivenbyEq.(17).OccupationnumbersN calcu-
−N P (N,J (N,h),h)=0, (23) µ i
i ij ν latedduringtheNLTEiterationstepareusedforacalculation
Xj6=i oftheradiativeforce.NotethatinEq.(17)numberdensitiesof
wherewe denotedallexplicitdependenciesof rate equations. individuallevelsareusedinsteadofthetotalnumberdensityof
The elements of vector N are level populations N , the ele- a superlevel. These number densities may be calculated from
i
mentsofvectorharefluidvariablesh=(n ,ρ ,T ,v ),and Eq. (13) for those levels which are explicitly included in the
e p e ri
the meancontinuumintensityJ dependonh and N via the NLTE rate equations. However,non-explicitlevels, i.e., some
ν
radiativetransfer equation.Clearly, it is notfeasible to obtain ofthehigherexcitedlevelsorexcitedlevelsofthehighestion
correctexplicitvalues of ∂N /∂h (where h denotes elements ofanatom,arealsoincludedinthecalculationoftheradiative
i
6 J.KrticˇkaandJ.Kuba´t:NLTEmodelsofline-drivenstellarwindsI.
force.Forthesenon-explicitlevelsweuseeithertheBoltzmann The calculation of a given model below the critical point
equilibriumat the wind temperaturein the case of metastable consistsofthreenestediterationcycles(seeFig.1).Inthein-
levelsortheradiativeexcitationequilibriumofAbbott&Lucy nermost iteration cycle the NLTE rate equations (1) together
(1985)forotherlevels.Theinfluenceofthesenon-explicitlev- with the radiative transfer equation both in lines (14) and in
elsontheradiativeforceisverysmall,however. continuum are solved. For this innermost iteration cycle the
The term ∂grad/∂(dv /dr) is calculated explicitly after fluidvariables(i.e.densities,velocitiesandtemperaturesofin-
i ri
(25) where the variable σ depends on a velocity gradient via dividualwindcomponents)arekeptfixed.
(20).Thistermisincludedinthecriticalpointcondition(gen- Inthemiddlecycleweiteratethehydrodynamicalstructure
eralised CAK condition, KKII, Eq. (57) therein). Derivatives ofthemodelforfixedlevelpopulations.Thisiterationcycleis
of numberdensitiescalculated by meansof Eq. (24) are used essentially the same asthe iterationcycle ofthree-component
to obtain derivatives of the radiative force with respect to nonisothermalwindmodelsdescribedbyKKII.To savecom-
fluid variables(densities,velocities,andtemperatures)forthe putingtime,forlargerelativechangesoffluidvariables(greater
Newton-Raphsoniterations. than 10−2) we perform only 3 iterations of NLTE rate equa-
Oscillatorstrengthsnecessaryforthecalculationofthera- tions,butwhentherelativechangesoffluidvariablesarelower,
diativeforceareextractedfromtheVALDdatabase(Piskunov and we may expect that we are closer to the correct solution,
et al. 1995, Kupka et al. 1999). This set of data consists of weperformmoreiterationsofNLTErateequations(upto50).
about180000linetransitionsinthewavelengthinterval250– Weassumethatmodeliswellconvergediftherelativechanges
10000A˚. ofbothfluidvariablesandoccupationnumbersarelowerthan
10−5.
In the outermost iteration cycle we solve for the value of
2.6. Radiativecoolingandheating
the base density for which a model smoothly passes through
thecriticalpoint.Again,theseiterationswerealreadydescribed
The radiative cooling and/or heating term is calculated using
by KKII.Forthe base density(whichdeterminesthe valueof
thethermalbalanceofelectrons(Kuba´tetal.1999).Forthecal-
mass-lossratedM/dt)anaccuracyof1%issufficient.
culationofthistermweincludealltheexplicitradiativebound-
After the base wind density is known, we calculate wind
freeandcollisional(bothbound-boundandbound-free)transi-
model above the critical point by a sequence of iterations of
tionsandfree-freetransitionsofhydrogenandhelium.Again,
fluid variables.Each of these iterations is precededby the set
inordertoobtainawell-convergedmodel,derivativesofnum-
of NLTE iterations. Boundary conditions of the model above
berdensitiescalculatedbymeansofEq.(24)areusedtoobtain
thecriticalpointarespecifiedatthecriticalpointandaretaken
derivatives of the radiative cooling and/or heating terms with
fromthemodelbelowthecriticalpoint.Themodelbelowthe
respecttofluidvariables.Thecalculatedderivativesofthe ra-
critical point is fixed during the solution above the critical
diativecooling/heatingtermareinsertedintotheJacobimatrix
point.
of the Newton-Raphson iteration of the model hydrodynamic
Splitting of a wind solution into two parts may cause
structure. Detailed expressions for the derivatives of heating
problems since the radiation emitted from the outer region
andcoolingratesarepresentedinAppendixC.
(i.e. above the CAK point) may influence the inner region.
However,testcalculationswherewechangedtheouterbound-
2.7. Electron density ary of inner solution showed that this is not the case, i.e. that
the inner wind solution dependsonly marginallyon the outer
The charge of the passive component is calculated from ion-
solution.Thereasonforsuchabehaviourisprobablythestrong
izationfractionsofbothhydrogenandhelium.Thisenablesus
dependenceofthewindcontinuumopticaldepthonfrequency.
to calculate the correctvalueof the electron density,whichis
Calculations show that the stellar wind above the CAK point
importantforthepropersolutionofNLTErateequations.
is essentially optically thin for continuum radiation with fre-
quenciesbelowtheHeII ionizationlimit,thusradiationisnot
2.8. Iteration scheme forthe modelcalculation absorbedinthecontinuumabovetheCAKpointanddoesnot
influence the stellar wind belowthe CAK point. On the other
TheinitialestimateoffluidvariablesfortheNewton-Raphson hand,stellarwindabovetheCAKpointisopticallyverythick
iterations(namelytheso-calledbeta-lawforvelocitiesandden- forradiationwithfrequenciesabovetheHeII ionizationlimit.
sity obtained from the continuity equation) was described by Again,duetothislargeopticalthicknesstheradiationemitted
KKII.Fortheinitialvalueoftheintensityweuseanoptically intheregionsabovetheCAKpointdoesnotinfluencetheinner
thin approximation, i.e., the radiation field is defined by the stellar wind.Theinfluenceoflinesabovethecriticalpointon
lowerboundaryconditionatthebeginningofourcalculation. thesubcriticalregionhasbeenalsoneglected.
The computational procedure involves two similar basic
steps,namelysolutionbelowandabovethewindcriticalpoint.
3. Convergencebehaviourofmodels
This splitting into two steps is necessary and is related to the
critical point condition. In the first step the value of the base Asanexamplewestudytheconvergencebehaviourofamodel
density, which allows the model to smoothly pass the critical correspondingtoαCam(HD30614)withparametersgivenin
point, is obtained. Afterwards, the model above this critical Table2.InFig.2weplottedrelativechangeofpopulationnum-
pointiscalculated. bers, and hydrodynamical variables (i.e. densities, velocities,
J.KrticˇkaandJ.Kuba´t:NLTEmodelsofline-drivenstellarwindsI. 7
M, R, T , z
eff
model below the critical point model above the critical point
NLTE equations NLTE equations
radiative transfer radiative transfer
convergency? convergency?
hydrodynamics hydrodynamics
convergency? convergency?
model above
critical the critical point
condition?
output model:
radial variation of densities,
model below
velocities, temperatures and
the critical point
level populations
Fig.1.Iterationschemeforthemodelcalculation
temperaturesand chargesof individualwind components)for weconcludethattheoscillatorybehaviouriscausedbyweakly
themodelbelowthecriticalpointandforthemodelabovethe populated levels and that their influence on the convergence
criticalpoint. behaviourofthehydrodynamicalstructureisnegligible(upper
rightpanelofFig.2).Thisoscillatorybehaviourisnotatypical
First we discuss the convergence behaviour of the model featureofconvergenceofourmodels.
belowthecriticalpoint.Fortheplotofrelativechangeoflevel
populations (upper left panel of Fig. 2) we combine NLTE The plot of hydrodynamicalstructure convergence(upper
iterations performed during individual hydrodynamical itera- right panel of Fig. 2) shows very fast convergence which is
tionstepsforfixedbasewinddensity(correspondingtothefi- slightly slowed down by the oscillatory behaviour aroundthe
nalCAK solution).Thus,inthisplotthepeakscorresponding relativechange10−3.Thisiscausedbythechangeofthemax-
to the hydrodynamicaliteration steps occur. They are caused imumallowednumberofNLTEiterationsteps(seeSect.2.8).
by a reaction of NLTE populations to a significant change of Althoughduetothischangetheconvergenceisslightlyslowed
the fluid variablesafter the calculationof the new estimate of down,thetotalcomputingtimeislessbecausewedonotlose
the hydrodynamic structure. For the relative change of about time with calculation of precise occupation numbers for hy-
10−5 relative changesoscillate. However,if we plot only rel- drodynamical structure which is far from the converged one.
ative changesof populationsof levelswith largestoccupation Finally,notethatthe relativelysmallinitialrelativechangeof
numbers(levelswithN >0.01),nooscillationsappear.Thus, fluid variables (of order 10−2) of the hydrodynamical struc-
i
8 J.KrticˇkaandJ.Kuba´t:NLTEmodelsofline-drivenstellarwindsI.
-1
1e-01 e 10
all levels ur
n major levels ct
o u
opulati 1e-02 mic str 10-2
of level p 1e-03 ydrodyna 10-3
mum relative change 111eee---000654 m relative change of h 1100--54
xi u
a m
m xi
1e-07 ma 10-6
0 100 200 300 400 500 600 0 2 4 6 8 10 12 14
iteration number iteration number
1
1 e 10
ur
n ct
o -1 u
mum relative change of level populati 1111100000-----65432 m relative change of hydrodynamic str 11110000-----154321
xi 10 u 10
a m
m xi
10-7 ma 10-6
0 100 200 300 400 500 600 700 0 5 10 15 20 25 30 35 40
iteration number iteration number
Fig.2.Theconvergenceofthe modelsbelow(upperpanels)and above(lowerpanels)the criticalpoint.Upperleft panel:The
convergenceof levelpopulationsbelow the critical point. Note that NLTE iterations in consecutivehydrodynamicaliterations
arecombinedinthisplot.Thisexplainstherelativelylargenumberofiterationsandpeaksthatoccurafterthehydrodynamical
iteration step. The oscillating behaviour for the relative change of order 10−5 is due to the weakly populated levels, because
therearenooscillationsintheplotofconvergenceofrelativechangeofmajorlevels(i.e.levelswithN > 0.01).Upperright
i
panel:Theconvergenceofhydrodynamicalstructure(maximumrelativechangeofallhydrodynamicalvariables,i.e.densities,
velocities,temperaturesandchargesofindividualwindcomponents)ofthemodelbelowthecriticalpoint.Theslightlyoscillating
behaviouraroundtherelativechange10−3 iscausedbythechangeofmaximumallowednumberofNLTEiterationsteps(see
thetext).Lowerleftpanel:Theconvergenceoflevelpopulationsofthemodelabovethecriticalpoint.Again,NLTEiterations
in differenthydrodynamicaliterations are put together and peaks occur after the hydrodynamicaliteration step. Note the fast
convergence of NLTE population between successive hydrodynamicaliterations. Lower right panel: The convergence of the
hydrodynamicalstructureofthemodelabovethecriticalpoint.
ture is caused by the way the critical point is calculated, ie. namicstructureyieldssomethinglikea”forest”ofiterationsin
that models with similar base densities are calculated. Thus, theplot,whichiscausedbyalargechangeofoccupationnum-
the initial change of fluid variables of the final model below bers after a new hydrodynamicstructure was calculated. The
the criticalpointis verysmall. The convergencebehaviourof convergenceof a hydrodynamicalstructure(lowerrightpanel
previous models (i.e. models which do not correspond to the ofFig.2)ishoweverslowerthanbelowthecriticalpoint.This
criticalsolution)issimilar,howevertheinitialchangeoffluid ispartiallycausedbylowermaximumchangeofwindparam-
variablesishigher. etersallowedduringNewton-Raphsoniterationstep.Thiscon-
straintofthemaximumchangeofwindparametersisnecessary
Theconvergencebehaviourofthemodelabovethecritical forthestabilityofthesolution(seeKrticˇka2003).
point(lowerpanelsofFig.2)issimilartothatbelowthecriti-
calpoint.Theconvergenceofoccupationnumbers(foragiven The convergenceof the model temperature is very fast in
hydrodynamicalstructure)isfasterbecausethestellarwindis thecaseofamodelbelowthecriticalpoint.Therelativechange
opticallythinner(notethatthisfastconvergenceofNLTEiter- of temperature is typically ten times lower than the relative
ationscombinedwith larger numberof iterations of hydrody- changeoftherestofthehydrodynamicalvariablesinthiscase.
J.KrticˇkaandJ.Kuba´t:NLTEmodelsofline-drivenstellarwindsI. 9
4.1. Terminalvelocities
2600
2400 Wind terminal velocities of galactic O stars can be measured
with a typical uncertainty of 100kms−1. Thus, in principle,
ved) 2200 they may provide a relatively strict test for any wind model.
er
bs Unfortunately, because terminal velocities depend on surface
L, o 2000 escapevelocities(seeCastoretal.1975)thepreciseprediction
S
L
-1m s] ( 1800 ostfelwlairndmatesrsmainndalravdeiluosc.iCtyornesleieqsueonntlyth,ethperleocnisge-sktannodwinlegdpgreoobf-
[k∞ 1600 lem of“massdiscrepancy”inO stars(cf.Herreroetal. 1992,
v Lanzetal.1996)complicatesastraightforwardcomparisonof
1400 calculatedandobservedterminalvelocities.Thisisespecially
the case of stars HD30614, HD188209, and HD218915, for
1200
1400 1500 1600 1700 1800 1900 2000 2100 2200 2300 2400 2500 whichweusedthesamestellarparameters,butwhichhavedif-
v∞ [km s-1] (calculated) ferentobservedterminalvelocities.
Fig.3. Comparison of calculated and observed (taken from Thecomparisonofobserved(takenfromLSLandBianchi
LSL,BG)terminalvelocities. & Garcia 2002, hereafterBG) and calculatedterminalveloci-
tiesformodelstarsshowninFig.3impliesthatourNLTEwind
modelsareabletopredictrelativelycorrectvaluesofwindter-
3.2 minal velocities. The fit is improved compared to our previ-
ouscalculationspresentedinKKII(seeTable2).Notethatour
3 calculations consistently solved the problem of high theoret-
ical terminal velocities obtained by LSL using the so called
2.8 “cookingformula”ofKudritzkietal. (1989), see also Table2
andFig.12inKKII.However,sinceourmodelsincludeseveral
/v∞esc 2.6 simplifyingassumptions,itisnecessarytotestthisconclusion
v
againstmoreadvancedmodels.
2.4
It remains unclear why LSL obtained too high values of
terminal velocities when theoretical calculations (KKII) with
2.2
the same force multipliers yielded terminal velocities consis-
tent with observations. Note that Pauldrach et al. (2001) also
2
predicted correct terminal velocities, however for a different
30000 31000 32000 33000 34000 35000 36000 37000 38000 39000 40000
T [K] setofstars.
eff
Fig.4.Comparisonofcalculated(crosses)andobserved(taken The calculatedterminalvelocitiesare, in fact, nottoo dif-
fromLSL,BG, errorbars)ratioof the terminalvelocityto the ferent from those calculated by KKII, although the terminal
escapevelocityv /v . velocities presented by KKII were obtained simply using the
∞ esc
force multipliers of Abbott (1982). This relatively surprising
conclusioncanbeattributedtothefactthattheαparametercal-
culated properlyusing line statistics (Puls et al. 2000) is very
close to the valuescalculatedby Abbott(1982). And because
Ontheotherhand,theconvergencebehaviourofamodelabove
theterminalvelocitiesdependontheαvalue,terminalveloci-
thecriticalpointisdominatedbythetemperaturessincetherel-
tiescalculatedinthispaperandthoseusingAbbott’smultipli-
ativechangesoftemperaturearethehighest.
ers (KKII) do not differ significantly. However, properly cal-
culatedk parameters(Pulsetal.2000) aresignificantlylower
thanthosecalculatedbyAbbott(1982).Indeed,ourmass-loss
rates (which are, in fact, influenced by the k parameter) are
4. Comparisonofcalculatedandobserved
systematically lower than those of KKII and are close to the
quantities
observedvalues(seealsothenextsectionforthediscussionof
mass-lossrates).
Because in future we intend to calculate mainly wind models
of B stars, to test the reliability of ourassumptionswe calcu- The stellar wind theory predicts a strong correlation be-
lated NLTEwindmodelsof starswith spectraltyperelatively tweenwindterminalvelocityandstellarescapevelocity.Thus,
close to B, i.e. later O stars. For this purpose we selected O wehavealsocomparedtheratiooftheobservedandcalculated
stars with a spectral type O6 or later from the set of Lamers terminal velocity to the escape velocity v /v (see Fig. 4,
∞ esc
etal.(1995,hereafterLSL).We includedonlystarsforwhich v is calculated using stellar parameters from Table 2). Our
esc
reliabiledeterminationof theirmass-lossrate existsinthe lit- derivedmeanratiov /v ∼2.6isinagoodagreementwith
∞ esc
erature.Stellarparameters(masses,radiiandeffectivetemper- themeanobservedv /v ratio2.5derivedbyLSLforstars
∞ esc
atures)takenfromLSLarelistedinTable2. withaneffectivetemperaturehigherthan25000K.
10 J.KrticˇkaandJ.Kuba´t:NLTEmodelsofline-drivenstellarwindsI.
Table 2. StellarandwindparametersofselectedO starsfromourtest sample.Stellar parameters(columnsM, R , T ) were
∗ eff
adoptedfromLamersetal.(1995).Observedmasslossratesweretakenfrom[1]–Leitherer(1988),[2]–Scuderietal.(1992),
[3] – Lamers& Leitherer(1993), [4] – Puls et al. (1996), [5] – Crowther & Bohannan(1997), [6] – Scuderiet al. (1998), [7]
–Lamersetal.1999,[8]–BG(seethetext),predictedvaluesofM˙ andv aredeterminedusingourcode.Notethatterminal
∞
velocitiescalculatedbyKKIIwereobtainedusingsimplifiedforcemultipliersafterAbbott(1982).
Star HD Sp. Stellarparameters MasslossratesM˙ Terminalvelocitiesv∞
number type M R∗ Teff Hα radio UV model LSL,BG KKII model
[M⊙] [R⊙] [K] [10−6M⊙yr−1] [kms−1]
αCam 30614 O9.5Ia 43.0 27.6 30900 5.2[4];2.9[6] 2.5[6] 2.4 1500±200 1450 1470
34656 O7II 30.0 9.9 38100 0.49[1] 0.54 2100±100 2590 2410
λOriA 36861 O8III 30.0 12.3 36000 0.81[7] 0.66 2200±300 2230 2260
V973Sco 151804 O8Iab 70.0 34.0 34000 12.0[5] 10.0[3] 18.0 1500±200 1370 1700
152424 O9.7Ia 52.0 33.4 30500 3.8[3] 5.5[3] 6.5 1500±100 1350 1540
63Oph 162978 O8III 40.0 16.0 37100 0.83[1] 2.5 2200±200 2090 2160
163758 O8III 50.0 20.1 37200 5.0[8] 5.7 2100±200 1890 2100
HR7589 188209 O9.5Ib 43.0 27.6 30900 3.4[2] 2.4 1600±100 1450 1470
190864 O7III 42.0 14.0 39200 1.5[4] 2.8 2300±200 2340 2370
210809 O9Ib 38.0 21.4 32000 4.0[4] 1.0 2000±200 1700 2200
λCep 210839 O6Iab 51.0 19.6 38200 5.3[4] 2.1[3] 1.5[8] 4.7 2100±200 1930 2150
218915 O9.5Iab 43.0 27.6 30900 1.1[1] 2.4 1800±100 1450 1470
10−4 10−4
♦ Pulsetal.(1996)
(cid:3)+ SSccuuddeerriieettaall..((11999928)) 00))
× Leitherer(1988) 20
△ Lamers&Leitherer(1993) (
˙M−1[Myr](observed)⊙ 1100−−65 ⋆cb LBCaGromwetrhseert×&a⋆l.B(o1♦h9a9n9n)an×(cid:3)++♦×(1♦997)△♦b b△△ △c −1yr](calculated,Vinketal. 1100−−65 ♦♦ ♦ ♦♦♦ ♦♦♦ ♦
⊙
M
[
˙M
10−7 10−7
10−7 10−6 10−5 10−4 10−7 10−6 10−5 10−4
M˙ [M⊙yr−1](calculated) M˙ [M⊙yr−1](calculated,thispaper)
Fig.5.Comparisonofcalculatedandobservedmass-lossrates. Fig.6.Comparisonoftheoreticalmass-lossrates.
O star winds for which unsaturated wind lines originate only
4.2. Mass loss rates
fromtraceelements.However,precisemodellingcanovercome
The credibility of observed mass-loss rates in most stars is thesedifficulties(BG).
much worse than their wind terminal velocities. Although in We concludethatthere is generallya goodagreementbe-
generalthe wind models are able to predictcorrectmass-loss tween predicted and observed mass-loss rates. However, sig-
rates (cf. Vink et al. 2000), there is a large scatter for indi- nificant differenciesexist for individualstars. Similarly, com-
vidual stars. This situation is demonstrated in Fig. 5 where parisonofmass-lossratescalculatedinthispaperandbyVink
wecomparecalculatedandobservedmass-lossrates.Wecom- etal.(2000)showsrelativegoodagreementbetweenthesepre-
piled observedmass-lossrates basedon the Hα measurement dicted mass-loss rates. This possibly supports the conclusion
(Leitherer1988,Scuderietal.1992,Pulsetal.1996,Crowther of Puls (1987) that multiple scattering (included into Vink et
& Bohannan 1997, Scuderi et al. 1998, Lamers et al. 1999) al. (2000) models) does not significantly influence mass-loss
andonthe radioemission(Lamers&Leitherer1993,Scuderi ratesofweakerstellarwinds.
et al. 1998) together with mass-loss determinations obtained An important caveat has to be mentioned. The observed
using UV lines by BG. Generally, UV line-diagnostic is be- mass-loss rates may be influenced by wind clumping (e.g.
lieved to be badly affected by the unknown wind ionization Crowtheretal.(2002),Massaetal.(2003),Bouretetal.(2003))
fractions(Lamersetal.1999).Thisisespeciallythecaseofthe andthus,theobservedmass-lossratesshallbelowerinsucha