Evolution of Line-Force Multiplier Parameters in Radiation Driven Winds of Massive Stars
EEvolution of Line-Force Multiplier Parameters inRadiation Driven Winds of Massive Stars
Alex Camilo Gormaz Matamala
Instituto de F´ısica y Astronom´ıaFacultad de CienciasUniversidad de Valpara´ısoPrograma de Doctorado en Astrof´ısica a r X i v : . [ a s t r o - ph . S R ] F e b A la memoria de mi abuela,Br´ıgida del Rosario Vega Palma.To the memory of my granny,Br´ıgida del Rosario Vega Palma.
This thesis is solely my own composition,except where specifically indicated in the text.Total or partial reproduction, for scientific or academic purposes,is authorised including a bibliographic reference to this document.Alex Camilo Gormaz MatamalaNovember 2019Valpara´ıso, Chile knowledgements “The most terrifying fact about the Universe is not that it is hostile but that it is indifferent;but if we can come to terms with this indifference and accept the challenges of life within theboundaries of death – however mutable [hu]man may be able to make them – our existenceas a species can have genuine meaning and fulfilment. However vast the darkness, we mustsupply our own light.”
Stanley Kubrick. ”All we have to decide is what to do with the time that is given to us.”
Gandalf.The present doctoral work would not have been possible to exist without the great helpand support from many people.Firstly, I would like to express my sincere thanks to my advisor Prof. Michel Cur´e for thecontinuous support of my Ph.D study and related research, for his patience, motivation, andimmense knowledge. His guidance helped me throughout these years of research. I also saythanks to the members of my former commission of my PhD Thesis Project defence, Prof.Radostin Kurtev and Prof. Francisco Najarro, who trusted in the success of this work. ToProf. Lydia Cidale and Roberto Venero, for their helpful support and feedback from beyondthe Andes. I sincerely thank J. Puls for helpful discussions that improved this work andfor having put to our disposal his code FASTWIND. My sincere thanks also goes to Prof.D. John Hillier and Prof. Jose Groh, who provided me the opportunity to work with them asvisitor at University of Pittsburgh (USA) and Trinity College Dublin (Ireland) respectively.And to Prof. Alex Lobel, who accepted me as an intern at the Royal Observatory of Belgium.Besides, I would like to say thanks to each LOC and SOC of the International Workshopon Wolf-Rayet Stars 2015, the Potsdam Astrophysical Summer School 2016, the 10th IAUSymposia 329
Lives and Death-throes of Massive Stars
Ap-plication of radiative transfer to stellar and planetary atmospheres o ontents Aknowledgements 7Abstract 111 Introduction to Massive Stars 13 k, α, δ ) . . . . . . . . . . . . . 383.2 Calculation of the M ( t ) factor . . . . . . . . . . . . . . . . . . . . . . . . . . 413.2.1 Selection of atomic database . . . . . . . . . . . . . . . . . . . . . . . 413.2.2 Ionisation equilibrium . . . . . . . . . . . . . . . . . . . . . . . . . . . 423.2.3 Radiation field . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 453.2.4 Determination of the optical depth range . . . . . . . . . . . . . . . . 463.2.5 Iterative procedure . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 473.2.6 A single-step test . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 483.3 m-CAK Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 483.3.1 Self-consistent calculations . . . . . . . . . . . . . . . . . . . . . . . . 483.3.2 Range of validity for line-force parameters . . . . . . . . . . . . . . . . 5290 CONTENTS W -Lambert equation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 664.2 Lambert-procedure . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 674.2.1 Calculation of line-acceleration g line ( r ) . . . . . . . . . . . . . . . . . . 694.2.2 Solution of equation of momentum . . . . . . . . . . . . . . . . . . . . 714.2.3 Convergence of models . . . . . . . . . . . . . . . . . . . . . . . . . . . 744.3 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 754.3.1 Initial converged Lambert-hydrodynamics . . . . . . . . . . . . . . . . 764.3.2 Mass-loss rate as free parameter . . . . . . . . . . . . . . . . . . . . . 794.4 Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 824.4.1 Fitting observational spectra . . . . . . . . . . . . . . . . . . . . . . . 824.4.2 Comparison with previous self-consistent studies . . . . . . . . . . . . 864.5 Conclusions for Lambert-procedure . . . . . . . . . . . . . . . . . . . . . . . . 90 bstract Massive stars expell strong stellar winds which are described by the theory of radiation-drivenwind. Accurate mass-loss rates are needed to properly describe the stellar evolution acrossthe Hertzsprung-Russel Diagram.We present two self-consistent procedures that couple the hydrodynamics with calculationsof the line-force in the frame of radiation wind theory. These procedures give us the line-force parameters, the velocity field, and the mass-loss rate. The first one is based on theso-called m-CAK theory. Such computations contemplate the contribution to the line-forcemultiplier from more than ∼ ,
000 atomic transitions, an NLTE radiation flux from thephotosphere and a quasi-LTE approximation for the occupational numbers. A full set ofline-force parameters for T eff ≥ ,
000 K and surface gravities higher than 3.4 dex for twodifferent metallicities are presented, along with their corresponding wind parameters (terminalvelocities and mass-loss rates). Here, we find that the already known dependence of line-force parameters on effective temperature is enhanced by the dependence on log g . Terminalvelocities present a steeper scaling relation with respect to the escape velocity, this mightexplain the scatter values observed in the hot side of the bistability jump. For the caseof homogeneous winds (i.e., without clumping) comparison of self-consistent mass-loss ratesshows a good agreement with empirical values. We also consider self-consistent wind solutionsthat are used as input in FASTWIND to calculate synthetic spectra. By comparison withthe observed spectra for three stars with clumped winds, we found that varying the clumpingfactor the synthetic spectra rapidly converge into the neighbourhood region of the solution.Therefore, this self-consistent m-CAK procedure significantly reduces the number of freeparameters needed to obtain a synthetic spectrum.The second procedure (called Lambert-procedure ) provides a self-consistent solution be-yond m-CAK theory and its approximations, and line-acceleration is calculated by the fullNLTE radiative transfer code CMFGEN. Both the mass-loss rate and the clumping factor areset as free parameters, hence their values are obtained by spectral fitting after the respectiveself-consistent hydrodynamics is calculated. Since performing the Lambert-procedure requiressignificant computational power, the analysis is made only for the star ζ -Puppis. It is foundthat fitted wind-parameters are close to those predicted by the m-CAK prescription. Thissuggests that both methodologies providing a lower clumping effect on the wind that those112 CONTENTS suggested by previous authors.We illustrate the future potential of the self-consistent m-CAK prescription, showing thefirst results of two ongoing works: the spectral fitting for a set of high resolution spectraobserved by
Hermes and the development of new evolutionary tracks with the Geneva evolu-tive code using self-consistent mass-loss rates. The promising results gives a positive balanceabout the future applications for the self-consistent solutions presented on this thesis. hapter 1
Introduction to Massive Stars
The study of massive stars (i.e., stars with M ∗ > M (cid:12) ) is a relevant topic in the frame-work of stellar astrophysics, because these stars exhibit some of the most extreme physicalconditions, such as the hottest temperatures, the highest outflows of matter and a complexnucleosynthesis.Strong outflowing stellar winds of massive stars eject high amounts of matter that con-tribute to the chemical enrichment of the interstellar medium in a relatively short timescale.Moreover, it has been found that differences on a factor of two in the mass-loss rate affectsconsiderably the final fate of a massive star (Meynet et al., 1994; Smith, 2014). Therefore,a better understanding about massive stars and their evolution strongly requires accuratedetermination of their fundamental parameters, with the amount of matter released beingthe most relevant (Kudritzki & Puls, 2000; Puls et al., 2008). Subsequently, it is necessary tounderstand more in detail the mechanism responsible for driving the wind on massive starsin order to predict more accurately their mass-loss rates. The motivation for the presentthesis is then, to have a better understanding about the physics involved in the generation ofthe strong stellar winds on massive stars, in order to perform new prescriptions capable toquantify their mass-loss rates for future issues such as the already mentioned stellar evolutionand chemical enrichment. Stars are giant spheres of gas at high temperatures emitting energy (as electromagnetic radia-tion mainly, although, as we will see later, this is not the only way) to the interstellar medium.According with Prialnik (2009, Section 1.1), a star can be defined as a body satisfying thefollowing two conditions: • It is bound by self-gravity. From this condition it is undergone that stars must have aspherical shape because of gravity, or spheroidal in the case of the existence of axisym-134
CHAPTER 1. INTRODUCTION TO MASSIVE STARS metric forces such as rotation. • It radiates energy supplied by an internal source. This source is normally thermonu-clear energy, although sometimes gravitational potential energy may play a role due tocontractions and collapse.For the purposes of the present thesis, we will focus on the thermonuclear energy produced inthe interior of stars by thermonuclear fusion: reactions where atomic nuclei are transformedinto another species releasing the excess of mass as energy according with the Einstein’sequation E = mc . The most important thermonuclear reactions that happen in a star arethose where hydrogen is burnt in: the proton-proton process (pp) and the CNO cycle (carbon-nitrogen-oxygen). Energy produced by nuclear fusion, in form of photons (electromagneticradiation), passes through all the stellar structure from the nucleus until the surface to befinally released from there to the space.The stars are formed by the gravitational collapse of an interstellar gaseous cloud (beingthese nebulae, supernova remnants or molecular complexes). Because of the perturbationsproduced by shock waves from nearby supernova explosions or collisions with other clouds(Prialnik, 2009, Section 12.2) or simply because random matter movements, the gas cloud(originally thought as homogeneous) starts to form different regions with over-densities. Ifthe density is high enough, it will collapse on these regions in a process called fragmentation .Compression due to the disruption causes temperature to increase: it increases until condi-tions for nuclear reactions are reached. The radiation pressure , produced by the triggeredenergy, counteracts the gravitational collapse and each point of concentred matter reachesthe equilibrium again. So, it is generated a sphere in hydrostatic equilibrium whose photonsproduced in its inner parts because of nuclear fusion will be liberated into the interstellarmedium: a star is born.Since our initial cloud was not homogeneous on density, the different points where thecloud collapses do not concentrate the same quantity of matter. Some agglomerations willbe bigger and the others smaller. This leads to, the mass of the incipient stars be varied:some of them with a mass of one tenth of Sun mass only (hereafter solar-mass, or M (cid:12) ), andothers fifty or up to hundred times more massive than the Sun. It has been found that thedistribution of these masses is not homogeneous but it follows a distribution known as initialmass function ξ ( m ) ∝ m − α , being α a factor around ∼ .
35 for stars with masses greaterthan m = 1 M (cid:12) (Salpeter, 1955; Kroupa, 2001). In other words, whereas the more massive astar is, the less abundant in the Universe is.More mass for a star implies a greater compression in its core; a greater compressionin the core implies a higher temperature; and a higher temperature implies in turn more Proton-proton process consists on four hydrogen atoms fusing themselves to generate one atom of helium,together with energy as gamma radiation. This way of hydrogen-burning is predominant is low-mass stars.On the other hand,
CNO cycle consists in a cyclic chain of thermonuclear reactions that uses
CNO elements ,carbon, nitrogen and oxygen as catalysers. This is the predominant process in more massive stars. .2. MASSIVE STARS ∼ Kelvin) does not make the rate of nuclear reactions toincrease only, but also makes new kind of nuclear reactions appear whose existence would notbe possible in a not so extremely hot environment (such as ∼ × K). Moreover, thereis one more element playing a role: not only photons are emitted from the star, but alsoan outflow of particles called stellar wind . Because massive stars release more energy to theinterstellar medium (hereafter ISM), even stellar wind will be stronger in these stars. Hence,a star with a great mass will experiment an evolutionary path completely different from asolar like star.We will go deeper into the evolution followed by a massive star in the next sections, butwe will bring more theory about massive stars and stellar winds at first.
As we previously mentioned, massive stars are those with a stellar mass ten or more timesthe solar mass M (cid:12) .In accordance with their spectral classification they are stars type O and type B (typicallyreferred then as just OB stars ), which corresponds to the hottest spectral types, together withbeing the most luminous stars (Gray, & Corbally, 2009, Section 3.1). In their spectra lines ofionised helium, neutral helium and hydrogen are mainly observed, together with the so-called”metals”. Their properties are summarised in the Table 1.1.Initial Mass (cid:38) M (cid:12) Temperature (cid:38)
10 kKLuminosity (cid:38) L (cid:12) Mass-loss rate (cid:38) − M (cid:12) yr − Lifespan ∼ yrsTable 1.1: Properties of massive stars. CHAPTER 1. INTRODUCTION TO MASSIVE STARS
Figure 1.1:
Hertzsprung-Russell diagram. Location of massive stars appears delimited by the ma-genta ellipse.
Given their high brightness and temperature, massive stars are located in the top-leftregion of the Hertzsprung-Russell diagram (diagram that organises stars as function on theirluminosity and temperature, see Figure 1.1). Their peaks of emission (range of their spectrumwhere radiation reaches it maximum value) are in the ultraviolet region, the reason why wesee these stars as blue-coloured.As it has been previously mentioned, massive stars are also characterised by presenting astrong stellar wind which makes them lose a big amount of matter during their lifetimes (i.e.,a high mass-loss rate). This feature is crucial for their future evolutionary stages, reason whywe need to have a better understanding about what stellar wind consists to later go deeperinto evolution of stars with high mass. .2. MASSIVE STARS Previously, we explained why the massive stars have reserved a different future compared totheir smaller siblings. The main consequences of the big amount of mass that will determinelater the evolutionary track are, as we know: the greater mass-loss rate due to the strongerstellar wind and a more complex nucleosynthesis in the hotter core. Concerning this lastpoint, it will be relevant only at the moment of the final fate of the star, when it explodesas a supernova and eventually becomes a black hole. The evolutionary track is then, mostlyaffected by the high value of mass-loss rate.Hereunder, we will describe the evolutionary track for a 60 M (cid:12) star presented by Maeder &Meynet (1987) and available in the book Introduction to Stellar Winds (Lamers & Cassinelli,1999, Section 13.2).The most critical effect produced by the mass-loss rate upon a star is the ”dismantling”of this one, i.e., the stellar wind destroys the outer layers of the star letting the inner layersexposed. It also produces instability: a very massive star will never become a red supergiant(as low and intermediate-mass stars) because its great mass-loss rate prevents to reach equi-librium when core starts to burn helium and it must expand. Instead of that, there will be asa result an unstable star, variables and capable to sent shocks of matter towards the space:so-called Luminous Blue Variable stars (or simply LBV).Dismantling will produce later that, once the star have consumed all the hydrogen in itscore and begins to burn helium, the remnants of the hydrogen-burning processes will appear inthe surface of the star (helium due to the proton-proton process and nitrogen due to the CNO-cycle process mainly). As consequence, these elements (initially hidden in the inner layersof the star) turn to be observable in the stellar spectrum. Given that the star will exhibitan extended atmosphere , it will be seen in the spectrum broad emission lines of helium andnitrogen: we will observe a Wolf-Rayet star (WR star). These WRs are considered the finalstage in the life sequence of a massive star, previous to the final explosion as supernova andsubsequent stage as black hole.It is important to remark, however, that the previous description is a general screenshot,and it neglects many details that makes evolutive scenario more complex. Some of theseissues are: • The accurate constraint in the initial mass, to delimitate whether the massive star willreach the LBV and WR stages. Current acceptable values are given by (Crowther,2007) for solar-metallicity and non-rotating stars . M ∗ (cid:38) M (cid:12) : O V → WNh → LBV → WN → WC → SN Ic The atmosphere of a star is the boundary between the stellar interior and the ISM. Photosphere, thesurface of the star, is then the most inner layer of the atmosphere where photons can finally escape from theinterior, then, spectral lines are formed in these zones (Lanz, 2000). Nomenclature are: LBV for
Large Blue Variables , WN for
Nitrogen Wolf-Rayet stars , WC for
CarbonWolf-Rayet stars and SN for supernovae . CHAPTER 1. INTRODUCTION TO MASSIVE STARS M ∗ ∼ − M (cid:12) : O V ( → LBV) → WN → WC → SN Ic M ∗ ∼ − M (cid:12) : O V → (LBV)/RSG → WN ( → WC) → SN Ib M ∗ ∼ − M (cid:12) : O V → RSG → WN → SN II/Ib M ∗ ∼ − M (cid:12) : O V → RSG → BSG → SN II • Related with the first point: what is the impact of metallicity and rotational effectsupon the evolution of a star?Over the last decade, different studies have performed more detailed analysis about evo-lutionary tracks for a wide range of masses, rotational speeds and metallicities (Ekstr¨omet al., 2012; Georgy et al., 2012, 2013; Groh et al., 2019), all of them developed using theGeneva evolutive code (
Genec , see Chapter 6 for details). All these studies have been agreat contribution in order to understand the whole picture about evolution of massive stars.Yet mass-loss rates employed by them comes from a prescription for the stellar wind whichis not self-consistent. On the following chapter of this thesis we present a new prescriptionthat provides new theoretical values for the mass-loss rate, undergone from a self-consistentcalculation for the stellar wind. However, before giving the details, it is necessary to give abrief general picture of stellar winds.
We call stellar wind to the outflow of particles which, as same as the photons, are releasedfrom the photosphere of the star towards the interstellar medium.The main mechanism that explains the existence of stellar winds is the fact that in thephotosphere of the star the forces making hydrostatic equilibrium, total pressure from theinner part (generated by the radiation and by the gas of the star) and gravity are not longerin equilibrium at all. Pressure force coming from the interior wins over gravity and, due tothis imbalance, an outflow of matter is produced (Lamers & Cassinelli, 1999). This explainspartially the fact that in hotter stars, where the radiation pressure is higher, the stellar windis stronger.The two main parameters of the stellar wind, which can be determined by spectral analysis,are: • Terminal velocity: ( v ∞ ), understood as the asymptotic velocity reached by the particlesof the wind at large distances, measures in km s − . • Mass-loss rate: ( ˙ M ), corresponding to the amount of matter released by the star perunit time, measured in M (cid:12) yr − .Both terms allow us, for instance, to know the amount of energy and momentum releasedinto the interstellar medium. .3. STELLAR WIND ON MASSIVE STARS β -law,and it is expressed as: v ( r ) (cid:39) v ∞ (cid:16) − r r (cid:17) β , (1.1)with r the radial coordinate and: r = R ∗ (cid:34) − (cid:18) v v ∞ (cid:19) /β (cid:35) , being v the wind velocity in the photosphere of the star, i.e., v ( R ∗ ) = v . Here, β is a factorindicating how steep the increment in velocity along the path is: the higher the β value, theless pronounced the increase in speed will be (Fig. 1.2).Figure 1.2: Velocity field with different values of β (Lamers & Cassinelli, 1999). Nevertheless, it is important to remark that β value is just an empirical value whichdescribes an approximate behaviour of the wind and agrees with observed spectra with quiteacceptance. Kudritzki & Puls (2000) argued that the use of the parameter β to describe thevelocity field is only justified a posteriori once the fit is achieved. However, it is also possiblethe existence of velocity fields that can be determined (always with v ∞ as asymptotic limit)without a specific value of β . We will discuss this point with more details in the followingchapters, where we will perform our own self-consistent velocity fields.For the Sun-like stars, mass-loss rates are in the order of ˙ M ∼ − M (cid:12) yr − , which isten thousand times less intense than the minimum that massive stars exhibit (see Table 1.1).Hence, it is possible to see that the previously mentioned dismantling due to the stellar windis not significant for the ordinary stars. However, for massive stars it will play a key role thatwill condition the future evolutionary stages. In order to figure out how these evolutionary0 CHAPTER 1. INTRODUCTION TO MASSIVE STARS stages are determined by the features of the stellar wind, it is required to understand howthis wind is produced, how wind accelerates and what values can we derive for the stellarwind parameters.
In order to accurately characterise stellar winds and theoretically predict their parameters,it is necessary to analyse the physical processes behind. This is the motivation behind thisthesis study. For that reason, our first main objective is the obtention of self-consistent so-lutions (i.e., acceleration of the wind and hydrodynamics must be in agreement) for stellarwind parameters (mass-loss rate and terminal velocity) given different initial set of stellarparameters (effective temperature, mass, radius, metallicity and abundances mainly). Theself-consistent characteristic of the wind properties to be determined implies that they corre-spond to a unique solution for a given set of stellar parameters, and then they do not dependof a priori assumptions such as a β -law for velocity profile. Besides, we have to evaluate theinfluence on the final self-consistent solution of different approximations such as treatmentfor atomic populations and radiation field. And finally, we proceed to explore the potentialfuture works derived from the results for the self-consistent solutions.The process responsible for driving the wind for hot massive stars outwards is called line-driven , because it is produced by absorption and reemission of photons by the matter of thewind, and it will be extensively explained in Chapter 2. In Chapter 3, we will employ them-CAK theory for line-driven winds to calculate self-consistent solutions for the accelerationof the wind and their velocity and density profiles (i.e., the hydrodynamics of the wind) fora set of hot massive stars. Results, comparison with observed wind parameters and the newsynthetic model spectra obtained from the prescription are also presented. In Chapter 4,we calculate self-consistent hydrodynamically solutions for the stellar wind beyond m-CAKtheory under a full non-LTE scenario, the so-called Lambert-procedure . A complete analysison the differences with m-CAK prescription and their consequences are also included. InChapter 5, we perform several synthetic spectra for a set of massive stars from the self-consistent solutions calculated under the m-CAK prescription. In Chapter 6, we use themass-loss rates derived from the self-consistent hydrodynamics to perform new evolutionarytracks for standard non-rotating massive stars. Finally, summary and conclusions of our workis presented on Chapter 7. hapter 2
Line-driven Winds
Through this chapter, we will discuss in detail the mechanism that allows the wind on massivestars accelerate outwards. Most part of the content presented here can be also seen in theChapter 8 of the book
Introduction to Stellar Winds (Lamers & Cassinelli, 1999), togetherwith the study made by Puls et al. (2000) and the reviews from Kudritzki & Puls (2000) andPuls et al. (2008). 𝜽 𝜸 in 𝜸 out metal ion Figure 2.1:
Schema showing the absorption of a photon coming from the stellar photosphere, andtheir further reemission in a random direction.
The key to explain both the big amount of matter released to the space by a massivestar by means of its stellar winds and the high acceleration reached by this outflow lies212
CHAPTER 2. LINE-DRIVEN WINDS in the capacity of the atoms and ions along the wind of absorbing and reemiting photonscoming from the photosphere of the star. Because photons are coming from a specific regionwhereas the reemission is released in any arbitrary direction (see Figure 2.1), in average theions gain momentum generating then a force: we called to this process line-driving . Due toDoppler effect, line-driving process is not limited only to a specific rest frequency ν wherethe transitions take place but also occurs in a wide range of radial velocities where the relativefrequency matches with the rest frequency by means of: ν = ν (cid:18) v ( r ) c (cid:19) . (2.1)As consequence, the effect of the absorption and further reemission is enhanced along differentparts of the wind, resulting in a force larger than the gravitational one.Notice that the gain of momentum due to the previously described line-driving processesis applied over the individual ions instead the entire fluid. This leads to the fact that the ionshaving more lines (i.e., metal ions) accelerate more than hydrogen and helium. However, thehigher momentum gained by metal ions is shared with the more abundant and with lighterelements hydrogen and helium through Coulomb collisions. This scenario is fulfilled when thetimescale necessary to transfer momentum due to collisions is small enough to decelerate themetal ions before letting them escape. For stellar winds with high densities, this conditionis easily reached and then the acceleration from the line-driving process is transferred to allthe plasma, whereas for atmospheres with low mass-loss rates and large terminal velocitiesthe ionic runaway effect (i.e., ions that escape without sharing all their momenta) becomesmore relevant (Springmann & Pauldrach, 1992).The force generated by means of the line-driving mechanism will play a key role in thecalculation of the mass-loss rate: the line-acceleration g line . Because the line-driving processesinvolves a large number of features from the wind, it is necessary to take them into accountunder different levels of relaxation in order to properly calculate the line-acceleration. Butbefore analysing how to determine g line , it is important to examine how wind parameters areobtained from a specific value for line-acceleration. We call wind hydrodynamics to the coupled density ρ ( r ) and velocity v ( r ) profiles character-ising the wind of a star. Both fields are related each other by means of the isothermal andnon-rotating stationary equation of momentum on spherical coordinates : v dvdr = − ρ dPdr − GM ∗ (1 − Γ e ) r + g line , (2.2) Hereafter, every spatial equation is considered in spherical coordinates because of the geometry of thestar. Besides, all of them are solved assuming spherical symmetry, reason why only relevant spatial coordinateis the radius r . This assumption is done in order to have consistency with codes such as FASTWIND andCMFGEN, which solves their equations in 1D. .1. HYDRODYNAMICS OF THE WIND M = 4 πρ ( r ) r v ( r ) = constant , (2.3)with dP/dr being the pressure gradient, M ∗ the total mass of the star and Γ e the Eddingtonfactor. Γ e = σ e L ∗ πcGM ∗ . (2.4)Equation 2.2 shows that there will be a positive acceleration when left-hand side is greaterthan zero, i.e., the radiative components of the acceleration must be greater than the gaspressure and gravitational components above the photosphere: g rad = GM ∗ Γ e r + g line > ρ dPdr + GM ∗ r . (2.5)Besides, from Eq. 2.3 is clearly seen that both fields determine the wind parameters: mass-loss rate ˙ M and terminal velocity v ∞ . Therefore, calculation of wind hydrodynamics meanscalculation of these wind parameters which are later constrained by observations in the stellarspectra.In order to solve equation of momentum, we can assume isothermal conditions in order tointroduce the equation of state for an ideal gas: P = a ρ ( r ) , (2.6)with a being the isothermal sound speed: a := (cid:115) k B T eff µm H (2.7)and with k B being the Boltzmann’s constant, µ the mean particle mass and m H the hydrogenatom mass. In this case, equation of momentum becomes the equation of motion : (cid:18) − a v (cid:19) v dvdr = 2 a r − GM ∗ (1 − Γ e ) r + g line . (2.8)As consequence of the fact that we are using a constant temperature equivalent to the effectivetemperature (isothermal wind), sound speed is considered as a constant. Besides, notice thefact that Eq. 2.8 does not longer depend explicitly on density ρ ( r ). Actually, dependence It is important to remark that, sometimes in the literature the equation of momentum is written using theterm g rad instead g line . In that case, the radiative acceleration due to the continuum (i.e., not produced bythe line-driving but by the photons doing Thomson scattering with the electrons of the wind) is not includedin the gravitational term and then the equation of momentum reads as: v dvdr = − ρ dPdr − GM ∗ r + g rad . Unfortunately there is not a consensus about the strict name of the equations of motion and momentum.On this Thesis we are using the names assigned by Puls et al. (2008), but Lamers & Cassinelli (1999) callsEq. 2.2 as equation of motion and Eq. 2.8 as equation of momentum. CHAPTER 2. LINE-DRIVEN WINDS on density is implicitly included inside term for g line , one of the reasons why it is necessaryprovide a hydrodynamics to solve line-acceleration self-consistently. However, as we willexamine later in Chapter 4, for the case where g line is directly obtained from the solution ofthe radiative transfer equation and its dependence on density is not directly known, we canconsider ρ (and therefore mass-loss rate ˙ M ) as a free parameter.Treatments employed to solve Eq. 2.8 depends on the formulation used to calculate theline-acceleration and what variables were considered for its calculation; therefore, we will focusthe discussion through this chapter into knowing how g line is determined under the ’classical’line-driving theory performed in the decade of the 70s: the CAK (and later m-CAK) theory. Lucy & Solomon (1970) described the mechanism that drives the strong stellar winds observedin hot stars: the so-called radiation driven winds. The process of absorption and furtherre-emission of photons and Coulomb interactions previously described at the beginning ofthis chapter is the mechanism responsible to give momentum to the wind of hot stars, thenproducing an outwards line-force. According to these authors, the effectivity of line-drivingmechanism lies in the fact that the most part of the atomic transitions involved come fromthe ultraviolet resonance lines, which in turn is where the peak of radiation field for hotstars is located. The foundation of the theory of radiation driven winds was later developedby Castor et al. (1975, hereafter CAK theory), who, based on the Sobolev and the point-star approximations, modelled the line-acceleration analytically in terms of the accelerationproduced by electron scattering times a force multiplier factor. This factor represents thecontribution of absorption and re-emission processes depending on the optical depth only,and it was parametrised by two constant parameters through the wind, namely k and α . In order to understand the theoretical bases of line-driven winds, let us analyse first the case ofgaining momentum from a single line. Assuming the wind is optically thick for this transition,implies that all the photons coming from the photosphere which could be absorbed by theatomic transition will do. Therefore, due to Doppler effect there is not a unique frequency ν to be absorbed by the wind, but a range going from the rest frequency (in the region of thewind where v = 0) to the external parts of the wind where v = v ∞ (with ν given by Eq. 2.1).As consequence, the total photospheric radiation being absorbed by the wind (per unit oftime) is given by: L ∗ , abs = 4 πR ∗ (cid:90) ν (1+ v ∞ /c ) ν F ν dν , (2.9) .2. THE M-CAK THEORY F ν being the flux at the line frequency: F ν = L ν πr , satistying: L ∗ = (cid:90) ∞ L ν dν . (2.10)The associated momentum is then given by: p line ,i = L ∗ , abs c = 4 πR ∗ c (cid:90) ν (1+ v ∞ /c ) ν F ν dν . (2.11)The total momentum gained by the wind will be then equivalent to the sum of all thelines where the wind is fully optically thick: p line = 4 πR ∗ c N (cid:88) i =1 (cid:90) ν (1+ v ∞ /c ) ν F ν dν . (2.12)Hence, line-acceleration could be obtained once an accurate calculation of all the possibleoptically thick lines that take part in the line-driving process. However, full optically thicklines are an idealisation for illustrative purposes. Instead, each line will present a differentvalue for its opacity depending on the atomic properties of its associated transition and thezone of the wind where the absorption takes place. Then, it is necessary to define the massabsorption coefficient κ l for a single line: κ l ( r ) = πe m e c f l N l ρ ( r ) (cid:18) − N u N l g l g u (cid:19) , (2.13)with N l and N u being the number density of the ion for the lower and upper excitation levelsrespectively, g l and g u being the respective statistical weights and f l being the oscillatorstrength of the atomic transition . Considering that the photon generated by this atomictransition has an energy of hν , absorption coefficient can be written in terms of frequency: κ ν = κ l φ (∆ ν ) , (2.14)with φ (∆ ν ) being a normalised profile function describing the range on the frequency domainwhere the transition occurs: φ (∆ ν ) = 1 √ π ∆ ν G e − (∆ ν/ ∆ ν G ) , (2.15)with ∆ ν G being the Gaussian width of the profile, determined by the thermal (see Eq. 2.42)and turbulent motions of the wind.∆ ν G = ν c (cid:114)
23 ( (cid:104) v th (cid:105) + (cid:104) v turb (cid:105) ) . (2.16)Because of the dependence on density and the influence of Doppler effect, absorptioncoefficient depends then not only on the frequency of the photon and the atomic information The usage of the subindex l seems to be misleading, meaning ”lower” for the atomic densities and thestatistical weights, and meaning ”line” for the oscillator strengths. However, since physically a spectral line isproduced due to an atomic transition from a specific lower level of excitation to an upper one (or viceversa),it is possible to use the subindex l without leading into errors. CHAPTER 2. LINE-DRIVEN WINDS
Figure 2.2:
Schema showing the velocity field v z = v ( r ) cos θ along the line of sight. A photonemitted by the photosphere at frequency ν p will be absorbed by the wind at a specific frequency ν determined by both the velocity field v ( r ) and the angle θ . Image taken from Lamers & Cassinelli(1999). for the involved ion, but also depends implicitly on the point of the wind there the absorptionwill take place. For the case of considering the wind for a star with finite disk (i.e., the starhas a specific radius R ∗ and it is not assumed as a point source), we can parametrise thelocation on the wind there the transition occurs in terms of the line of sight z : z = r cos θ = (cid:112) r − p , (2.17)with θ being the angle between the radial direction and the line of sight, and p being the impact parameter , perpendicular to the line of sight (see Fig. 2.2). Notice that it is alwayssatisfied that p ≤ R ∗ .Location of this z point represents the place where the peak of the profile function φ (∆ ν )is located. Since the neighbourhood around z represents the region where the radiationcoming from the photosphere is absorbed, the atmosphere beyond this zone becomes opaquefor photons at that range of frequencies (although, as we explained before, this opacity is notinfinite). In order to represent mathematically this situation, we introduce the optical depth defined as: τ ν p ( z ) = (cid:90) ∞ z κ ν ( z ) ρ ( z ) dz , (2.18)and using Eq. 2.14 and Eq. 2.13: τ ν p ( z ) = πe m e c f l (cid:90) ∞ z N l ( z ) (cid:18) − N u ( z ) N l ( z ) g l g u (cid:19) φ (∆ ν ) dz . (2.19)The interval covered by ∆ ν depends on the wind velocity at the point z , which is givenby Doppler effect: ∆ ν ( z ) = ν p (cid:18) − v z ( r ) c (cid:19) − ν = ν p (cid:18) − zr v ( r ) c (cid:19) − ν (2.20) .2. THE M-CAK THEORY z , which implies moreover to know the behaviour of our atomic populations infunction of the radius. However, as we have pointed out previously, the absorption takes placeinside a region determined by the velocity of the wind. If we consider a steep velocity gradient dv/dr , the ∆ v where the transition occurs leads into a even narrower ∆ r ; and moreover, anarrow Gaussian profile (see Eq. 2.16), will also lead into a narrow region on radius. Giventhese scenarios, it is possible to take this assumptions in order to simplify the calculation ofthe optical depth. The consideration of the absorption region as small enough in length (and therefore z ) inorder to reduce the zone to a single point, is called Sobolev approximation (Sobolev, 1960).In this limit, the profile φ (∆ ν ) becomes a delta-function and so: τ ν p ( z ) = πe m e c f l (cid:90) ∆ ν ( z →∞ )∆ ν ( z ) N l ( z ) (cid:18) − N u ( z ) N l ( z ) g l g u (cid:19) φ (∆ ν ) (cid:18) dzd ∆ ν (cid:19) d ∆ ν , = πe m e c f l N l ( r s ) (cid:18) − N u ( r s ) N l ( r s ) g l g u (cid:19) (cid:18) dzd ∆ ν (cid:19) r s , with z ≤ r s , = κ l ( r s ) ρ ( r s ) (cid:18) dzd ∆ ν (cid:19) r s . (2.21)Thanks to Sobolev approximation, optical depth can be expressed as a step function wherethe absorption coefficient (see Eq. 2.13) is evaluated in a single point (called Sobolev point ).In order to obtain the final expression for τ s (where s also means Sobolev), we solve thederivative ( dz/d ∆ ν ) r s using Eq. 2.20 and Eq. 2.17: (cid:18) dzd ∆ ν (cid:19) r s = cν p (cid:20) sin θ v ( r ) r + cos θ dvdr (cid:21) − r s . (2.22)For simplicity, we can define µ = cos θ = z/r . Then, Sobolev optical depth in terms ofthe rest frequency where the absorption takes place is written as: τ s,ν = κ l ( r s ) ρ ( r s ) cν p (cid:20) (1 − µ ) v ( r ) r + µ dvdr (cid:21) − r s , = κ l ( r s ) ρ ( r s ) cν (cid:20) v ( r ) r + µ (cid:18) dvdr − vr (cid:19)(cid:21) − r s , = κ l ( r s ) ρ ( r s ) cν (cid:16) vr (cid:17) − (cid:20) µ (cid:18) rv dvdr − (cid:19)(cid:21) − r s , = κ l ( r s ) ρ ( r s ) λ (cid:16) vr (cid:17) − (cid:2) µ σ (cid:3) − r s , (2.23)where we have introduced the variable σ : σ = d ln vd ln r − rv dvdr − . (2.24)8 CHAPTER 2. LINE-DRIVEN WINDS
Besides, we have used the approximation ν p (cid:39) ν valid for not relativistic wind velocities.Decomposing the absorption coefficient the Sobolev optical depth gives: τ s,ν = πe m e c λ f l N l ( r s ) (cid:18) − N u ( r s ) N l ( r s ) g l g u (cid:19) (cid:18) r/v µ σ (cid:19) r s . (2.25)Hence, optical depth finally depends on the atomic and wind conditions on the point whereit is being evaluated. This point comes from the fact that, under Sobolev approximation,we are considering the interaction region where the atomic transition takes place as beinginfinitely narrow. However, this assumption is not real at all because of random motion(thermal and turbulence velocities), expressed in the term ∆ ν G inside Eq. 2.15. In spiteof that, Sobolev approximation works well as far as the velocity field v ( r ) is larger enoughcompared with v th and v turb (both usually in the order of ∼ −
30 km s − ) and whether thevelocity dv/dr is large enough to keep the absorption region as narrow as possible (in order tokeep the particles density N l and N u almost constant). This means, the region should havea width of: ∆ r (cid:39) v th ( dv/dr ) ≡ L s , (2.26)with L s being the Sobolev length . Previous conditions are easily satisfied downstram fromsonic point, with sound speed in the order of ∼
25 km s − and a high acceleration on thewind, outwards. Thus, we establish the range of validity for Sobolev approximation (andtherefore m-CAK line-acceleration) to be from the sonic point to infinite. This discussion willbe retaken in Chapter 3. Once we have obtained an analytical expression for optical depth, we proceed to derive anexpression for the acceleration due to line-driving.Previously, we had derived a temptative expression for the momentum gained by thewind given the ideal case of being absorbing all the radiation at that range of frequencies.In reality, the amount of momentum to be gained is proportional to the opacity of the windfor that line at that point of the wind, which implies a complex problem of radiative transfersince we need to know how much flux of energy from the photosphere reaches to the point r .However, thanks to the Sobolev approximation this problem is easily reduced to consider thatthe photon emitted by the photosphere will interact with the wind at that transition only inthe Sobolev point r s . Thus, the amount of radiation accelerating the wind at point r willdepend only on the intensity coming from the stellar photosphere and the local conditionsaround r s and then determined by the optical depth τ ν ( r s ).If we assume that radiation coming from the photosphere is homogeneous, the intensityof radiation reaching the Sobolev point is given by: I ν p ( µ ) = I ∗ ν p e − τ νp ( µ ) , (2.27) .2. THE M-CAK THEORY τ ν p we adapted the expression given by Eq. 2.18 to run from the photosphere tothe radius r where intensity is being evaluated : τ ν p ( µ ) = (cid:90) r phot κ ν p ρd l , = (cid:90) r phot κ l ρφ (∆ ν ) d l , = κ l ( r ) ρ ( r ) (cid:18) d l d ∆ ν (cid:19) (cid:90) ∆ ν ( r )∆ ν (phot) φ (∆ ν ) d (∆ ν ) , = τ ν ( r ) (cid:90) ∆ ν ( r )∆ ν (phot) φ (∆ ν ) d (∆ ν ) . (2.28)Notice that absorption coefficient and density have been extracted from the integral usingthe Sobolev approximation. The integral can be defined as:Φ(∆ ν µ ) = (cid:90) ∆ ν ( r )∆ ν (phot) φ (∆ ν ) d ∆ ν , (2.29)with ∆ ν µ being ν p − ν (1 + µv ( r ) /c ), analogous to Eq. 2.20. Besides, once again we assumethat the intensity emitted by the photosphere is almost constant in the interval covered by∆ ν µ , so we can set the photospheric intensity I ∗ ν p (cid:39) I ∗ ν where the last term is the continuumintensity at the rest frequency.However, we need to consider the radiation coming from all the angles of the stellarphotosphere, and for that reason we focus in the mean intensity J ν p , interpreted as theintensity per steradian of radiation at the frequency ν p . This is equivalent at the half of theintegration of the intensity I ν p over µ = cos θ , from the minimal possible value µ ∗ on thephotosphere of the star ( µ ∗ = cos θ ∗ = (cid:112) − ( R ∗ /r ) ) to the maximal value at θ = 0. J ν p = 12 (cid:90) µ ∗ I ν p dµ = I ∗ ν (cid:90) µ ∗ e − τ ν Φ(∆ ν µ ) dµ . (2.30)Since φ (∆ ν ) is a normalised function, we can later calculate the mean intensity integratedover all frequencies as:¯ J = I ∗ ν (cid:90) µ ∗ (cid:90) ∆ ν µ =+ ∞ ∆ ν µ = −∞ φ (∆ ν µ ) e − τ ν Φ(∆ ν µ ) d (∆ ν ) dµ , = I ∗ ν (cid:90) µ ∗ (cid:90) e − τ ν Φ(∆ ν µ ) d Φ(∆ ν ) dµ , = I ∗ ν (cid:90) µ ∗ − e − τ ν τ ν dµ , ¯ J ( r ) = I ∗ ν (cid:90) µ ∗ − e − τ ν τ ν dµ . (2.31)Eq. 2.31 gives us the amount of energy per second (power) per unit of surface and persteradian gained at radius r by means of the absorption at rest frequency ν . Since this We use the variable d l to determine an infinitesimal segment instead dl , in order to not confuse with l for’line’. CHAPTER 2. LINE-DRIVEN WINDS energy is absorbed from the radiation coming from the stellar photosphere, the momentumgained by the wind at radius r is equal to absorption coefficient κ ν times the flux (= 4 π ¯ J )divided by the speed of light c : g line = 2 πc (cid:90) µ ∗ (cid:90) ∞−∞ κ ν p (∆ ν ) I ν p ( µ ) d (∆ ν ) µdµ , = 2 πc I ∗ ν (cid:90) µ ∗ (cid:90) ∞−∞ κ ν p (∆ ν ) e − τ ν Φ(∆ ν ) d (∆ ν ) µdµ , = 2 πc κ l I ∗ ν (cid:90) µ ∗ − e − τ ν τ ν µdµ . (2.32)The integral can be solved if we consider the limit case that the star is a single point,and therefore the only valid possible value for µ ∗ is 1. This simplification is called the pointsource limit , and it was introduced by Castor et al. (1975), and it constitutes a fundamentalpart of the CAK theory. Hence, in the point source limit with µ ∗ = 1 we have: g line = 2 πc κ l I ∗ ν − e − τ s τ s , = 2 πc κ l I ∗ ν − e − τ s τ s . (2.33)We have eliminated the subindex 0, because hereafter all frequencies are only referred torest frequency (or rest wavelength). For the same reason, we substitute the subindex ν for theoptical depth and we replace it for s , to emphasise that we are using Sobolev approximation.If µ = 1, Eq. 2.23 becomes: τ s ( r ) = κ l ρλ (cid:16) vr (cid:17) − (cid:20) rv dvdr (cid:21) − r s , = κ l ( r ) ρ ( r ) λ (cid:18) dvdr (cid:19) − . (2.34)Therefore, for line-acceleration: g line = 2 πc κ l I ∗ ν (1 − e − τ s ) dv/drκ l ρλ , = 2 πc νI ∗ ν (cid:18) dvdr (cid:19) − e − τ s ρ . (2.35)Finally, considering the relation between the intensity and the luminosity: I ∗ ν = L ν πr , we obtain the following analytical expression: g line ( r ) = L ν ν πr c (cid:18) dvdr (cid:19) − e − τ s ρ ( r ) . (2.36)This expression shows that line-acceleration, besides the classical dependence on the fre-quency ν and the luminosity L ν coming from the photosphere, has a important dependence .2. THE M-CAK THEORY r . This last dependence is perhaps the most important, because it showsthat line-acceleration will present a different behaviour depending on the strength of τ s . Forexample, for the case of lines with small τ s (called optically thin lines ), we can approximate e − τ s (cid:39) − τ s and then: g thin line = L ν ν πr c (cid:18) dvdr (cid:19) τ s ρ ( r ) , = L ν ν πr c (cid:18) dvdr (cid:19) κ l λdv/dr , = L ν πr c κ l , (2.37)= κ l c F ν . (2.38)This result shows that for optically thin lines, acceleration is given mainly by the pho-tospheric flux and is intrinsically dependent on density by means of absorption coefficient,but it is independent on velocity gradient. On the other hand, if the line has τ s (cid:29) optically thick lines ): g thick line = L ν ν πr c ρ ( r ) (cid:18) dvdr (cid:19) , (2.39)where line-acceleration is independent on absorption coefficient, but also on photospheric andhydrodynamic conditions (Puls et al., 2000, 2008). Line ensemble
However, we are interested in evaluate the resulting acceleration produced by all the lines involved in the line-driving process, being them optically thick or thin. In order to obtainthat expression, the work of Castor et al. (1975) consisted in the search of an expression easyto sum and analyse. Combining Eq. 2.13 and Eq. 2.34, the full expression for Sobolev opticaldepth is: τ s ( r ) = πe m e c λf l N l (cid:18) − N u N l g l g u (cid:19) (cid:18) dvdr (cid:19) − . (2.40)This value will vary from line to line, because each spectral line carries its own informa-tion about atomic populations and statistical weights. But the hydrodynamical components(density and velocity gradient) will be the same for all lines because they depend on r only.Hence, it is convenient to separate both components, in order to define a new optical depthindependent on atomic information. For that purpose, Castor et al. (1975) have introducedthe new variable t , the independent optical depth or optical depth for an expanding atmo-sphere (Abbott, 1982), defined as: t ≡ σ e v th ρ ( r ) (cid:18) dvdr (cid:19) − , (2.41) Hereafter and during all the CAK procedure, we refer t simply as optical depth only, omitting the word independent . To avoid confusions, we will explicitly specify when we refer to the classical optical depth τ . CHAPTER 2. LINE-DRIVEN WINDS with σ e = 0 .
325 cm g − (Castor et al., 1975; Abbott, 1982) being the electron scatteringopacity and v th the mean thermal velocity of the protons: v th = (cid:114) k B T eff m H . (2.42)The inclusion of thermal velocity is important, because random thermal movements playsa role enhancing the range of frequencies to be absorbed by means of Doppler effect. We candefine the Doppler enhancement due to thermal motions as:∆ ν D = ν v th c = v th λ . (2.43)No less important, the component of the full optical depth τ depending on atomic infor-mation only is read as: η line = πe m e c λf l N l (cid:18) − N u N l g l g u (cid:19) σ e v th ρ , = πe m e c ρσ e ∆ ν D f l N l (cid:18) − N u N l g l g u (cid:19) . (2.44)This term η line represents the ratio between line to electron scattering opacity, and it isfulfilled that τ line = t η line . With these new definitions, line-acceleration can be rewrittenfrom Eq. 2.36 to: g line ( r ) = L ν ν πr c (cid:18) dvdr (cid:19) − e − η l t ρ ( r ) , = L ν ν πr c σ e v th − e − ηt t , = σ e F ν c ∆ ν D − e − ηt t . Here, we have used the relation between luminosity and flux from Eq. 2.10 and we haveomitted the subindex l for η . This expression for line-acceleration has the advantage of beingwritten in a very similar way to the standard acceleration due to radiation pressure (i.e., thatacceleration produced by electron scattering interactions): g e = σ e F c = σ e L ∗ πr c . (2.45)Therefore, line-acceleration can be expressed as the radiative acceleration due to electronscattering, multiplied by some factor representing the contribution of all the lines involved inthe line-driving process. This factor was defined by Castor et al. (1975) as the force multiplierfactor as: M ( t ) = g line g e = (cid:88) lines ∆ ν D F ν F − e − ηt t . (2.46)Notice the fact that the force multiplier is now defined not as function of radius, but asoptical depth t . The great advantage of this procedure is, M ( t ) can be easily parametrisedby a simple power law. M ( t ) = k t − α , (2.47) .2. THE M-CAK THEORY k and α the line-force [multiplier] parameters .This was the most important result from the revolutionary work done by Castor, Abbottand Klein in 1975, and for that reason is called CAK theory . The authors calculated anacceptable force multiplier following Eq. 2.46, by means of the sum of the spectral lines ofthe ions of carbon. The posteriori challenge was, the inclusion of more precise atomic data(oscillator strengths, statistical weights, excitation energies) for all the individual ions involvedin the line-driving process, together with the calculation of an accurate thermodynamicaltreatment in order to calculate accurate atomic populations N l and N u .This pioneer study opened the door to the possibility of obtaining a solution for equationof motion (Eq. 2.8), and thereafter every study dedicated to stellar wind on massive stars istotally or partially based on CAK theory. However, later studies included relaxations to someof the main assumptions of the CAK theory, together with other considerations not takeninto account by Castor et al. (1975). These new improvements leaded to the generation ofthe modified CAK (m-CAK) theory . Details about these changes are given in the followingsection. Seven years after the introduction of CAK theory, Abbott (1982) performed a detailed cal-culation of the line-force multiplier M ( t ) taking into account the contribution of a full setof atomic line transition data for elements from hydrogen to zinc. Moreover, the line-forcemultiplier was calculated over a fixed grid of optical depths t and also for different valuesof diluted electron densities N e /W , for a wide range of stellar temperatures. The most re-markable result from this study, was the inclusion of an extra exponential dependence on thediluted electron density N e /W for M ( t ): M ( t ) = k t − α (cid:18) N e, W ( r ) (cid:19) δ , (2.48)with N e, being the electron density in units of 10 cm − , δ our third line-force parameterand W ( r ) the dilution factor, i.e, the function showing how radiation is ’diluted’ through thewind: W ( r ) = 12 (cid:32) − (cid:114) − R ∗ r (cid:33) . (2.49)This expression can be easily obtained in the limit τ ν (cid:28) M ∝ N δe, can be obtained inspecting Fig. 2.3, where is clearly seenthat force-multiplier increases with electron density in an almost exponential fit, and its ex-planation lies in the relationship between the number of allowed transitions and the ionisationstage. Higher electron densities leads into lower stages of ionisation, as it is shown by Sahaequation and its version for ionisation equilibrium in expanded atmospheres given by Mihalas4 CHAPTER 2. LINE-DRIVEN WINDS
Figure 2.3:
Figures a from Abbott (1982), showing the dependence of the force multiplier M ( t ) onelectron density. (1978, Eq. 5-46): N i +1 N i = 2 (cid:18) πm e k B h (cid:19) / T R √ T e N e /W U i +1 U i e − EikBTe , (2.50)being m e the electron mass, T R and T e the radiative and electron temperatures respectively, U i the partition functions of each ionisation stage and E i the ionisation energy. Lower ionisedstages have more lines, which allows the absorption of more radiation by line-driving.Due to the point-star approximation ( µ ∗ = 1 for Eq. 2.32) the derived hydrodynamicalvalues for mass-loss rates given by Castor et al. (1975) and Abbott (1982) were overestimated.The explanation for this disagreement lies in the fact that all incoming photons of Fig. 2.1in reality enter not always with θ = 0 (as in the point source approximation) but a range ofvalues for the angle which will reduce the effective value for cos θ . Pauldrach et al. (1986) andFriend & Abbott (1986) relaxed this point source approximation and considered the finitedisk shape of the star. This modification consisted in the inclusion of a finite disk correctionfactor over the point-source multiplier factor. D f = M disk ( t ) M point ( t ) , = (1 + σ ) α +1 − (1 − σµ ∗ ) α +1 (1 − µ ∗ )( α + 1) σ (1 + σ ) α +1 , (2.51)with α being the line force parameter. With this modified theory (hereafter m-CAK) theysolved the equation of momentum and obtained improved theoretical results, in better agree-ment with the observed mass-loss rates. .2. THE M-CAK THEORY β –law velocity profile.Indeed, from the m-CAK simulations performed by Pauldrach et al. (1986) it was determineda velocity profile following the β − law (Puls et al., 2008): v ( r ) = v ∞ (cid:18) − R ∗ r (cid:19) β , β = 0 . β = 0 . β = 2 . β had to be relaxed in order to fit the spectra, so these theo-retical values are not longer valid. This simplified description of the velocity field is widelyused as input in radiative transfer and NLTE model-atmosphere codes such as FASTWIND(Santolaya-Rey et al., 1997; Puls et al., 2005) or CMFGEN (Hillier, 1990b; Hillier & Miller,1998; Hillier & Lanz, 2001) to calculate NLTE synthetic spectra. In this procedure, stellarand wind parameters (terminal velocity and mass-loss rates) are treated as free and are de-termined by varying them to adjust synthetic profiles to observed ones. In the particularcase of CMFGEN, the final solution considers a full NLTE treatment and it also provides aradiative acceleration calculated beyond the simplifications undergone from m-CAK (such asSobolev approximation), but this radiative acceleration is not consistent with the β –law setas input for hydrodynamics. Kudritzki & Puls (2000) argued that the usage of β − law for thevelocity field is only justified a posteriori once the fit is achieved. However, there are otherapproaches that coupled the hydrodynamics with comoving frame radiative transfer, see e.g.Sander et al. (2017) or Krtiˇcka & Kub´at (2010, 2017), that do not use a β − law velocityprofile.Nevertheless, in spite of the disadvantages associated to m-CAK theory (such as Sobolevapproximation and corrections over point source assumptions) and efforts trying to obtaina self-consistent solution beyond it, this remains being a valid reference for the calculationof line-acceleration. Moreover, m-CAK theory provides us the great chance to perform afast self-consistent solution beyond assumptions from β –law, and for that reason we havechosen this regime to execute our analysis for stellar winds in Chapter 3. In spite of that, afull parallel analysis of self-consistent solutions in a complete NLTE regime beyond m-CAKprescription will be done in Chapter 4.6 CHAPTER 2. LINE-DRIVEN WINDS hapter 3
Self-consistent Solutions in theframe of m-CAK Theory
In this chapter we go in details about the solutions obtained for line-driven winds in the frameof the m-CAK theory described on Section 2.2. Solution of equation of motion is obtainedusing the line-acceleration g line with the derived values of the line-force multiplier parameters k , α and δ . This procedure has the enormous advantage of saving a lot of computational effort,and therefore it allows the execution of a large set of models in a short time. However, theprice to be paid is the adoption of several assumptions whose consequences are discussed. Dueto scarce works involving NLTE (non-local thermodynamic equilibrium) calculations of theline-force parameters (Pauldrach et al., 1986; Puls et al., 2000; Kudritzki, 2002; Pauldrach,2003; Noebauer & Sim, 2015), it was difficult to obtain from the m-CAK hydrodynamicsthe velocity profiles and mass-loss rates, thus, the community of massive star researchersstarted to use the β -law velocity profile instead of the proper hydrodynamics. This simplifieddescription of the velocity field is widely used as input in radiative transfer and NLTE model-atmosphere codes such as FASTWIND (Santolaya-Rey et al., 1997; Puls et al., 2005) orCMFGEN (Hillier, 1990a; Hillier & Miller, 1998; Hillier & Lanz, 2001) in order to calculatesynthetic spectra. In those procedures, stellar and wind parameters (terminal velocity andmass-loss rates) are treated as free parameters and are determined by adjusting them to fitsynthetic line profiles with observed ones. However, there are other approaches that coupledthe hydrodynamics with comoving frame radiative transfer, see e.g. Sander et al. (2017) orKrtiˇcka & Kub´at (2010, 2017), that do not use a β -law velocity profile.Calculations of line-force wind parameters coupled with hydrodynamics are necessary toderive self-consistent values of the velocity profiles and the mass-loss rates. Moreover, theseline-force parameters depend non-linearly on the stellar parameters, chemical abundances,and atomic data via the wind driven mechanism. Therefore, to obtain the line-force parame-ters it is necessary to calculate the total acceleration produced by the contribution of hundreds378 CHAPTER 3. SOLUTIONS IN THE FRAME OF M-CAK THEORY of thousands lines involved in the absorption and re-emission processes (i.e., line-acceleration, g line ) which requires reliable atomic data, as they are essential to perform line-statistics cal-culations.The number of contributing lines to the line-driven acceleration depends on the windopacity and it is strongly coupled to the wind density and velocity profiles. To solve thishighly non-linear problem an iterative procedure is required to satisfy both: line-statisticsand m-CAK hydrodynamics.In this chapter, we calculate self-consistent solutions to obtain accurate m-CAK line-forceparameters ( k, α, δ ) and wind properties of hot massive stars. The hydrodynamics is providedby the code HydWind (Cur´e, 2004), whereas abundances have been adopted from Asplundet al. (2009). Final self-consistent line-force values must correspond to an unique solutionobtained when line-force parameters, velocity profile and mass-loss rate simultaneously con-verge. Hence, we present here a new set of m-CAK self-consistent line-force parameters for T eff ≥
32 kK and log g ≥ . v inf /v esc on log g , which was not previously known.We have to mention that this chapter corresponds to the manuscript Gormaz-Matamalaet al. (2019), also referred as Paper I. ( k, α, δ ) Using the expression for line-acceleration from the force-multiplier M ( t ) (Eq. 2.48), we re-write (Eq. 2.2) as: v dvdr = − ρ dPdr − GM ∗ (1 − Γ e ) r + g e kt − α (cid:18) N e, W ( r ) (cid:19) δ , (3.1)being k , α and δ the already mentioned line-force parameters (see Eq. 2.48). Assuming anisothermal ideal gas P = a ρ , equation of momentum is transformed into equation of motion(Eq. 2.8): (cid:18) − a v (cid:19) v ( r ) dvdr = 2 a r − GM ∗ (1 − Γ e ) r + g e k ( σ e v th ρ ) − α (cid:18) N e, W ( r ) (cid:19) δ (cid:18) dvdr (cid:19) α . (3.2)Rewriting dv/dr as v (cid:48) , we can define momentum equation as a function F ( r, v, v (cid:48) ) = 0.Besides, we use Eq. 2.45 to modify g e and equation of continuity Eq. 2.3 to express density .1. EQUATION OF MOTION WITH LINE-FORCE PARAMETERS ( K, α, δ ) 39in terms of the mass-loss rate. Thus, F ( r, v, v (cid:48) ) reads: F ( r, v, v (cid:48) ) = (cid:18) − a v (cid:19) r vv (cid:48) + GM ∗ (1 − Γ e ) − a r − σ e L ∗ πc k (cid:32) σ e v th ˙ M πr v (cid:33) − α (cid:18) N e, W ( r ) (cid:19) δ (cid:18) dvdr (cid:19) α , = (cid:18) − a v (cid:19) r vv (cid:48) + GM ∗ (1 − Γ e ) − a r − σ e L ∗ πc k (cid:18) πσ e v th ˙ M (cid:19) α (cid:18) N e, W (cid:19) δ (cid:0) r vv (cid:48) (cid:1) α , = (cid:18) − a v (cid:19) r vv (cid:48) + GM ∗ (1 − Γ e ) − a r − h ( r ) C (cid:0) r vv (cid:48) (cid:1) α = 0 , (3.3)with the constant C being: C = σ e L ∗ πc k (cid:18) πσ e v th ˙ M (cid:19) α , and the function h ( r ) being: h ( r ) = (cid:18) N e, W (cid:19) δ = (cid:34) N e × (cid:32) − (cid:114) − R ∗ r (cid:33)(cid:35) δ . Despite the fact that the ionisation density h ( r ) = N e, /W is not strictly constant, h ( r ) C was assumed as constant under the classical m-CAK formulation because the exponent δ hasvalues typically in the order of (cid:46) .
15 (Abbott, 1982; Pauldrach et al., 1986). This is theso-called fast solution . However, for cases when δ takes higher values (in the order of (cid:38) . h ( r ) is not longer constant and a new hydrodynamical solutions may arise: the so-called δ -slow solution (Cur´e, 2004; Cur´e & Rial, 2004).Due to the exponent α , Eq. 3.3 is a non-linear differential equation in which we will lookfor a monotonically increasing function v ( r ). This solution must be unique, and starting fromsubsonic velocities ( v < a ) close to the stellar photosphere reaching a supersonic asymptoticterminal value v ∞ (cid:29) a at large radius. Castor et al. (1975) demonstrated that there aredifferent cases for v ( r ) satisfying Eq. 3.3, but no one of them fulfil the conditions previouslymentioned. However, it is possible to couple one of the solutions starting from subsonicregion to one starting from the infinite, which match on a specific point called critical or singular point (see Section 8.7.2 of Lamers & Cassinelli, 1999, for details). This criticalpoint r c is determined once the eigenvalue (mass-loss rate) of the equation of motion isnumerically calculated, being then the respective eigenfunction, the formal solution for thehydrodynamics of the wind (Friend & Abbott, 1986; Pauldrach et al., 1986; Kudritzki et al.,1989). For this calculation, it is necessary to provide a location for this critical point as anstandard method to solve the equation. On this work however, hydrodynamics is calculatedusing the code HydWind (Cur´e, 2004), where the equation of motion is calculated by meansof finite-difference method, modified to handle singular points (Nobili & Turolla, 1988). Thiscode has the advantage of obtaining, according to the initial trial solution, different solutionswith other critical points in addition to the standard one. Here, we use the letter h instead g in order to avoid confusions with acceleration. CHAPTER 3. SOLUTIONS IN THE FRAME OF M-CAK THEORY
Although in Section 2.2 we introduced k , α and δ merely as constant parameters to fit theforce-multiplier, it is possible to find physical meanings for the resulting line-acceleration andhydrodynamics depending on their final values. Then, physical interpretation of the line-forceparameters (see, e.g., Puls et al., 2000) are: • The k parameter, which takes values between 0 and 1, is directly proportional to theeffective number of driving lines, and is related to the fraction of the photospheric fluxwhich would have been blocked by all lines if they were optically thick and overlappingeffects were not considered. Higher values of k are obtained at higher densities and,therefore, higher mass-loss rates. In addition to the dependency on ρ ( r ), k presentsalso a strong dependence with metallicity and temperature due to the large number ofdriving lines: a lower temperature implies lower ionization stages, and thus more lines,therefore a higher M ( t ). More lines (above a given threshold line-strength) are alsopresent for higher metallicities.The overlapping of two or more spectral lines produces an overestimation in the calcu-lated value of k . On the other hand, k is underestimated when multi-scattering effectsare not taken into account (i.e., the summation in M ( t ) considers only direct photo-spheric radiation, and not radiation reprocessed in the wind). As was pointed out byPuls (1987), the inclusion of both effects might cancel, at least for O stars, and the effective k becomes moderately reduced. In this work, we have not considered theseeffects, therefore, our k values should be maximum. • The α parameter, which usually takes values between 0.45 and 0.75, is related to theexponent of the line-strength distribution function, and quantifies also the ratio of theline acceleration from optically thick lines to the total one (for details, see Puls et al.,2008). Higher values of α implies both high mass-loss rates and terminal velocities inthe resulting hydrodynamics. • The δ parameter represents the change in the ionisation throughout the wind. Accordingwith classical literature (Lamers & Cassinelli, 1999), it takes lower values, rarely higherthan 0.1. However, it has been found that, high values of δ ( (cid:38) .
25) makes the wind’slow’, yielding a different wind solution (Cur´e et al., 2011) because C from Eq. 3.3cannot be longer treated as constant. Besides, according with Puls et al. (2000) δ takesan ”exact” value of 1/3 for neutral hydrogen as a trace element.These interpretations are coming from previous studies; hence, our work will consist inthe analysis about whether these statements agree with our results or not, and of what newconsiderations we can include for discussion. Besides, some studies have pointed out thatthe line-force parameters are a function of radius (Schaerer & Schmutz, 1994) or can beconsidered in a piecewise constant structure (Kudritzki, 2002). Nevertheless, in this workwe will consider k , α and δ as constants throughout all the wind because their variation is .2. CALCULATION OF THE M ( T ) FACTOR M ( t ) factor When we talk about the calculation of the line-force parameters k , α and δ , we actually meanthe calculation of the line-force multiplier M ( t ) using Eq. 2.46. To do that, we perform ascript called Alfakdelta27 , including the following improvements:i) a larger atomic line list,ii) a quasi-NLTE approach for the ionisation equilibrium,iii) a NLTE radiative stellar flux andiv) an optical depth range consistent with the wind structure.Then we test it for one single-step (i.e., without iterations) first and after we execute thewhole iteration procedure until the convergence of line-force parameters, velocity profile andmass-loss rate is achieved.
To calculate the line-acceleration and obtain a proper value of M ( t ), Abbott (1982) estab-lished that it is necessary to sum the contribution of hundreds of thousands of spectral linesparticipating in the line-acceleration processes. Indeed, that work was pioneer in the inclusionof a larger atomic database taking into account more ions than those previously consideredby Castor et al. (1975). Therefore, aiming to get the most accurate value of M ( t ), we decidedto employ around ∼
900 000 line transitions, whose atomic data were obtained (and modifiedin format) from the atomic database list used by the code CMFGEN (Hillier, 1990a; Hillier& Miller, 1998). We have chosen this database because it is the most complete one available,specially containing an extensively complete number of atomic transitions for heavy elementslike iron and nickel (which contributes with ∼
80% of all the spectral lines). Specifically, wehave extracted information related to energy levels, degeneracy levels, partition functions andoscillator strengths f l , which are necessary to calculate the absorption coefficient η line (seeEq. 2.44) of each line in terms of lower ( l ) and upper ( u ) level populations n l and n u , andtheir statistical weights g l and g u .Before continuing, it is important to remark the differences on notation and definitionscompared with other studies. The most important atomic information about transitions forthe calculation of M ( t ) is the oscillator strength f l , defined as the dimensionless quantitythat express the probability of absorption of reemission between the both energy levels cor-responding to the line in question. However, because a change in the notation some authors Atomic data used here are that one which were updated by DJH in 2016 ( http://kookaburra.phyast.pitt.edu/hillier/cmfgen_files/atomic_data_15nov16.tar.gz ). CHAPTER 3. SOLUTIONS IN THE FRAME OF M-CAK THEORY
Elem. Ions N o lines Elem. Ions N o linesH I 599 He I − II 1 342Li I − III 273 Be I − IV 76B I − V 85 C I − IV 25 421N I − VI 8 691 O I − VI 6 851F I − VI 187 Ne I − VI 30 880Na I − VI 8 138 Mg I − VI 7 136Al I − VI 5 613 Si I − VI 2 839P I − VI 3 331 S I − VI 15 455Cl I − VI 534 Ar I − VI 27 376K I − VI 287 Ca I − VI 37 556Sc I − VI 322 Ti I − VI 791V I − VI 920 Cr I − VI 779Mn I − VI 688 Fe I − VI 278 923Co I − VI 489 Ni I − VI 492 341
Table 3.1:
Atomic elements and ionisation stages used to calculate M ( t ). Range of frequency of thespectral lines goes from UV to IR range. express the factor η as: η line = πe m e c ρσ e ∆ ν D f l g l (cid:18) N l g l − N u g u (cid:19) , = πe m e c ρσ e ∆ ν D ( gf ) (cid:18) N l g l − N u g u (cid:19) . (3.4)Where the former oscillator strength is now multiplied by the lower statistical weight, creatingthen the so called gf -factor (see, for example, notation used in Puls et al., 2000). Moreover,some other authors use the Einstein A coefficient instead of the oscillator strength, which arerelated by the formula: A ul = 8 π ν e (cid:15) m e c g l g u f [s − ] . (3.5)Even when Einstein A coefficient do not appear explicitly on m-CAK notation, some authorssuch as Noebauer & Sim (2015) have derived their oscillator strength from them (privatecommunication in 2018).Elements and ionisation stages considered in this work are listed in Table 3.1. FollowingAbbott (1982), we consider ions up to ionisation stage VI only. The total number of lines perelement is also specified. Line-acceleration is calculated over the contribution of numerous transitions for every elementat every ionisation stage present in the wind. Abbott (1982) determined the ionization degrees .2. CALCULATION OF THE M ( T ) FACTOR (cid:18) N i +1 N i (cid:19) LTE = 2 (cid:18) πm e k B h (cid:19) / T R √ T e N e /W U i +1 U i e − EikBTe , (3.6) being T R , T e the radiation and electron temperatures, respectively, and E i the ionisationenergy from stage i to i + 1. More precise treatment called approximate NLTE (hereafterquasi-NLTE) has been used by, e.g., Mazzali & Lucy (1993) and Noebauer & Sim (2015). Herethe ionisation balance is determined by the application of the modified nebular approximation(Abbott & Lucy, 1985). Following this treatment, the ratio of number densities for twoconsecutive ions can be expressed in term of its LTE value, multiplied by correction effectsdue to dilution of radiation field and recombinations: N i +1 N i ≈ (cid:18) N e W (cid:19) − [ ζ i + W (1 − ζ i )] (cid:114) T e T R (cid:18) N i +1 N e N i (cid:19) LTE , (3.7) where ζ i represents the fraction of recombination processes that go directly to the groundstage. Eq. (3.7) is an alternative description to the one given by Puls et al. (2005), who in-cluded a different radiative temperature dependence in the wind, which is specially importantin the far UV region of the spectrum that is not optically thick.Modifications in the treatment of atomic populations X i , being i the excitation level, arealso based on the work of Abbott & Lucy (1985). It is necessary to make distinction betweenmetastable levels (with no permitted electromagnetic dipole transitions to lower energy levels)and all the other ones: (cid:18) X i X (cid:19) = (cid:16) X i X (cid:17) LTE metastable levels, W ( r ) (cid:16) X i X (cid:17) LTE others.Atomic partition functions, U i (necessary for Saha’s equation and the calculation of atomicpopulations), are calculated following the formulation of Cardona et al. (2010), i.e.: U i = U i, + G jk e − ε jk /T + m n − e − ˆ E n ∗ jk /T , (3.8)where U i, are the constant partition functions, ˆ E n ∗ jk is the mean excitation energy of thelast level of the ion, n is the maximum excitation stage to be considered, while G jk , ε jk and m are parameters tabulated by Cardona et al. (2010).The advantage of this treatment is that it provides values for atomic partition functionsexplicitly as function of temperature and implicitly of electron density, giving a more accurateionization balance. Following Noebauer & Sim (2015), the temperature will be treated as aconstant ( T R = T e = T eff ). Then, for a specific value of ( T eff , N e ), the ratio between numberdensities of ionization stage j and i (for a specific Z -element) is calculated by a matrix(hereafter ionization matrix) given by: D Z,i,j = N j N i = (cid:89) i ≤ k Final value of N e /W ( r ) as function of stellar radius even when N e, is set as constantinput (black solid line), after one iteration (single dashed line), after two iterations (dashed-dottedline) and after five iterations (red solid line). In reference to the abundances of the different chemical elements, these were adopted fromthe solar abundances given by Asplund et al. (2009). However, these can be easily modifiedto evaluate stars with non-solar metallicity (see Section 3.3).At this point, it is necessary to remark that previous authors (Abbott, 1982; Noebauer& Sim, 2015) have considered the diluted-electron density N e /W as constant throughout thewind. Nevertheless, to calculate δ , M ( t ) must be evaluated considering changes in the ionisa-tion stages, and therefore, N e ( r ) /W ( r ). Since, the calculation of electron density depends onthe ionisation stages of each specie which in turn are functions of N e , we deal with a couplednon-linear problem. To obtain a solution, we use the following formula to calculate (as initialvalue) the electron number density: N e, = ρ ( r ) m H X H + 2 X He X H + 4 X He , (3.10)being m H the hydrogen atom mass, and X H and X He the abundances of hydrogen and helium,respectively.We used this initial electron density to start calculating the ionisation matrix and tore-calculate both atomic populations and electron density iteratively: N e ( r ) = (cid:18) X H D , , D , , + X He ( D , , + 2 D , , )1 + D , , + D , , (cid:19) × ρ ( r ) X H + 4 X He . (3.11)Convergence of N e is easily obtained after just a few iterations (see Fig. 3.1). It is impor-tant to remark, that alternatively we can start iterations using N e, /W as a constant valuefollowing Abbott (1982) and Noebauer & Sim (2015) instead of starting using Eq. 3.10, andanyway the final converged value for N e ( r ) is the same. .2. CALCULATION OF THE M ( T ) FACTOR Together with an accurate treatment of atomic populations and electron density, Eq. 2.46requires as an input the radiation field in the terms of F ν /F .Simplest expression for the radiation field comes from the black-body Planck’s law: B ν ( T ) = 2 hν c e − hν/k B T − . (3.12) × × × × × × Ang W / ( m s r ) Figure 3.2: Comparison of radiation fields F ν obtained by black-body Planck’s radiation law (cyan,dashed) and those ones obtained from the models of Kurucz (red) and Tlusty (dark blue), imple-mented for a star with T eff = 40 000 K and log g = 4 . But this back-body scenario assumes a full LTE, which is not valid for stellar atmospheresbecause of the transport of energy and matter; it is required then to solve the equation ofradiative transfer in order to incorporate effects due to opacity and hence obtain an accurateradiation field. Since this calculation is beyond the scope of the present work, we proceed tosimply employ the already calculated flux fields for different stellar models. Some of the mostcommon radiation field models used by the stellar wind community are those performed byKurucz (1979, in LTE only) and the more modern Tlusty models (Hubeny & Lanz, 1995;Lanz & Hubeny, 2003, in both LTE and NLTE). Differences among these models are presentedin Fig. 3.2, where black-body radiation Planck’s law is shown together with a Kurucz and a Tlusty model.The usage of already performed stellar models for the flux field was also implementedbefore. For example, Abbott (1982) used the radiation fields from Kurucz’ models, whereasNoebauer & Sim (2015) from a black-body. In this work, we use the radiation field calculatedby the NLTE line-blanketed plane-parallel code Tlusty .6 CHAPTER 3. SOLUTIONS IN THE FRAME OF M-CAK THEORY The overlap effects among ten of thousands of spectral lines are not considered when wesum the contributions to the force-multiplier M ( t ) across the wind. However, line blanketingeffects are partially considered as we are using Tlusty radiation field in the calculations of M ( t ). α k δ - - - - - - l og ( M ) v ∞ Figure 3.3: Left: Values of α , k and δ as a function of the iteration number, starting from differentinitial values. Different initial values (iteration 0, not shown) converge to the same final self-consistentsolution. Right: idem as left, but for the mass-loss rate and terminal velocity. Previous studies by Abbott (1982) and Noebauer & Sim (2015) have considered a fixed rangefor the optical depth t to fit the force multiplier (Eq. 2.47). However, given the definition of t (Eq. 2.41), it is clear that the optical depth range is constrained by the physical propertiesof the stellar wind (density and velocity profiles). For this reason, calculations presented in .2. CALCULATION OF THE M ( T ) FACTOR t .Because m-CAK theory is based upon Sobolev approximation (Sobolev, 1960, see alsoSection 2.2.2 of the present Thesis) in this work we will use as upper and lower limits for theoptical depth t , its values at the sonic point and at infinity (usually r ∼ R ∗ ), respectively.It is important to remark that although t decreases outwards it never reaches zero . Therefore,it is always possible to define a proper t range. Velocity profile and mass-loss rate from hydrodynamics are required to calculate the line-acceleration g line . At the same time, line-force parameters fitted from g line , are necessary tosolve the m-CAK hydrodynamic equations and obtain the mass-loss rate and velocity profile.Therefore a self-consistent iterative procedure must be implemented to solve this couplednon-linear problem.Our procedure is the following:i. Using a β − law profile with a given mass-loss rate, initial values for the line force pa-rameters ( k , α , δ ) are calculated.ii. A numerical solution of the equation of motion (Eq. 3.2) is obtained with HydWind ,getting an improved hydrodynamics: v ( r ) and ˙ M .iii. A new force multiplier is calculated.iv. New line-force parameters ( k i , α i , δ i ) are fitted from M ( t )v. Steps ii - iv are iterated until convergence.Convergence is usually obtained after ∼ − i , i + 1) get a value for (cid:107) ∆ p (cid:107) = (cid:107) p i +1 − p i (cid:107) ≤ − , where p is a line-force parameter andthis condition should be satisfied for each one of these parameters.Right panel of Fig. 3.3 shows the convergence of the mass-loss rate (top panel) and theterminal velocity (lower panel) as function of the procedure’s iterations. Both values dependnon-linearly on the stellar and line-force parameters. See also Section 8.4 of Lamers & Cassinelli (1999, Introduction to Stellar Winds ) This condition lies in the fact that, at larger distances, both density ρ and velocity gradient dv/dr decreaseas ∼ r − , cancelling each other. This code solves the m-CAK equation of motion with an eigenvalue that depends on the mass-loss rate.At the location of the singular point, both solution branches (singular point to stellar surface and singularpoint to infinity) are smoothly merged to obtain the velocity profile, see Pauldrach et al. (1986); Friend &Abbott (1986); Cur´e (2004) for details. CHAPTER 3. SOLUTIONS IN THE FRAME OF M-CAK THEORY previous studies present work T eff N e /W δ k α k α [kK] [cm − ]A 30 1 . × . × . × . × . × . × . × . × . × . × Table 3.2: Comparison of k and α parameters from Abbott (A) and Noebauer & Sim (N), with ourone single-step results. To compare our line-force parameters with the results obtained by Abbott (1982) and Noe-bauer & Sim (2015), we use just one single-step. Following these authors, δ and N e /W areset as input and the optical depth range is fixed between − < log t < − 1. The selection ofa fixed interval of log t does not require any velocity field structure. Furthermore we haveconsidered Kurucz’ and black-body fluxes to reproduce Abbott (1982) and Noebauer & Sim(2015) calculations, respectively. Then, starting from a β –law and a ˙ M , we calculate k and α (single-step). These results are shown in Table 3.2. The coefficients of determination, R -Squared, for α and k (respectively) between previous and our single-iteration results are:i) R α = 0 . 87 and R k = 0 . 93 for T eff ≥ 40 000 K;ii) R α = 0 . R k = 0 . 81 for T eff ≥ 30 000 K.We conclude that our calculations reproduced previous results, now using a modern atomicdatabase and abundances. This section is focused on the resulting values obtained for both line-force and wind parametersfollowing the methodology previously described. The following results are computed self-consistently with the methodology detailed in Sec-tion 3.2. .3. M-CAK RESULTS Teff = = ℳ ( t ) Teff = = ℳ ( t ) Teff = = ℳ ( t ) Teff = = - - ℳ ( t ) Figure 3.4: Force-multiplier M ( t ) as function of t for some stellar models presented on Table 3.3with T eff = 45 000 K and log g = 4 . T eff = 40 000 K and log g = 3 . T eff = 36 000K and log g = 3 . T eff = 32 000 K and log g = 3 . t where the fits for ( k, α, δ ) have been adjusted. CHAPTER 3. SOLUTIONS IN THE FRAME OF M-CAK THEORY log g = = = = k α 32 34 36 38 40 42 440.020.040.060.080.100.120.140.16 T eff δ Figure 3.5: Behaviour of line-force parameters ( k, α, δ ) as a function of the effective temperature(in kK), for different surface gravities and metallicities. Circles represent models with log g = 4 . g = 3 . 8, stars: log g = 3 . 6, and triangles: log g = 3 . 4. Black dashed lines are for modelswith solar metallicity and grey dashed lines for Z = Z (cid:12) / .3. M-CAK RESULTS g from 3 . . t -range, and the fitted m-CAK line-force. In addition, wecalculated the corresponding wind solution using HydWind , and their error margins werederived considering variations of ∆ T eff = ± g = ± . 05, and ∆ R ∗ = ± . R (cid:12) in thestellar radius, keeping constant the line-force parameters.Convergence has been checked for each solution. Figure 3.4 shows the final resulting M ( t )given by the last iteration for different four models from Table 3.3 at their respective rangesof t . Due to the quasi-linear behaviour of the logarithm of the force-multiplier, parameters k and α are easily fitted and their values can be considered constant throughout the wind (seeSect. 3.3.2). To fit δ in the M ( t )– N e /W plane, it is necessary to perform an extra calculationof M ( t ) using a slightly different value for the diluted-electron density. Last column of thistable shows the ratio between our mass-loss rate and the one calculated using Vink’s recipe(Vink et al., 2001), with v ∞ /v esc = 2 . M SC / ˙ M Vink = 0 . ± . 2. As we have not included in our proceduremulti-line nor line-overlapping processes, we support Puls (1987) conclusion that these effectsare somewhat canceled, because we do not observe relevant discrepancies in the mass-lossrates when a comparison with Vink’s recipe is performed.In Fig. 3.5, we observe clear trends for the behaviour of the ( k, α, δ ) parameters with T eff , log g , and Z . While k increases, δ decreases as function of the effective temperature,for both metallicities. It is interesting to remark the influence of the surface gravity on theresulting line-force parameters, values for k and δ decrease as the gravity decreases. Noticethat globally our line-force parameter results are similar to the values obtained in previousworks (Puls et al., 2000; Kudritzki, 2002; Noebauer & Sim, 2015). However, we found animportant dependence on log g as a result of the hydrodynamic coupling in the self-consistentprocedure.On the other hand, the behaviour of α depends on the metallicity, it increases with effectivetemperature for solar abundance, but for low abundance and low gravities, it slowly decreaseswith temperature. Moreover, the change in α is more significant for log g than for T eff : adifference in ∆ log g ± . α ∼ . 04, whereas variations on ∆ T eff = ± α ∼ . HydWind (hereafter Abbott’s procedure). We found that ˙ M increases with effectivetemperature and metallicity and decreases with gravity. This behaviour is similar to theone obtained using Abbott’s procedure, but the self-consistent calculated mass-loss rates are2 CHAPTER 3. SOLUTIONS IN THE FRAME OF M-CAK THEORY about 30% larger.From the mass-loss results tabulated in Table 3.3, a simple linear relationship for solar-likemetallicities (with a coefficient of determination or R –squared, R = 0 . M Z =1 . =10 . × log (cid:18) T eff (cid:19) − . × log g + 0 . × ( R ∗ /R (cid:12) ) − . , (3.13)and for metallicity Z/Z (cid:12) = 0 . R = 0 . M Z =0 . =11 . × log (cid:18) T eff (cid:19) − . × log g + 0 . × ( R ∗ /R (cid:12) ) − . , (3.14)where ˙ M are given in 10 − M (cid:12) yr − .These relationships could be considered analogous to that given by Vink et al. (2000) toobtain theoretical mass-loss rates for solar-like metallicities. However, the advantage of ourdescription is that it depends only on stellar parameters and we do not need to consider thevalue of v ∞ /v esc . It is important to remark however, that this formula has been derived forthe following ranges: • T eff = 32 − 45 kK • log g = 3 . − . • M ∗ /M (cid:12) ≥ . v ∞ is almost constant with respect to the effective temperature, but itdecreases as a function of log g and Z . On the other hand, Abbott’s procedure results do notshow the same behavior and exhibit a maximum in the T eff interval. It is important to remember that the range of optical depths used to calculate our self-consistent line-force parameters is defined along almost all the atmosphere of the star, i.e.,downstream from the sonic point. This procedure improves the criterion used by Abbott .3. M-CAK RESULTS - - - - - - l og ( M ) 32 34 36 38 40 42 44 - - - - - - eff l og ( M ) V ∞ 32 34 36 38 40 42 4415002000250030003500 T eff V ∞ Figure 3.6: Left: Behaviour of mass-loss rate as a function of effective temperature (in kK) fordifferent abundances and gravities. Top panel is for self-consistent calculations and bottom panelis for Abbott’s procedure, now including the finite-disk correction factor. Symbol description is thesame as Fig. 3.5. Right: same as left panel, but for terminal velocities. (1982), who arbitrarily defined the parameters at t = 10 − . This value sometimes laysoutside the optical depth range here defined, as it was shown in Fig. 3.4.To analyse the change on the line-force parameters due to the selection of the t -range,we define four different intervals inside the whole range of t , and compute these parametersin each range. Table 3.4 summarises these calculations. Regarding the uncertainties of ourprocedure in the terminal velocities, these are of the same order as the uncertainties due to theerrors in the determination of the stellar parameters in the range 32 000 K < T eff < 40 000K, while, the uncertainties in ˙ M are much lower than the ones produced by variations ofstellar parameters. These small uncertainties indicate that it is a good approximation toconsider line-force parameters as constants throughout the wind. Due to the fact that theentire t -range represents the physical conditions of almost all the wind, we recommend to usethe complete optical depth range to derive the line-force parameters.For T eff < 30 000 K, we found that log M ( t ) is not longer linear with respect to log t and4 CHAPTER 3. SOLUTIONS IN THE FRAME OF M-CAK THEORY the corresponding line-force parameters can be approximated to a linear piecewise description.Due to this reason, we establish that our set of self-consistent solutions describes stellarwinds for effective temperatures and log g in the range 32 000 − 45 000 K and 3 . − . In order to know whether our calculations reproduce realistic physical features observed inhot stars, we calculated synthetic spectra for three O-type stars using FASTWIND. Weselected some stars in the range of the considered T eff , trying to cover the extreme casesof temperature and log g . We chose first the O4 I(n)fp star ζ -Puppis (HD 66811) because ithas been extensively studied (Puls et al., 1996; Repolust et al., 2004; Puls et al., 2006; Sotaet al., 2011; Bouret et al., 2012; Noebauer & Sim, 2015). Because HydWind allows the optionto include a rotational velocity, self-consistent solutions for ζ -Puppis consider a v sin i = 210km s − which is a value in agreement with previous authors. Mentioned authors have alsoadopted independently different set of stellar and wind parameters, which are summarisedin Table 3.5. Here, the wind parameters were determined by Repolust et al. (2004). Pulset al. (2006) has used Repolust’s parameters and derived clumped mass-loss rates from H α ,IR and radio, using analytical expressions for the corresponding opacities, whereas Bouret etal. (2012) used CMFGEN. Both calculations include clumping, so these results correspondto a clumped mass-loss rate. On the other hand, the mass-loss rate given by Noebauer &Sim (2015) was obtained using their Monte-Carlo radiation hydrodynamics (MCRH) methodassuming a homogeneous media ( f cl = 1 . ζ -Puppis:i) the ”conventional” ( d = 460 pc) andii) the one given by Sahu & Blaauw (1993, d = 730 pc).We examine here the ”conventional” case with R ∗ /R (cid:12) = 18 . 6. We can observe fromTable 3.5 (last row), that our new calculated mass-loss rate agree quite well with the valuefrom Puls et al. (2006).Figure 3.7 shows the observed spectra (kindly provided by D. J. Hillier) and the result-ing synthetic spectra for ζ -Puppis. Stellar parameters are taken from Puls et al. (2006,see Table 3.5) and wind parameters from our self-consistent procedure ( ˙ M SC = 4 . × − FASTWIND uses the clumping factor f cl ≥ f cl = 1 representing the smooth limit), where f cl = 1 /f if the inter-clump medium was void (Sundqvist, & Puls, 2018). On the other hand, CMFGEN-clumping isrepresented by the so-called volume filling factor f , which scales homogeneous and clumped mass-loss ratesunder the relationship ˙ M hom = ˙ M clump / √ f (notice that this f takes values between 0 and 1). .4. SYNTHETIC SPECTRA λ ( Ang ) N o r m a li se dF l u x H α λ ( Ang ) N o r m a li se dF l u x H β λ ( Ang ) N o r m a li se dF l u x H δ λ ( Ang ) N o r m a li se dF l u x He II 4200 λ ( Ang ) N o r m a li se dF l u x He II 4541 λ ( Ang ) N o r m a li se dF l u x He II 4686 λ ( Ang ) N o r m a li se dF l u x He I 4026 λ ( Ang ) N o r m a li se dF l u x He I 4471 λ ( Ang ) N o r m a li se dF l u x He I 4922 Figure 3.7: Resulting FASTWIND spectra for ζ -Puppis with T eff = 39 kK, log g = 3 . R ∗ /R (cid:12) =18 . M = 4 . × − M (cid:12) yr − . Clumping factors are f cl = 1 . f cl = 5 . f cl = 9 . M (cid:12) yr − ). We calculated three synthetic spectra with different clumping factors: f cl = 1 . f cl = 5 . f cl = 9 . 0. The best fit is for f cl = 5 . 0, which is the sameclumping factor found by Puls et al. (2006) with their ˙ M = 4 . × − M (cid:12) yr − . Moreover,we also include the synthetic spectra obtained with the self-consistent solution (see Fig. 3.8),calculated using the stellar parameters given by Bouret et al. (2012, see Table 3.5) and Na-jarro et al. (2011). The best fit is achieved when we use a clumping factor of f cl = 5 . 0. Theseresults suggest that the real stellar parameters lies in the neighbourhood given by Puls et al.(2006) and Najarro et al. (2011).The observed spectrum for HD 163758 (O9 I) has been obtained from UVES-POP database .We calculated the synthetic spectra for this star (see Fig. 3.9) using stellar parameters fromBouret et al. (2012) and wind self-consistent parameters (see Table 3.6) with different clump-ing factors, the best fit is for f cl = 6 . CHAPTER 3. SOLUTIONS IN THE FRAME OF M-CAK THEORY λ ( Ang ) N o r m a li se dF l u x H α λ ( Ang ) N o r m a li se dF l u x H β λ ( Ang ) N o r m a li se dF l u x H δ λ ( Ang ) N o r m a li se dF l u x He II 4200 λ ( Ang ) N o r m a li se dF l u x He II 4541 λ ( Ang ) N o r m a li se dF l u x He II 4686 λ ( Ang ) N o r m a li se dF l u x He I 4026 λ ( Ang ) N o r m a li se dF l u x He I 4471 λ ( Ang ) N o r m a li se dF l u x He I 4922 Figure 3.8: Resulting FASTWIND spectra for ζ -Puppis with T eff = 40 kK, log g = 3 . R ∗ /R (cid:12) =18 . M = 5 . × − M (cid:12) yr − . Clumping factors are f cl = 1 . f cl = 5 . f cl = 9 . Last spectrum corresponds to the O3.5 V star HD 164794, also obtained from UVES-POPdatabase. Stellar parameters were extracted from Krtiˇcka et al. (2015), as shown in Table 3.6.Contrary to previous cases, the best fit is obtained for the homogeneous model ( f cl = 1 . 0, seeFig. 3.10).In view of these first results, our self-consistent iterative procedure takes us quickly intothe neighborhood of the solution that reproduces the observed wind spectra for O-type stars. We have developed a self-consistent methodology to calculate the line-force parameters andderived consequently mass-loss rates and velocity profiles. We found that mass-loss rateis about 30% larger than the one obtained using Abbott’s procedure (non self-consistentcalculation). .5. DISCUSSION λ ( Ang ) N o r m a li se dF l u x H α λ ( Ang ) N o r m a li se dF l u x H β λ ( Ang ) N o r m a li se dF l u x H δ λ ( Ang ) N o r m a li se dF l u x He II 4200 λ ( Ang ) N o r m a li se dF l u x He II 4541 λ ( Ang ) N o r m a li se dF l u x He II 4686 λ ( Ang ) N o r m a li se dF l u x He I 4026 λ ( Ang ) N o r m a li se dF l u x He I 4471 λ ( Ang ) N o r m a li se dF l u x He I 4922 Figure 3.9: Resulting FASTWIND spectra for HD 163758 with T eff = 34 . g = 3 . R ∗ /R (cid:12) = 21 . M = 3 . × − M (cid:12) yr − . Clumping factors are f cl = 5 . f cl = 6 . f cl = 7 . It is well known the scaling relation for the terminal velocity in the frame of the line-drivenwind theory. This relation (Puls et al., 2008) reads: v ∞ ≈ . (cid:114) α − α v esc . (3.15)This is an approximation of the formula found by Kudritzki et al. (1989, their Eqs. 62 to 65).In Fig. 3.11 we have plotted v ∞ /v esc versus (cid:112) α/ (1 − α ) using the results from Table 3.3.Contrary to the expected result (Eq. 3.15) for solar abundances, we find a different linearbehaviour, which strongly depends on the value of log g . This is a new result that comesfrom applying our self-consistent procedure. The m-CAK equation of momentum shows aninterplay between the gravity (log g ) and the line-force term. This balance of forces definesthe location of the singular point and therefore fixes the value of ˙ M . As a consequence,the velocity profile depends also on the value of log g . This result cannot be obtained from8 CHAPTER 3. SOLUTIONS IN THE FRAME OF M-CAK THEORY λ ( Ang ) N o r m a li se dF l u x H α λ ( Ang ) N o r m a li se dF l u x H β λ ( Ang ) N o r m a li se dF l u x H δ λ ( Ang ) N o r m a li se dF l u x He II 4200 λ ( Ang ) N o r m a li se dF l u x He II 4541 λ ( Ang ) N o r m a li se dF l u x He II 4686 λ ( Ang ) N o r m a li se dF l u x He I 4026 λ ( Ang ) N o r m a li se dF l u x He I 4471 λ ( Ang ) N o r m a li se dF l u x He I 4922 Figure 3.10: Resulting FASTWIND spectra for HD 164794 with T eff = 43 . g = 3 . R ∗ /R (cid:12) = 13 . M = 2 . × − M (cid:12) yr − .Clumping factors are f cl = 5 . f cl = 2 . f cl = 1 . Eq. 3.15 which is an oversimplification of this non-linear coupling. However, Eq. 3.15 presentsa fair fit when Z = Z (cid:12) / 5, where the dependence of the slope on log g is weak since the radiationforce is driven by few ions.The dependence of v ∞ /v esc on log g yield that stars with solar abundances present anintrinsic variations of v ∞ /v esc in the range 2 . − . 7, as shown in Fig. 3.11. This range mightexplain the scatter observed in the hot side of the bi-stability jump shown by Markova & Puls(2008, in their Fig. 12). In this section we want to compare our theoretical results with the ones obtained from line-profile fittings for homogeneous (unclumped) winds with a β –law, and the mass-loss (recipe)from Vink et al. (2000).Table 3.7 shows our results for two O-type star reported by Bouret et al. (2005): HD 96715, .5. DISCUSSION ( α /( - α )) V ∞ / V e sc Figure 3.11: v ∞ /v esc versus (cid:112) α/ (1 − α ). For each set of log g values there is a linear dependencefor Z (cid:12) . Slope 2.25 of Eq. 3.15 is also displayed. For sub-solar abundance there is a unique linearrelationship (see text for details). Symbol description is the same as in Fig. 3.5 T eff = 43 . g = 4 . 0, and HD 1904290A, T eff = 39 kK, log g = 3 . 6. These results wereobtained for the self-consistent solution together with the ones after just one iteration startingfrom a β –law. It is observed that models starting from a β –law largely overestimate theterminal velocity and slightly underestimate the mass-loss rate. Self-consistent calculationsfind a fairly good agreement to both: the observed mass-loss rate and terminal velocity. Forthe mass-loss rate in this Figure, we have included the result calculated using Vink et al.(2000) recipe. It is clear that our self-consistent method gives values of ˙ M much closer to theobserved ones.We also apply our self-consistent procedure to objects analysed by means of FASTWINDadopting unclumped winds. For that purpose, we also examine some field Galactic O-typestars from Markova et al. (2018). Table 3.8 summarises our results, where we found a fairagreement between observed and calculated mass-loss rates (see Fig. 3.12). These resultsconfirm that our methodology delivers the proper mass-loss rate for the ranges in T eff andlog g given above. Below these thresholds, mass-loss rates present larger values comparedwith both: observational and Vink’s theoretical values. This is probably due to the fact thatthe line-force multiplier is not longer a linear function of t (in the log-log plane, see Fig. 3.4),and the line-force parameters are not constant throughout the wind.However, it is important to remark that uncertainties of ∆ T eff ∼ ± g ∼± . CHAPTER 3. SOLUTIONS IN THE FRAME OF M-CAK THEORY - - - - - l og ( M ) 30 32 34 36 38 40 42 44100015002000250030003500 T eff ( kK ) V ∞ Figure 3.12: Comparison of mass-loss rates (upper panel) and terminal velocities (lower panel) as afunction of the effective temperature. Blue stars correspond to results from this work, black disks toBouret et al. (2005) and Markova et al. (2018) results, and red triangles to theoretical values fromVink et al. (2000). The same colour code but with modified symbols (inverted blue stars, unfilledblack circles and inverted red triangles) are used to represent Markova’s stars with the same effectivetemperature but higher surface gravity. .6. CONCLUSIONS FOR SELF-CONSISTENT M-CAK SOLUTIONS M ,these good results are strongly dependent on the assumed stellar parameters. In the present Chapter we have presented a method to calculate a self-consistent line-forceparameters coupled with the hydrodynamics in the frame of the radiation driven wind theory.Thanks to this procedure, we achieve a unique well-converged solution that does not depend onthe chosen initial values. This is important because it reduces the number of free parameters(now β , v ∞ and ˙ M are no more input parameters) to be determined by fitting syntheticspectra against observed ones.Our calculations contemplate the contribution to the line-force multiplier from more than ∼ 900 000 atomic transitions, a NLTE radiation flux from the photosphere and a quasi-LTEapproximation for the occupational numbers. We have to notice that for T eff > 30 000 K theline force parameters can be confidently used as constants throughout the wind.The set of solutions given in Table 3.3 differs from previous line-force parameter calcu-lations performed by Abbott (1982) and Noebauer & Sim (2015). With these new values,we found a different scale relation for the terminal velocity that is steeper than the usuallyaccepted one. This new relation might explain the observed scatter found in the terminalvelocity from massive stars located at the hot side of the bi-stability jump (Markova & Puls,2008).Concerning the wind parameters derived from modelling O-type stars with homogeneouswinds, our mass-loss rates are in better agreement with the predicted ones given by Vink etal. (2000) formula.For the calculation of synthetic spectra for O-type stars ( ζ -Puppis, HD 163758 and HD164794), we conclude that our procedure’s values for mass-loss rate and hydrodynamics re-produce the observed line-profiles when an adequate value for the clumping factor is chosen.Even knowing the limitations of the m-CAK theory, this remains an extremely usefulframework to get an approach about the real parameters of stellar winds on massive stars. Inspite of the approximations assumed under this theory, we obtain reliable values for mass-lossrates and self-consistent hydrodynamics in a short period of time with a great CPU timesaving (compare with big efforts made by, e.g., Mokiem et al. 2005 or Fierro-Santill´an et al.2018).Our new self-consistent procedure can be used to derive accurate mass-loss rates and:i. Build evolutionary tracks, where a high precision on terminal velocities is not required.ii. Derive truly clumping factors via line-profile fittings.2 CHAPTER 3. SOLUTIONS IN THE FRAME OF M-CAK THEORY T eff log g R ∗ /R (cid:12) Z / Z (cid:12) log t in log t out k α δ v SC ∞ ˙ M SC ˙ M SC / ˙ M Vink [kK] [km s − ] [10 − M (cid:12) yr − ]45 4.0 12.0 1.0 − . − . 53 0.167 0.600 0.021 3 432 ± 240 2 . ± . . − . − . 85 0.142 0.493 0.017 2 329 ± 210 0 . ± . . . − . 07 0.135 0.648 0.022 3 250 ± 300 6 . ± . . − . − . 28 0.114 0.545 0.014 2 221 ± 230 1 . ± . . − . − . 36 0.137 0.629 0.027 3 235 ± 300 3 . ± . . − . − . 73 0.108 0.534 0.019 2 313 ± 230 0 . ± . . − . 80 0.122 0.671 0.039 2 738 ± 230 11 ± . . − . 09 0.091 0.586 0.022 2 043 ± 200 3 . ± . . − . − . 97 0.164 0.581 0.027 3 300 ± 220 0 . ± . . − . − . 44 0.133 0.492 0.038 2 329 ± 160 0 . ± . . . − . 96 0.118 0.659 0.044 2 813 ± 290 6 . ± . . − . − . 40 0.091 0.572 0.025 2 116 ± 230 1 . ± . . . − . 14 0.099 0.715 0.094 1 548 ± 240 14 . ± . . . − . 50 0.073 0.650 0.047 1 334 ± 230 4 . ± . . − . − . 79 0.130 0.610 0.036 3 153 ± 240 1 . ± . . − . − . 28 0.091 0.542 0.033 2 473 ± 300 0 . ± . . − . − . 50 0.132 0.580 0.036 3 314 ± 200 0 . ± . . − . − . 97 0.101 0.517 0.068 2 402 ± 140 0 . ± . . − . − . 55 0.104 0.644 0.062 2 809 ± 240 2 . ± . . − . − . 09 0.071 0.581 0.033 2 534 ± 220 0 . ± . . . − . 77 0.091 0.686 0.116 1 708 ± 170 4 . ± . . . − . 21 0.072 0.607 0.048 1 566 ± 160 1 . ± . . − . − . 37 0.103 0.604 0.043 3 093 ± 210 0 . ± . . − . − . 94 0.069 0.555 0.028 2 791 ± 180 0 . ± . . − . − . 82 0.095 0.637 0.074 2 732 ± 180 1 . ± . . − . − . 46 0.058 0.590 0.031 2 642 ± 180 0 . ± . . . − . 30 0.078 0.675 0.159 1 653 ± 190 1 . ± . . − . − . 15 0.053 0.610 0.052 1 847 ± 140 0 . ± . . Table 3.3: Self-consistent line-force parameters ( k, α, δ ) for adopted standard stellar parameters,together with the resulting terminal velocities and mass-loss rates ( ˙ M SC ). Ratios between self-consistent mass-loss rates and Vink’s recipe values (re-scaled to match metallicity from Asplundet al., 2009) using v ∞ /v esc = 2 . ± 500 for effective temperature, ± . 05 for logarithmof surface gravity and ± . .6. CONCLUSIONS FOR SELF-CONSISTENT M-CAK SOLUTIONS T eff log g log t in log t out k α δ | ∆ v ∞ | | ∆ ˙ M | [km s − ] [10 − M (cid:12) yr − ]45 000 4.0 − . − . 03 0.099 0.686 0.037 780 0 . − . − . 87 0.107 0.650 0.029 600 0 . − . − . 71 0.120 0.638 0.027 420 0 . − . − . 55 0.167 0.600 0.021 0 040 000 4.0 − . − . 50 0.099 0.633 0.040 521 0 . − . − . 32 0.099 0.634 0.036 610 0 . − . − . 14 0.107 0.621 0.026 594 0 . − . − . 96 0.164 0.581 0.027 0 040 000 3.6 0.08 − . 44 0.095 0.666 0.090 247 0 . − . 28 0.098 0.680 0.075 75 0 . − . 12 0.101 0.692 0.067 323 0 . − . 96 0.118 0.659 0.044 0 036 000 3.6 − . − . 00 0.084 0.637 0.112 520 0 . − . − . 85 0.092 0.648 0.078 114 0 . − . − . 70 0.089 0.668 0.075 267 0 . − . − . 55 0.104 0.644 0.062 0 032 000 3.4 0 . − . 49 0.066 0.630 0.251 631 0 . . − . 43 0.075 0.636 0.221 457 0 . . − . 37 0.079 0.662 0.179 168 0 . . − . 31 0.078 0.675 0.159 0 0 Table 3.4: Influence of the optical depth interval on the line-force parameters for some referencemodels given in Table 3.3. Absolute values of the differences in the resulting wind parameters withrespect to the reference solution are presented. previous studies present workReference T eff log g R ∗ /R (cid:12) ˙ M v ∞ k α δ ˙ M SC v SC ∞ [kK] [10 − M (cid:12) yr − ] [km s − ] [10 − M (cid:12) yr − ] [km s − ] Noebauer & Sim (2015) 42 3.6 19.0 45 . . ± . . ± Bouret et al. (2012) 40 3.64 18.7 2 . . ± . . ± Puls et al. (2006) 39 3.6 29.7 8 . . ± . . ± . . ± . . ± Table 3.5: Stellar and wind parameters for ζ -Puppis from previous studies compared with ourself-consistent results. Line-force parameters are also listed. previous studies present workName T eff log g R ∗ /R (cid:12) ˙ M v ∞ k α δ ˙ M SC v SC ∞ [kK] [10 − M (cid:12) yr − ] [km s − ] [10 − M (cid:12) yr − ] [km s − ]HD 163758 34.5 3.41 21.0 1.6 2 100 0.087 0.679 0.112 3 . ± . . ± . ± . . ± Table 3.6: Idem Table 3.5, but for HD 163758 and HD 164794. Stellar and wind parameters arefrom Bouret et al. (2012) and Krtiˇcka et al. (2015) respectively. CHAPTER 3. SOLUTIONS IN THE FRAME OF M-CAK THEORY Model T eff log g R ∗ /R (cid:12) k α δ v ∞ ˙ M [kK] [km s − ] [10 − M (cid:12) yr − ]Self-Consistent 43.5 4.0 11.9 0.159 0.603 0.032 3 342 ± 240 1 . ± . . β = 1 . ± 290 1 . ± . . Self-Consistent 39 3.6 19.45 0.116 0.657 0.079 2 412 ± 210 5 . ± . . β = 0 . ± 570 4 . ± . . Table 3.7: Comparison of self-consistent with β –law (single-step) models for the two stars analyzedby Bouret et al. (2005). Self-consistent models reproduce better the line-fitted wind parametersobtained by these authors ( β =1: v ∞ = 3000 km s − , ˙ M = 1 . × − M (cid:12) yr − , and β =0 . v ∞ = 2300 km s − , ˙ M = 6 × − M (cid:12) yr − ). Field star T eff log g R ∗ /R (cid:12) k α δ v SC ∞ ˙ M ˙ M SC ˙ M obs ˙ M SC ˙ M Vink [kK] [km s − ] [10 − M (cid:12) yr − ]HD 169582 37 3.5 27.2 0.102 0.668 0.063 3 017 ± 700 7 . ± . . ± 540 1 . ± . . ± 470 0 . ± . . ± 580 3 . ± . . ± 460 0 . ± . . ± 350 5 . ± . . ± 300 3 . ± . . ± 300 0 . ± . . Table 3.8: Resulting self-consistent wind parameters ( v SC ∞ and ˙ M SC ) calculated for stars analyzedby Markova et al. (2018). Error margins presented here for wind parameters are undergone fromuncertainties of ± T eff and ± . g . Last two columns show the ratio between self-consistent and observed mass-loss rates and the ratio between self-consistent and Vink’s mass-lossrates. hapter 4 Self-consistent Solutions Under W -Lambert Procedure During the previous two Chapters, we have focused our analysis on self-consistent solutionsfor line-driving winds based on m-CAK theory, which has demonstrated to be a fast andconfident enough prescription to give theoretical values for wind parameters, specially formass-loss rates. However, self-consistent solutions calculated by Alfakdelta27 are basedon two main assumptions:i) Sobolev approximation which is a fundamental part of m-CAK theory andii) more important, the quasi-NLTE treatment given by atomic populations.Despite the fact that quasi-NLTE scenario provided us reliable results, a more completeanalysis is required for a complete NLTE treatment in stellar winds to make a parallel betweenboth approaches.The goal of this chapter will be finding a self-consistent solution for the stellar winds usingthe radiative transfer CMFGEN (Hillier, 1990a,b; Hillier & Miller, 1998) as a tool. This code,together with providing us with a synthetic spectrum for a specific set of stellar and windparameters taking into account all the statistical equilibrium relationships between all theatomic populations existing in the atmosphere, it also provides us with an output for radiativeacceleration g rad calculated from the solution of the radiative transfer equation. Since thisoutput can be expressed as a function of radius only g rad ( r ), it can be adopted later to solveEq. 2.2 (equation of momentum) analytically using the so-called W -Lambert equation. Bymeans of introducing this new velocity field calculated from equation of momentum to re-execute CMFGEN, it is possible to perform again an iterative procedure capable of reaching anew self-consistent solution, this time with the full NLTE line-acceleration given by CMFGENand the analytical solution for the equation of momentum. The price to pay in this case is abigger computational effort and hence more time, so the analysis made during this Chapteris limited only to the well-known star ζ -Puppis. Another price to pay is the fact that, as we656 CHAPTER 4. SOLUTIONS UNDER LAMBERT PROCEDURE mentioned in Section 2.1, since Eq. 2.8 does not have an implicit dependence on density (andthen mass-loss rate) and line-acceleration is obtained directly from the output of CMFGEN,mass-loss rate is an extra free parameter to be set, together with effective temperature, surfacegravity and abundances. Nevertheless, results presented here could be the basis for extendedstudies on other massive stars.The discussion of the results and analysis done in this Chapter will be compared withthe study implemented by Sander et al. (2017), who performed a similar procedure for a fullNLTE self-consistent solution using the radiative transfer code PoWR (Gr¨afener et al., 2002;Hamann & Gr¨afener, 2003). And later, the comparison will focus between this full NLTEtreatment with the quasi-NLTE performed for Gormaz-Matamala et al. (2019). W -Lambert equation As we previously mentioned, equation of momentum (Eq. 2.2) can be solved analyticallywith the help of a mathematical tool called the Lambert W -function (also named productlogarithm ), which is defined by the inverse of the function: z ( w ) = we w , (4.1)with z being any complex number. That means: W ( z ) e W ( z ) = z , (4.2)(Corless et al., 1993, 1996).Lambert W -function presents infinite solutions depending of all the possible non-zerovalues that z may take. If we limit our search only for real values of z , we find that W ( z )function is multivalued because f ( w ) = we w is not injective. Because of this reason, we splitLambert function into two sections, corresponding them to the branches W for W ( z ) ≥ − W − for W ( z ) ≤ − z : W (cid:18) − e (cid:19) = W − (cid:18) − e (cid:19) = − . (4.3)From Eq. 4.1 and Figure 4.1, it is clearly seen that the domains for z and the codomainsfor the two branches are : z ∈ (cid:20) − e , (cid:20) for W − ∈ [ − , −∞ [ , (4.4) z ∈ (cid:20) − e , + ∞ (cid:20) for W ∈ [ − , + ∞ [ . (4.5) Notation for open and closed intervals following definition given by the Encyclopaedia of Mathematics: .2. LAMBERT-PROCEDURE - - - 20 z W ( z ) Figure 4.1: Real branches of W k ( z )-Lambert function: k = 0 (dashed blue) and k = − This means, branch W − diverges when z approaches zero or, in other words, z is anasymptote when W − tends to −∞ . Because of this condition, W − function can be used tosolve any function with the same asymptotic behaviour as z , such as the velocity field of astellar wind. Then, the analytical solution for the equation of momentum with g line ( r ) given isobtained once we reformulate Eq. 2.2 in terms of the Lambert function, as it is demonstratedin Section 4.2.2. To evaluate whether Lambert-procedure can be successfully applied into hydrodynamic mod-els, we take as starting point a CMFGEN model for ζ -Puppis with the parameters shown inTable 4.1. These both (stellar and wind) parameters are chosen because they were used to fit ζ -Puppis (Bouret et al., 2012, Section 6.5) together with being also determined by Marcolinoet al. (2017) on their spectral analysis in infrared.In the previously mentioned Section 6.5 of Bouret et al. (2012), they also check the consis-tency of the line-force. Radiative acceleration for this initial model is plotted in Fig. 4.2. This g rad ( r ) was internally calculated by CMFGEN starting from the initial parameters shown inTable 4.1. However, the values tabulated there do not lead to recover the same wind param-eters; hence, it is not a self-consistent result. A truly self-consistent solution must satisfy theequation of momentum: g rad = v dvdr + 1 ρ dpdr + GM ∗ r . (4.6)8 CHAPTER 4. SOLUTIONS UNDER LAMBERT PROCEDURE T eff 41 000 Klog g R ∗ R (cid:12) ˙ M . × − M (cid:12) yr − v ∞ − β f ∞ Stellar and wind parameters of our initial CMFGEN model. R ∗ is considered for opticaldepth τ = 2 / 3, whereas parameter f ∞ corresponds to the infinite filling factor (see Section 4.3.2 fordetails). ( r / R * ) g r a d Figure 4.2: Radiative acceleration (arbitrary units) as a function on radius for the initial CMFGENmodel with β –law. Mass-loss rate M ∞ is linked with density and velocity profiles by means of the equationof continuity (Eq. 2.3)Figure 4.3 shows both left-hand and right-hand sides of Eq. 4.6 for the initial CMFGENmodel tabulated in Table 4.1, where it is clearly seen that they do not match. This discrepancyyields in the fact that velocity field is not calculated from the line-acceleration itself, but itis consequence of using a β –law. It would be possible to argue that this discrepancy mightbe a specific situation only, but it was previously noticed also by Bouret et al. (2012); hence,the lacking of consistency between hydrodynamics and radiative acceleration for CMFGENseems to be a general rule.As a consequence of this, we must perform a prescription capable to equalise both sidesof Eq. 4.2. Following Lambert-procedure described below in this section is done aiming to .2. LAMBERT-PROCEDURE g rad / g elec ( vdv / dr + ρ - dP g / dr + g )/ g elec ( r ) g / g e l ec Figure 4.3: Comparison of left-hand side of Eq. 4.6 (blue) and right-hand side (orange) for thehydrodynamics obtained from the initial CMFGEN model (Table 4.1). Accelerations were rescaledby g elec (Eq. 4.13) for illustrative reasons. evaluate:i) whether there is an hydrodynamic solution capable to properly couple line-accelerationand velocity/density fields,ii) whether this solution is stable or not andiii) whether this solution does reproduce a spectrum in agreement with the observations. g line ( r ) Despite the fact that line-acceleration used to solve equation of momentum is an outputof a converged CMFGEN model, it is important to describe in general terms how is thiscalculated inside the code. For that purpose, it is necessary to do a brief summary aboutradiative transfer, and how this is related with the resulting g line ( r ). To avoid confusions,we remark the fact that what CMFGEN gives us exactly, corresponds to the total radiativeacceleration g rad , i.e., the sum of both the acceleration produced by the line-driving processand that one produced by the continuum by means of Thomson scattering. However, thesetwo terms are easily separated, as it is shown in Eq. 2.5.Radiative transfer equation, which describes the gaining or losing of radiative energytrough a path due to emission and absorption, is given by the formula: dI ν ds = − χ ν I ν + η ν , (4.7)with χ ν = κ ν ρ being the opacity, given κ ν as defined in Eq. 2.14, and η ν the emissivity and0 CHAPTER 4. SOLUTIONS UNDER LAMBERT PROCEDURE ds is an infinitesimal element of path (Hillier, 1990a). Using the differential form of opticaldepth (Eq. 2.18), dτ ν = χ ν ds = κ ν ρds , we rewrite radiative transfer equation as: dI ν χ ν ds = − I ν + η ν χ ν ,dI ν dτ ν = − I ν + S ν , (4.9)with S ν being the source function , the ratio between the emission and absorption coefficients.Notice that the formal solution of this equation is: I ν ( τ ν ) = I ν (0) e − τ ν + (cid:90) τ ν e τ (cid:48) ν − τ ν S ν ( τ (cid:48) ν ) dτ (cid:48) ν , (4.10)which is the general expression for Eq. 2.27 where, thanks to the Sobolev approximation, wewere focused in the absorption only and then we neglected possible any emission source.Full solution for the intensity of radiation then, depends on the value of the source function S ν which in turn depends on the case if we are assuming LTE conditions or not. For localthermodynamic equilibrium, it is fulfilled that the source function is equal to the Planckfunction B ν (Carroll & Ostlie, 1996, Section 9.4). However, for expanding atmospheres thisis not the case, and thus S ν must be calculated taking into account all the atomic processesamong the different ions (collisions, recombinations, bound-free and bound-bound transitions,etc.). Since many of these interactions depend on the intensity of radiation, this becomes acoupled problem and then codes such as CMFGEN, FastWind or PoWR are necessaryto solve the statistical equilibrium equations. Henceforth, solving this coupled problem forthe radiative transfer it is possible later to calculate the flux mean opacity by means of theintegration over all the frequencies: 1¯ χ = (cid:82) ∞ χ ν ∂S ν ∂T dν (cid:82) ∞ ∂S ν ∂T dν . (4.11)Radiative acceleration is then evaluated using the formula: g rad ( r ) = ¯ χ f L ∗ πcρr , (4.12)where χ f is the flux mean opacity.This radiative acceleration given by CMFGEN corresponds to the total acceleration dueto radiative processes, i.e., it considers not only the effects of absorption and reemission ofphotons by line transitions, but also electron scattering. However, acceleration by electronscattering is implicitly included in the momentum equation by means of the Eddington factorΓ e . g elec = GM Γ e r . (4.13) We have presented here the simplest expression for the differential radiative transfer equation (i.e., alonga simple segment ds ) because of didactic reasons. The full expression introduced by Mihalas et al. (1975) fora spherical coordinate system is: µ ∂∂r I ν + 1 − µ r ∂∂µ I ν − ν v ( r ) cr (cid:20) − µ + µ (cid:18) d ln vd ln r (cid:19)(cid:21) ∂∂ν I ν = − χ ν I ν + η ν . (4.8)In this case, the variables χ ν , η ν and I ν are in function of the radius r and the angle µ = cos θ . .2. LAMBERT-PROCEDURE g line = g rad − g elec . (4.14) The equation of momentum for a stationary, one-dimensional, non-rotating, isothermal, out-flowing wind in spherical coordinates is given by: v dvdr = − ρ dpdr − GM ∗ (1 − Γ e ) r + g line ( r ) , (4.15)where p is the gas pressure, v is the wind velocity and GM ∗ (1 − Γ e ) /r is the gravitationaleffective acceleration.According with M¨uller & Vink (2008) and Araya et al. (2014), Eq. 2.8 can be expressedin a dimensionless way by making the following change of variables:ˆ r = rR ∗ , ˆ v = va and ˆ v crit = v esc a √ a (cid:115) GM ∗ (1 − Γ e ) R ∗ , (4.16)being a the isothermal speed of sound given by: a = k B T eff µm H + 12 v (4.17)This formula differs from those used by Sander et al. (2017), because we are consideringtemperature field T ( r ) and mean particle µ ( r ) as constants. Micro-turbulence velocity v mic is included, which is using to be in the order of ∼ 10 km s − . After that, defining thedimensionless line acceleration: ˆ g line ( r ) = R ∗ a g line ( r ) , (4.18)the equation of momentum reads:ˆ v d ˆ vd ˆ r = − ρ dpd ˆ r − ˆ v ˆ r + ˆ g line (ˆ r ) . (4.19)With the use of equation of state for an ideal gas ( p = a ρ ), Eq. 4.19 becomes independentof density (and therefore independent of mass-loss rate) equation of motion: (cid:18) ˆ v − v (cid:19) d ˆ vd ˆ r = − ˆ v ˆ r + 2ˆ r + ˆ g line (ˆ r ) , (4.20)By the integration of both sides along the atmosphere, Eq. 4.20 can be solved analytically inorder to obtain ˆ v (ˆ r ). Subsonic v/s supersonic region Assuming a monotonic behaviour of the velocity field throughout the atmosphere ( d ˆ v/d ˆ r > v = a . This condition is fulfilled at the sonic (or critical) point ˆ r c : − ˆ v ˆ r c + 2ˆ r c + ˆ g line (ˆ r c ) = 0 , (4.21)2 CHAPTER 4. SOLUTIONS UNDER LAMBERT PROCEDURE Because of the monotonic behaviour of v ( r ) (and therefore ˆ v ), the sonic point becomesa boundary between two regions. The first of them is where ˆ v < subsonicregion . The second one, where ˆ v > 1, is then called the supersonic region . Differentiationbetween both is not only made due to mathematical analysis but also physical reasons. In aone dimensional fluid moving at velocities below sound speed perturbations are propagatedboth inwards and outwards, whereas perturbations occurred on a supersonic fluid only arepropagated in the direction of the flux.This means that we actually are searching two solution branches that merged at the criticalpoint. If we focus at first in the supersonic region, going from the sonic point towards infiniteand where velocity field scales from sound speed to the asymptotic terminal velocity, theequation to be solved by Lambert W -function (branch W − ) is obtained by the integrationof Eq. 4.20. ˆ v e − ˆ v = − (cid:18) ˆ r c ˆ r (cid:19) exp (cid:20) − − v (cid:18) r − r c (cid:19)(cid:21) × exp (cid:34) − (cid:90) ˆ r ˆ r c ˆ g line d ˆ r (cid:35) , with ˆ r > ˆ r c . (4.22)Therefore, if we define the right-hand side of Eq. 4.22 as x (ˆ r ), expression for velocity fieldis directly calculated from the Lambert W -function:ˆ v (ˆ r ) = (cid:113) − W j ( x (ˆ r )) , (4.23)with j being the branch of the Lambert W -function (0 or − W -function, which are alsomerged in a specific point (see Fig. 4.1), it would be expectable to obtain an analyticalexpression for the subsonic region using the branch W . However, when this idea was executedon CMFGEN big errors were produced for the acceleration and velocity profiles on the lowestpart of the wind (hydrostatic region, below the photosphere). Therefore, we decide to useLambert-procedure to calculate hydrodynamics above the sonic point only, making later agood coupling between both supersonic and subsonic regions without introduce big errors inthe hydrostatic part. This coupling must ensure a smooth transition between subsonic andsupersonic zones, reason why those conditions must be fulfilled:lim r → a − v ( r ) = lim r → a + v ( r ) , (4.24)and lim r → a dv ( r ) dr (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) ra . (4.25)The second limit is easier to evaluate if we use the logarithm derivatives by means of theequivalence: d ln vd ln r = rv ( r ) dvdr . (4.26) .2. LAMBERT-PROCEDURE v ( r ) - - - - - - R * / r ( d l n v ) / ( d l n r ) Figure 4.4: Zoom over the sonic point (red point), where subsonic and supersonic regions arecoupled during the Lambert-procedure. Magenta solid line represents the new velocity profile for thesupersonic region obtained with Lambert W -function, whereas blue dashed line represents the oldvelocity profile for the subsonic region. Then: lim r → a d ln v ( r ) d ln r (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) ra . (4.27)As consequence of these conditions, subsonic region must be readapted too, in functionof saving a smooth transition and avoid instabilities. Fig. 4.4 shows the coupling of v ( r ) and d ln v/d ln r around the sonic point for a velocity profile after solving hydrodynamics withLambert W -function. The gaps displayed are produced because the sonic point given byfinding the root of Eq. 4.21 may not be equal to the old v ( r c ) = a given by the input velocityprofile (differences on these two sonic points were also previously referred by Sander et al.,2017, see their Fig. 1), and because Lambert W -function may produce a more or less steepwind close to the sonic point. In any case, our self-consistent solution must be rescaled in the4 CHAPTER 4. SOLUTIONS UNDER LAMBERT PROCEDURE Initial CMFGEN modelStellar parameters: T e↵ , log g , R ⇤ Wind parameters: ˙ M , f -law velocity profileNO Running CMFGEN g line (r)satisfies convergence criteria? v ( r ), ⇢ ( r ) Lambert W -FunctionYES Lambert-convergedCMFGEN modelCMFGEN synthetic spectrum Figure 4.5: Flowchart of the Lambert-procedure. Notice that stellar parameters are the same duringthe entire procedure, as well as mass-loss rate and clumping factor. subsonic region in order to cancel these gaps.Aiming this, it is important an accurate determination of features on the inner part of thewind, such as photospheric radius and the sonic point location itself. We do not use Lambert W -function to obtain a new velocity field for the subsonic region because CMFGEN is toosensitive to modifications close to the photosphere (the hydrostatic part of the wind). Steepmodifications affect largely features such as opacity and atomic populations, which leads tolarge errors in the re-calculation of line-acceleration. That is the main reason why we applyLambert W -function only to supersonic region, whereas wind below sonic point is slightlymodified by rescaling the velocity profile in order to satisfy accurately Eq. 4.24 and Eq. 4.27. The convergence of the Lambert-iterations can be checked by evaluating both either velocityfield v ( r ) or line-acceleration g line ( r ). In order to accept a Lambert-model as well converged,we impose as condition to satisfy the following relationships: (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) log v ( r ) n v ( r ) n − (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) ≤ . , (4.28)and: (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) log g line ( r ) n g line ( r ) n − (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) ≤ . , (4.29)where n is the n -th Lambert-iteration executed . This condition is applied for line-accelerationtoo. Reason to choice this value as threshold yields in the fact that no significative differences Notice that, if f ( r ) n = f ( r ) n − the value of the absolute expression would be zero, because the ratiobetween both functions would be 1. .3. RESULTS v ∞ (km s − ) β T41blaw01 . T41blaw02 . T41blaw03 . T41blaw04 . T41blaw05 . Hydrodynamic parameters of initial models. on resulting spectra are observed between models with hydrodynamic differences smaller thatthis value.The threshold previously presented is applied not only for external part of the wind (i.e.,terminal velocity v ∞ ) but for the entire range on radius r from the photosphere outwards.This, because the part which is most sensitive to changes is where the coupling is done (aroundthe sonic point). A good convergence around this zone allows us the obtention of a stableLambert-model with a smooth transition between subsonic and supersonic regions.Later, existence of a unique solution must be confirmed. Wind hydrodynamics expressedunder β –law depends on two parameters: terminal velocity v ∞ and the β exponent. If theLambert-procedure is able to generate a well converged hydrodynamic solution, it should bethe same for all the initial velocity profiles (i.e., initial v ∞ and β parameters) chosen.Table 4.2 shows hydrodynamic parameters for each initial model. All these models havethe same value for mass-loss rate ( M = 2 . × − M (cid:12) yr − ) and clumping filling factor( f ∞ = 0 . M and f ∞ are implemented.A scheme of Lambert-procedure is presented in Fig. 4.5. Typically, Lambert-procedureconverges after 5 or 6 Lambert-iterations, each one of them corresponding to a CMFGENmodel executed with 30 inner iterations. Through this section we will analyse the convergence of the Lambert procedure. First partis focused on the results starting from the initial model with the parameters given in Ta-ble 4.1, together with the check of the convergence of the alternative models tabulated inTable 4.2 over the same final solution. Second part is a more extensive analysis evaluatingthe consequences of using different initial values for mass-loss rate and clumping. Nomenclature for the name CMFGEN models executed in this procedure is as follows: T41 means a modelwith T eff = 41 kK, blaw means a model with β –law whereas lamb means a converged Lambert-model. CHAPTER 4. SOLUTIONS UNDER LAMBERT PROCEDURE - - - - - - R * / r v ( r ) - - - - - - R * / r v ( r ) Figure 4.6: Differences between velocity field with β –law (blue) and the resulting after Lambert-iterations (orange) for a CMFGEN model with ˙ M = 2 . × − and f ∞ = 0 . 1. A zoom over thetransition region is included, displaying the sonic point with a red dot. - - - - - R * / r g li n e Figure 4.7: Same as Fig. 4.6, but comparing line-accelerations instead velocity fields. Applying Lambert-procedure over the starting CMFGEN model with the parameters shownin Table 4.1, a new CMFGEN model has been created with a new hydrodynamics and a newline-acceleration. These both new features are related each other by means of the equationspreviously discussed in Section 4.2, i.e., wind hydrodynamics and line-acceleration are self-consistent between them. Comparisons between initial and converged self-consistent velocityfields are shown in Fig. 4.6, whereas comparisons between line-accelerations are shown inFig. 4.7. Both features were evaluated to satisfy the threshold condition from Eq. 4.28.Because this time the model is self-consistent, both left-hand and right-hand sides of Eq. 4.6 .3. RESULTS g rad / g elec ( vdv / dr + ρ - dP g / dr + g )/ g elec ( r ) g / g e l ec Figure 4.8: Radiative acceleration (divided by acceleration due to electron scattering) as a functionon wind velocity for the final CMFGEN Lambert-model (blue), compared with the expected valuewhen equation of momentum is solved (orange). are in agreement as it is observed in Fig. 4.8 (compare with previous Fig. 4.3).Concerning to the converged self-consistent new hydrodynamic, as first comment we ob-serve that resulting terminal velocity has increased in a factor of ∼ . − . Resulting line-acceleration is also higher for the Lambert-model, but this result maybe consequence on the chosen initial parameters ( T eff , log g , ˙ M ) only. Besides, due to therescaling in the subsonic region, there are not significant differences in the resulting velocityprofile below the sonic point, just above ∼ − 13 km s − . As we previously stated in Sec-tion 4.2, this is the most effective method to ensure a good coupling between two zones, andalso it estabilised the subsonic regions (which is very sensitive to sharp changes in CMFGEN).Moreover, it is interesting the fact that hydrodynamics can be characterised not onlywith a new terminal velocity, but also they can be fairly approximated with a new β factor(which may differ from the initially set as input). This effect is clearly seen in Fig. 4.9, wherevelocity profile with β = 0 . v ( r ), although that does not meanthat Lambert procedure gives a new value for β , because resulting hydrodynamics are notfitted but calculated from respective line-acceleration. However, it is possible to find a β value (hereafter quasi- β ) capable to closely fit the new velocity field. This approximate fitwith a β -law is expectable because the wind of ζ -Puppis lies in the range of the so-called fastsolutions (Cur´e et al., 2012; Gormaz-Matamala et al., 2019), but we cannot assure that wecould find the same quasi- β for stars at lower temperatures with δ -slow solutions (Cur´e et al.,2011).8 CHAPTER 4. SOLUTIONS UNDER LAMBERT PROCEDURE - - - - - - R * / r v ( r ) - - - - - - R * / r v ( r ) Figure 4.9: Velocity profile of the converged Lambert-model (red), fitted with a β -law using β = 0 . β = 0 . 95 (green) and β = 0 . .3. RESULTS - - - - - R * / r v ( r ) - - - - - R * / r v ( r ) Figure 4.10: Left panel: velocity profiles (in units of km s − ) corresponding to the models shownin Table 4.2: T41blaw01 (red), T41blaw02 (green), T41blaw03 (cyan), T41blaw04 (blue), T41blaw05 (magenta). Right panel: converged Lambert-profiles for each initial model. To be sure that this well converged Lambert-solution is independent on initial values for β and v ∞ , we proceed to repeat the process for every model signalled in Table 4.2, and whosegraphs are shown in Fig. 4.10: left panel represents the initial velocity profiles with β -law,whereas right panel shows the final converged Lambert velocity profiles (where the overlapdemonstrates all the initial models converging into the same solution). Previously, when we transform initial momentum equation, Eq. 2.8, into Eq. 4.20, the explicitterm for ρ ( r ) was replaced by v ( r ) by means of the equation of continuity, Eq. 2.3, and thuswe obtained a formula independent on density. As consequence, Lambert-procedure cannotcalculate a new mass-loss rate because ˙ M must be set as an input. For this reason, this windparameter is considered as free when we are calculating our Lambert-procedure.However, the dependence on density, and therefore the mass-loss rate, is implicitly in-cluded in Eq. 4.20 inside the g line term. From Eq. 4.12 remains clear that line-accelerationis inversely proportional to mass-loss rate ˙ M , and hence also inversely proportional to thedensity ρ ( r ). Because this density (implicitly included inside line-acceleration) is directly pro-portional to ˙ M which keeps constant during Lambert-procedure, these differences in g line ( r )leads to different final Lambert-hydrodynamics. Higher mass-loss rates produces slower line-accelerations (Fig. 4.11), which undergoes into velocity profiles with lower terminal velocities.Thus, the mass-loss rate is not the only wind parameter to be considered as free. Line-acceleration is affected not only by the general matter density but also for the small scaleinhomogeneities (i.e., clumping) present through the wind. Clumping is implemented byCMFGEN in terms of the volume filling factor f , which assumes a void interclump medium0 CHAPTER 4. SOLUTIONS UNDER LAMBERT PROCEDURE - - - - - - - - - R * / r g li n e Figure 4.11: Line-accelerations from initial CMFGEN models with different values for ˙ M . Everyother stellar and wind parameter is the same as in Table 4.1. and the clumps to be small compared to the photons mean free path. The filling factor issuch that ρ = ρ /f , where ρ is the homogeneous (unclumped) wind density (Bouret et al.,2005). Volume filling factor is defined in terms on the velocity field: f ( v ( r )) = f ∞ + (1 − f ∞ ) e − v ( r ) /v cl . (4.30)Usually, clumping can be expressed by the [infinite] filling factor f ∞ only. Any strongerclumping factor, expressed with a smaller filling factor f ∞ , gives a stronger line-acceleration(Fig. 4.12), which will also lead into faster hydrodynamics with higher terminal velocities.This effect is not a consequence derived from Eq. 4.12 because density profile keeps unchangedin all these cases. We could argue that the presence of the overdensities (which produce theoverestimation for mass-loss rates) is the responsible of the increasing on g line but as wepointed out previously, line-acceleration is inversely proportional to ρ ( r ) and hence withhigher clumping factor g line should be smaller. Therefore, the reason why line-accelerationbecomes higher when clump inhomogeneities intensifies remains not completely clear.Since neither mass-loss rate nor filling factor are altered when Lambert-iterations areexecuted, different combinations of ( ˙ M , f ∞ ) will lead into different converged new Lambert-hydrodynamics, with their corresponding new terminal velocities. Some of these new self-consistent hydrodynamics generated by Lambert-iterations are shown in Table 4.3. Thehomogeneous equivalence for mass-loss rate ( ˙ M / √ f ∞ ) is shown in order to make clearer theinfluence of chosen clumping into the resulting Lambert-hydrodynamics. The most importantresult derived from Table 4.3 (even when it could be considered an obvious consequence), isthe direct relationship between line-acceleration and the resulting terminal velocity for eachself-consistent solution. The higher g line , the larger final v ∞ , and therefore a self-consistent .3. RESULTS f ∞ = f ∞ = f ∞ = - - - - - - R * / r g li n e Figure 4.12: Line-accelerations from initial CMFGEN models with different values for f ∞ . Everyother stellar and wind parameter is the same as in Table 4.1. solution with an specific terminal velocity can be obtained using the adequate set of ˙ M and f ∞ .However, it is important to emphasise that this family of numerically well convergedLambert-hydrodynamics given different combinations for ˙ M and f ∞ are not valid for any possible mass-loss rate, but it is constrained within a range of plausible values. This comesfrom the fact that the error associated in the equation of momentum for a CMFGEN model,defined as: Err( r ) = 200 (cid:16) v dvdr + ρ dPdr + Pρ dTdr − g tot (cid:17)(cid:12)(cid:12)(cid:12)(cid:12) v dvdr (cid:12)(cid:12)(cid:12)(cid:12) + (cid:12)(cid:12)(cid:12)(cid:12) ρ dPdr (cid:12)(cid:12)(cid:12)(cid:12) + (cid:12)(cid:12)(cid:12)(cid:12) Pρ dTdr (cid:12)(cid:12)(cid:12)(cid:12) + | g tot | , (4.31)grows considerably when the initial hydrodynamics introduced departs too far away from theprevious model, specially in the so-called acceleration zone (where velocity field takes valuesaround ∼ − 100 km s − ). For example, when we consider mass-loss rates below 1 . × − M (cid:12) yr − , resulting hydrodynamic is so fast in the supersonic region that does not allow asmooth coupling; hence, we can consider this value as the lower plausible limit for Lambert-solutions. A similar scenario is seen when the introduced mass-loss rate is too large as morethan 6 . × − M (cid:12) yr − : resulting Lambert-hydrodynamics is so slow that is ”compressed”in the supersonic region. The existence of a minimum error should be used as a tracer to findthe most accurate value for the mass-loss rate; nevertheless this is not possible at all, becausethe erratic behaviour of Err( r ) does not allow the performance of a good constraint. Thisimplies, the best free mass-loss rate (and clumping filling factor) must be determined by thespectral fit.2 CHAPTER 4. SOLUTIONS UNDER LAMBERT PROCEDURE Name ˙ M clump f ∞ ˙ M hom v ∞ , Lamb ( M (cid:12) yr − ) ( M (cid:12) yr − ) (km s − ) T41lamb01 . × − . × − T41lamb02 . × − . × − T41lamb03 . × − . × − T41lamb04 . × − . × − T41lamb05 . × − . × − T41lamb06 . × − . × − T41lamb07 . × − . × − T41lamb08 . × − . × − T41lamb09 . × − . × − T41lamb10 . × − . × − T41lamb11 . × − . × − Converged models, given a set of different values for mass-loss rate and clumping fillingfactor for stellar parameters T eff = 41 kK, log g = 3 . R ∗ = 17 . R (cid:12) . Discussion about Lambert-solutions will be centred on their capacity to be combined withspectral fitting analysis in order to derive the best self-consistent wind parameters, and thecomparison of these new Lambert-solutions with previous self-consistent studies (see, e.g.,Sander et al., 2017, and Chapter 3). Despite the find of several plausible self-consistent hydrodynamics for (this case) ζ -Puppis, the real physical solution must be unique, and just one of these solutions can be the accurate one.For that reason, the following step is to evaluate the resulting spectra from all these Lambert-hydrodynamic solutions, checking which combination of ( ˙ M , f ∞ ) provides an accurate fit forthe observed spectrum. Observed spectra for ζ -Puppis (HD 66811) were taken with FEROS for the visible range and with IUE for the ultraviolet, and corresponds to those ones usedpreviously by Bouret et al. (2012).The two main criteria to be satisfied on the spectral fitting are: the blue side of the P-Cygni profile C IV 1548 for the ultraviolet range, and H α for the optical. C IV 1548 is agood indicator of v ∞ in O, B and WR stars (Prinja et al., 1990), and the blue side of thisP-Cygni profile can be reproduced using CMFGEN (Gormaz-Matamala et al., 2015). Because Fiber-fed Extended Range Optical Spectrograph: International Ultraviolet Explorer: https://archive.stsci.edu/iue/ .4. DISCUSSION N o r m a li s ed F l u x λ ( Ang ) Figure 4.13: Fit of T41lamb09 (see Table 4.3) over the observed FEROS spectrum for the visiblerange from 4 000 to 7 000 ˚ A . of the direct proportionality between line-acceleration and self-consistent terminal velocity, itis possible to find a set of ( ˙ M , f ∞ ) capable to fit v ∞ , which has a value around ∼ − . This, together with the constraint coming from the fitting of the H α profile (which isproportional to the homogeneous mass-loss rate, tabulated in Table 4.3) lead us to an uniquecombination of ˙ M and f ∞ capable to satisfy both criteria. Following this, we proceed then toanalyse the spectra of our Lambert-solutions, and it is found that the best fit (done by-eye)corresponds to the model T41lamb09 , with ˙ M = 3 . × − M (cid:12) yr − and f ∞ = 0 . ∼ . CHAPTER 4. SOLUTIONS UNDER LAMBERT PROCEDURE N o r m a li s ed F l u x λ ( Ang ) Figure 4.14: Fit of T41lamb09 (see Table 4.3) over the observed IUE spectrum for the ultravioletrange from 1 200 to 1 800 ˚ A . Name ˙ M clump f ∞ ˙ M hom v ∞ , Lamb ( M (cid:12) yr − ) ( M (cid:12) yr − ) (km s − ) T40lamb01 . × − . × − T40lamb02 . × − . × − T40lamb03 . × − . × − T40lamb04 . × − . × − T40lamb05 . × − . × − Converged models, given a set of different values for mass-loss rate and clumping fillingfactor for stellar parameters T eff = 40 kK, log g = 3 . 64 and R ∗ = 18 . R (cid:12) . using a value for clumping half than before. This implies an almost similar value for the”homogeneous” mass-loss rate ( ˙ M / √ f ∞ ): 7 . × − M (cid:12) yr − for us, 8 . × − M (cid:12) yr − for Bouret et al. (2012). In spite of this result, since this initial Lambert-converged solutionwas obtained starting from a set of stellar parameters different to those finally determinedfor the star, a more deep analysis must be done starting from the stellar parameters used byBouret et al. (2012) to fit the spectra of ζ -Puppis, in order to include a comparison of spectralfits, i.e., new initial stellar parameters are T eff = 40 kK, log g = 3 . 64 and R ∗ = 18 . R (cid:12) . .4. DISCUSSION N o r m a li se dF l u x λ ( Ang ) Figure 4.15: Fit of T40lamb04 from Table 4.4 (blue solid line) over the observed FEROS spectrumfor the visible range from 4 000 to 7 000 ˚ A . Former spectrum for T41lamb09 (red dashed line) isincluded behind, for comparison purposes. Summary of Lambert-converged CMFGEN models using this new set of stellar parametersare presented on Table 4.4. In this case, we found that the best model corresponds to T40lamb04 , whose spectra for optical and ultraviolet are presented in Fig. 4.15 and Fig. 4.16respectively. Because this new model presents a slightly better fit for some lines such as H α and He II 4684, we will consider this new self-consistent solution given by T eff = 40 kK,log g = 3 . 64 and R ∗ = 18 . R (cid:12) as stellar parameters to be our best fit for ζ -Puppis.6 CHAPTER 4. SOLUTIONS UNDER LAMBERT PROCEDURE N o r m a li se dF l u x λ ( Ang ) Figure 4.16: Fit of T40lamb04 from Table 4.4 (blue solid line) over the observed IUE spectrum forthe ultraviolet range from 1 200 to 1 800 ˚ A . Former spectrum for T41lamb09 (red dashed line) isincluded behind, for comparison purposes. The search for a full self-consistent solution (coupling line-acceleration, hydrodynamics andradiative transfer) for the wind of massive stars has been approached previously by otherstudies (Puls et al., 2000; Kudritzki, 2002). In particular, we emphasise the work done bySander et al. (2017), where another iterative procedure was implemented in order to obtain aself-consistent solution for the wind hydrodynamics of ζ -Puppis under the non-LTE regime,using the radiative transfer code PoWR (Gr¨afener et al., 2002; Hamann & Gr¨afener, 2003). Inthat study, radiative acceleration was also obtained from the output of the radiative transfersolution (PoWR for them, CMFGEN for us), and that acceleration was re-used to obtain anew hydrodynamic profile.Although their study shares a similar philosophy with ours, there are important differences.First, our new velocity profile is calculated by means of the Lambert W -function, whereas their v ( r ) is recalculated by updating the stratification of the wind. The advantage of W -Lambertprocedure yields in the fact that this is a mathematical tool which ensures the existence ofa unique solution at the end, i.e., final converged Lambert-model does not depend on theinitial velocity profile used ( v ∞ and β ). Besides, our procedure considers iterative changes .4. DISCUSSION M ), whereas Sander’s prescription does notallow a relaxation because every parameter is being recalculated. Finally, the methodologyintroduced by us is focused on the line-acceleration g line instead radiative acceleration g rad ,which gives the chance to compare these new results with previous formulations based onm-CAK framework (Chapter 3).Sander et al. (2017) Chapter 3 This workRT code PoWR FASTWIND CMFGENHydro method – HydWind Lambert W -function T eff (kK) 40.7 40 40log g R ∗ /R (cid:12) v ∞ (km s − ) 2 046 2 700 2 310˙ M ( M (cid:12) yr − ) 1 . × − . × − . × − f ∞ (= D − ∞ = f − ) 0.1 0.2 0.15˙ M hom ( M (cid:12) yr − ) 5 . × − . × − . × − Table 4.5: Summary of self-consistent solutions performed by Sander et al. (2017), Gormaz-Matamala et al. (2019) and this present study. These three models assume the same abundances asBouret et al. (2012). Radiative transfer code used to perform the respective synthetic spectra is alsoindicated, together with the code/methodology used to calculate the hydrodynamics. Equivalenciesbetween filling factor and clumping factor f cl (used by FASTWIND) and D ∞ (used by PoWR ) weremade assuming the interclump medium is void, following Sundqvist, & Puls (2018). For those reasons, we select these two previous studies to be compared with the resultsgiven by the Lambert-procedure. Comparison of these three self-consistent solutions for ζ -Puppis is summarised in Table 4.5. As an initial comment, we remark the discrepancybetween the wind parameters derived from Sander et al. (2017) and those derived by us,even when differences on stellar parameters are not so significative. First, for the case of theterminal velocity, their value lies below typical values obtained by spectral fitting: 2 250 kms − by Puls et al. (2006) and 2 300 km s − by Bouret et al. (2012). This is an importantdetail, because our self-consistent wind parameters under Lambert-procedure were calculatedsatisfying the criterium to find an accurate terminal velocity by fitting C IV 1548 line (i.e.,self-consistent mass-loss rate would have been different if we were looking for solutions with v ∞ = 2 046 km s − ). However, since it is derived from Tables 4.3 and 4.4 that mass-loss rate is8 CHAPTER 4. SOLUTIONS UNDER LAMBERT PROCEDURE inversely proportional to the resulting terminal velocity for a self-consistent hydrodynamics, itis clear that a solution with lower v ∞ it would have generated an even higher ˙ M than the onetabulated on Table 4.5. Second, for mass-loss rates we found our self-consistent value doublesthe Sander’s one. Although we employ a less deep clumping factor, our ”homogeneous” ˙ M isa ∼ 50% higher. Nevertheless, because the fit on H α (Fig. 4.15) is more accurate using ourwind parameters than those given by Sander et al. (2017, see their Fig. 9), we consider ourself-consistent solution as more reliable.The comparison with the self-consistent solutions under m-CAK theory performed byGormaz-Matamala et al. (2019) is more extended, because both studies share the same stellarparameters. Therefore, following analysis is focused in the differences between this currentmethodology using Lambert-function with former m-CAK methodology, which had threemain aspects. In Chapter 3 we: • used a simplified ”quasi-NLTE” scenario for the treatment of atomic populations, follow-ing formulations performed by Mazzali & Lucy (1993) and Puls et al. (2000); whereasin this work we solve the proper statistical equilibrium equations when CMFGEN isexecuted. • used a flux field calculated by Tlusty (Lanz & Hubeny, 2003), which uses the plane-parallel approximation; whereas our flux field is calculated within CMFGEN run. • did not consider effects from clumping upon the resulting g line ; whereas, from Fig. 4.12,it is clear that line-acceleration changes when mass-loss rate keeps constant but clumpingfactor is modified. • includes the effects of rotation inside the execution of HydWind , whereas Lambert-procedure have not included any influence of rotation.Some of these factors might explain the differences between wind parameters for ζ -Puppisgiven by Paper I and current study. However, it is important to remark again that currentresults are semi-theoretical: hydrodynamic is self-consistent with line-acceleration and veloc-ity profile is dependent on mass-loss rate, but this latter is set as a free parameter to beconstrained by the spectral fit.Differences on the self-consistent wind parameters between Chapter 3 and this presentchapter are related with the differences on the resulting self-consistent line-acceleration.Fig. 4.17 shows the comparison of both g line , the obtained for the self-consistent Lambert-solution T40lamb04 and the one obtained by the m-CAK self-consistent solution for ζ -Puppispresented on the Chapter 3 by means of: g line = σ e F c k t − α (cid:18) N e, W ( r ) (cid:19) δ , (4.32)with k = 0 . α = 0 . 655 and δ = 0 . 039 (see Table 5 of Gormaz-Matamala et al., 2019).Line-acceleration obtained with m-CAK prescription is higher than g line obtained with the .4. DISCUSSION ( r / R * ) g li n e Figure 4.17: Line-accelerations given by the self-consistent Lambert-procedure (continue blue),compared with the one given by the self-consistent . Lambert-procedure, which could be explain then the resulting mass-loss rate and terminalvelocity presented on the Paper I.Therefore, it is possible to claim that m-CAK theory predicts a higher line-accelerationthan a fully calculated hydrodynamics. The explanation may lie in the three points previouslystated: m-CAK prescription for the calculation of the line-force parameters (and hence g line )uses different atomic populations and different flux field, together with a different treatmentfor the clumping.This could lead us to conclude that differences between these two prescriptions (Chapter 3and this one) lies on the line-acceleration only, but the trend observed on Tables 4.3 and 4.4shows that this is not true. Even if the resulting g line calculated by CMFGEN with parameters˙ M = 5 . × − M (cid:12) yr − and f ∞ = 0 . 2, the resulting self-consistent hydrodynamics wouldhave a terminal velocity far below ∼ − . Hence, the wind parameters predicted bythe m-CAK prescription are not recoverable by the Lambert W -function.The influence of rotational effects is interesting, because it is well known that rotationenhances the values for mass-loss rate on the equator of the star. Hydrodynamics calculatedwith HydWind in Gormaz-Matamala et al. (2019) used a value for the normalised stellarangular velocity of Ω = 0 . 39 in order to reproduce the known value of v sin i = 210 km s − for ζ -Puppis (Bouret et al., 2012), whereas Lambert-procedure do not consider rotational effects Normalised stellar angular velocity is defined as:Ω = v rot v crit , with v crit as defined in Eq. 4.16. CHAPTER 4. SOLUTIONS UNDER LAMBERT PROCEDURE because CMFGEN presents a 1D geometry. According to Venero et al. (2016), mass-lossrates increase their values in a factor of ∼ . − . . 4, which leads us to think thatself-consistent value for ˙ M obtained on Chapter 3 would decrease to ∼ × − M (cid:12) yr − ifΩ = 0, closer to the ˙ M determined by the Lambert-procedure.However, it is important to remind that both m-CAK and Lambert prescriptions, despiteboth obtain self-consistent solutions, work under different philosophies. On Chapter 3, m-CAK prescription calculates its own self-consistent mass-loss rate, which is later tested byspectral analysis in order to check how near or far falls from the real solution. On the otherhand, Lambert-procedure presented on this work calculates a self-consistent solution for thewind hydrodynamics letting the mass-loss rate as a free parameter, which is later constrainedby the spectral fitting . Besides, let us consider the fact that we are working with two differentradiative transfer codes (FASTWIND and CMFGEN) which have discrepancies between them(Massey et al., 2013). Therefore, it is expectable that both methodologies do not find exactlythe same values, but plausible ones deserving to discuss them. In the present Chapter we have presented a methodology to calculate self-consistent hy-drodynamics beyond the m-CAK prescription presented on Gormaz-Matamala et al. (2019,Chapter 3). For that purpose, we have used the CMFGEN radiative transfer code and themathematical Lambert W -function. This function allows us to analytically solve the equa-tion of motion (Eq. 2.2), using the line-acceleration g line given by the solution of the radiativetransfer equation on CMFGEN, to provide us a new velocity profile. This procedure is iter-ated (Lambert-procedure) until the convergence. Hydrodynamic solution given by Lambert-procedure is valid only in the supersonic region of the wind, whereas subsonic region needsto be rescaled in order to obtain a continuous solution.Lambert-procedure has proved to converge into the same hydrodynamic solution, inde-pendent of the initial velocity profile (i.e., terminal velocity and β -value) chosen (Fig. 4.10).On the other hand, because Eq. 2.8 is explicitly independent of density, mass-loss rate is notrecovered by the self-consistent iterations and then it needs to be set as a free parameter.However, dependence on density (and hence dependence on mass-loss rate too because ofequation of continuity Eq. 2.3) is implicitly included in the modified equation of momen-tum by means of the line-acceleration term (Fig. 4.11). Line-acceleration also depends onthe clumping factor (Fig. 4.12), thus the final self-consistent hydrodynamic obtained by theLambert-procedure depends on the initial values for ˙ M and f ∞ introduced to the CMFGENmodels. Therefore, given a specific set of stellar parameters (effective temperature, surfacegravity, stellar radius and abundances), a range of self-consistent hydrodynamics is found fordifferent mass-loss rates and different clumping effects . This does not mean that there are different possible hydrodynamics for a specific set of stellar parameters. .5. CONCLUSIONS FOR LAMBERT-PROCEDURE λ α (indicator for mass-loss rate). Following this recipe, we have found a self-consistent hydrodynamic for the set of stellar parameters presented on Table 4.1, togetherwith another self-consistent hydrodynamic for the set of stellar parameters determined byBouret et al. (2012). Best-fit spectrum found for ζ -Puppis using a self-consistent solution, isobtained for this stellar parameters.Compared with Sander et al. (2017), where all stellar and wind parameters are recalculatedinside the iterative process letting then no free parameters, Lambert-procedure presented inthis work has the advantage of letting some parameters such as the stellar ones and mass-lossrate as free. Hence, these parameters can be tuned independently in order to look for theself-consistent hydrodynamics fitting better the spectral observations. Therefore, followingthe studied correlations found for self-consistent solutions, such as the inverse proportionalitybetween initial mass-loss rate and final terminal velocity (see Table 4.3 and Table 4.4), it ispossible to find the best combination of ( ˙ M , f ∞ ).As a result of the comparison with Chapter 3, we find that both our mass-loss rate andour terminal velocity are lower, despite stellar parameters and the atomic information arethe same for both cases. Even though these differences could be partially explained by thedifferences on the self-consistent line-accelerations found (indicating then that differences onfinal wind parameters would lie on the difference of the methodologies to calculate g line ), amore detailed analysis is needed in order to support this hypothesis.Finally, despite the big computational effort required in an iterative loop involving CMF-GEN, Lambert-procedure has demonstrated to be a confident methodology to find hydrody-namically self-consistent solutions for a stellar wind. Follow-up research would be focused onthe expansion of this prescription for more stars beyond ζ -Puppis. We reproduce a range of solutions for different combinations of ( ˙ M, f ∞ ) because Lambert-procedure is notcapable to calculate a proper density for the wind, but it is clear that only one mass-loss rate and only oneclumping factor are the correct ones for the initial set of stellar parameters. CHAPTER 4. SOLUTIONS UNDER LAMBERT PROCEDURE hapter 5 Fitting Spectra of Massive Starswith First Self-consistentSolutions In this chapter, we will apply the self-consistent methodology, developed, calculated andanalysed using the m-CAK prescription in Chapter 3. To test the accuracy of these spectralfits is a necessary step, because one of the goals of using a new prescription capable of self-consistently describing the hydrodynamics of hot massive stars is to determine both stellarand wind parameters with the help of spectral fitting.For that purpose, we use the observed spectral data from a set of hot massive starsobtained with the Hermes spectrograph . Particularly, we analyse the O type stars HD57682, HD 195592, 9 Sge and HD 192639. This set of stars have been previously used inprevious studies (Grunhut et al., 2017; de Becker et al., 2010; Martins et al., 2015; Bouretet al., 2012), who had already constrained some of the stellar parameters such as effectivetemperature and stellar mass; therefore, they represent a starting point in our search of stellarand wind parameters. However, different from Figures 3.7 to 3.10, where the spectral fit wasfocused on the search for the best clumping filling factor keeping the stellar parameters fixed,here we will proceed modifying all the stellar parameters in order to find the best fit to eachspectrum. In some cases, the new found stellar parameters may lie close to those found inprevious studies, whereas in others they are remarkably different. The consequences andconclusions are presented at the end of this chapter.Most of this work was done during a two months internship at the Royal Observatory ofBelgium in Brussels, under the supervision of Dr. Alex Lobel, as part of the Physics Of High-Efficiency and high-Resolution Mercator Echelle Spectrograph: https://fys.kuleuven.be/ster/instruments/the-hermes-spectrograph CHAPTER 5. SPECTRA WITH SELF-CONSISTENT SOLUTIONS Extremely Massive Stars (POEMS) project . It is important to remark that, this researchwork is still in progress. For future work, many other stellar spectra should be included inour study, for strengthening confidence in our self-consistent m-CAK procedure. Firstly, it is important to emphasise the differences among the parameters that will be con-strained by the spectral fit. Because of the self-consistent procedure, the wind parameters(mass-loss rate and terminal velocity) are depending on the set of initial stellar parametersand cannot be individually modified. • Stellar parameters ( T eff , log g , R ∗ and abundances) are set at the beginning of the m-CAK prescription, and they directly determine the final self-consistent values for mass-loss rate and terminal velocity. Therefore, their modifications are the most complicatedpoint because it requires the execution of a new iterative procedure from the beginning.For this reason, we will call them first-order modifications (FOM) and they are the mostimportant parameters to be fine-tuned during the spectral fitting. • Parameters such as the clumping factor f cl and the turbulence velocity v turb are set atthe beginning of the execution of FASTWIND but are not included in the self-consistentprocedure. We will call them second-order modifications (SOM), and they are modifiedafter FOMs once the self-consistent solution has been achieved. • Final parameters such as macro-turbulence v macro and the rotational velocity v rot areset at the end once the output of FASTWIND is obtained. We will call them third-ordermodifications (TOM). Inside this group we should include the radial velocity v rad , whichfits the core of individual spectral lines.The full picture of the self-consistent procedure plus the execution of the FASTWINDcode is shown in Fig. 5.1. Notice that the mass-loss rate and terminal wind velocity are notlonger independent, because their values are determined by the FOMs at the beginning. Asimilar situation is valid for the line-force parameters represented in Fig. 5.1: final values for k , α and δ are a consequence of the modification of FOMs. Hence, tuning the line-force andwind parameters first requires a fine-tuning of the stellar parameters.It is important to remark, however, the order we presented is related to the nature of theparameters which does not represent an order step for doing the fit. The order step resultsfrom the empirical effects on the lines; for example, the macro-turbulence velocity is onlymodified at the end (TOM), because its value is directly proportional to the width of He Ilines, and it is therefore one of the first parameters to be constrained. http://stelweb.asu.cas.cz/~kraus/POEMS/Secondments.html .1. METHODOLOGY T e↵ , log g , R ⇤ ,abundances HydWind k , ↵ , v ( r ), ⇢ ( r ) M ( t ) = P lines ⌫ D F ⌫ F e ⌘t t self-consistent ˙ M , v f , v turb FASTWIND v macro , v rot convolution FASTWINDsyntheticspectrum Figure 5.1: Scheme of the self-consistent m-CAK procedure combined with FASTWIND, showingthe stage where the different initial parameters are set. For each star, we use parameters derived from previous studies as starting values. Incase of lacking reliable information about a particular star, we use the calibration of stellarparameters by Martins et al. (2005). The most important stellar parameters modified in theself-consistent procedure (effective temperature and stellar radius) are tuned simultaneouslyin order to keep the stellar luminosity constant, according with Stefan-Boltzmann definitionof luminosity: (cid:18) L ini L fin (cid:19) = (cid:18) R ini R fin (cid:19) (cid:18) T ini T fin (cid:19) = 1 . (5.1)Moreover, we mention that inside the item abundances for the FOM, we consider metal-licity, He/H ratio and the individual element abundances. With regards to the last item, weconsider as default the solar abundances of Asplund et al. (2009). Due to the fact that spectraof O type stars are dominated by lines of hydrogen, helium, carbon and nitrogen mainly, it iscurrently not possible to individually fine-tune elements such as the Iron group (iron, cobaltand nickel). However, heavier elements. Elements contributing more to the line-acceleration g line are those having large number of permitted atomic transition lines. Because the irongroup (hereafter FeG) represents ∼ 80% of the total line-force (see Table 3.1) it is impor-tant to investigate if slight changes in the abundances of these elements can affect our finalself-consistent solution. After running many self-consistent models using different abundance6 CHAPTER 5. SPECTRA WITH SELF-CONSISTENT SOLUTIONS values, it is possible to state that the self-consistent mass-loss rate is directly proportional toabundance of iron, cobalt and nickel: ˙ M ∝ Iron-group ∗ Iron-group (cid:12) . (5.2)Hereafter, we will consider the simultaneous fine-tuning of the three elements Fe, Co andNi together, represented as the modification of FeG ∗ /FeG (cid:12) . The group of lines calculated by FASTWIND is presented in Table 5.1.H α β γ δ (cid:15) Lines included in FASTWIND. Blue lines are available for HD 192639 only. As a first stage, we decide to use H and He lines only for the spectral fitting in the opticaland infrared wavelength range, whereas our observational data consist of optical spectra only.Thus, the lines to be fitted correspond to those ones available both in FASTWIND and inthe observational data. In the specific case of our first three stars, we had available only ninespectral lines for the analysis, whereas HD 192639 was studied using fifteen spectral lines inFASTWIND (see Table 5.1).Given the large number of free parameters to fine-tune the spectral fitting (even whendue to the self-consistent procedure we eliminated three of them, see Fig. 5.1), a deeperanalysis is necessary to determine how different elements and different lines are affected bythe modification. After many empirical tests and comparisons letting only one parametervary while the other are held fixed, we have observed some general trends that help us todetermine the best model fit: • Macro-turbulent velocity v macro is fitted to the width of He I lines, whereas v rot mustbe fitted to the shape of the wings of He II λ λ • T eff is fitted to He II λ λ • log g is fitted to the He I lines, particularly He I λ ∗ /FeG (cid:12) also affects the intensity of helium lines, specially the He I lines It would be possible to include lines of silicon and nitrogen, but that requires rebooting the compilationof FASTWIND, together with a big computational effort to run these models. For practical purposes, we willlimit this section to the usage of hydrogen and helium only. .2. RESULTS T eff ), henceboth parameters are fitted simultaneously. • v turb is fitted to He I λ • f cl is fitted to the Balmer lines, specially fitting the peak of H α .The radial velocity, could be obtained from the SIMBAD database , but because we ignoreif the spectrographs were calibrated for v rad we fit it manually to make the core of the linesmatch in wavelength. We present here the spectra of the best fit spectrum for each analysed star. Normalised spectrawere provided by Dr. Alex Lobel. Comments about individual lines are also included. As a first step, we model the wind of the O 9.5 V star HD 57682. Grunhut et al. (2012) havestudied this star, and found line profile variability (LPV) probably due to magnetic fields.We aim to fit spectra self-consistently, in order to check whether or not can reproduce partof the line profiles.Based on Grunhut et al. (2012), we start with the following stellar parameters: T eff = 35kK, log g = 4 . R ∗ = 7 R (cid:12) . In this particular case, the self-consistent solution obtainedfrom this set of stellar parameters quickly leads to a good fit. The corresponding spectrumis presented in Fig. 5.2, whereas the obtained parameters are shown in Table 5.2.Methodology described on Section 5.1 helped us to find good fits to the helium lines,especifically to He II λ λ λ λ λ λ β , but even very high valuefor the clumping factor does not help to reproduce the abnormal shape of the line cores.This is because it is thought that the line cores show the presence of magnetic fields in HD57682 (Grunhut et al., 2017), which are beyond the current capacity of the self-consistent fitprocedure.Compared to Grunhut et al. (2017), we find that our v rot is close to their v sin i (8 kms − ) but their macro-turbulence is too large (65 km s − ). http://simbad.u-strasbg.fr/simbad/ CHAPTER 5. SPECTRA WITH SELF-CONSISTENT SOLUTIONS λ ( Ang ) N o r m a li se dF l u x H α λ ( Ang ) N o r m a li se dF l u x H β λ ( Ang ) N o r m a li se dF l u x H δ λ ( Ang ) N o r m a li se dF l u x He II 4200 λ ( Ang ) N o r m a li se dF l u x He II 4541 λ ( Ang ) N o r m a li se dF l u x He II 4686 λ ( Ang ) N o r m a li se dF l u x He I 4471 λ ( Ang ) N o r m a li se dF l u x He I 4713 λ ( Ang ) N o r m a li se dF l u x He I 4922 Figure 5.2: Best fit for the star HD 57682, using v rad = 10 km s − . We proceed obtaining self-consistent parameters for O9.7Ia (Sota et al., 2011) star HD 195592.The initial mass M ∗ = 30 M (cid:12) , radius R ∗ = 30 M (cid:12) and luminosity L ∗ = 3 . × L (cid:12) weretaken from de Becker et al. (2010), who derived them from Martins et al. (2005). This, andthe following stars, corresponds to late O supergiants, which are near the lower threshold ofvalidity of our self-consistent fit procedure ( T eff ≥ 30 kK and log g ≥ . 4, see Chapter 3).Only for this particular case, we present two fits with different sets of parameters (Fig. 5.3),in order to show the relevance of the modification of the abundances for the iron group alreadymentioned. Parameters used are tabulated in Table 5.3. Fit to helium lines is made by fine-tuning FOM parameters in the same way as the previous case, whereas the hydrogen lines alsorequire to fit clumping factor. The model with solar iron abundance has a higher mass-lossrate, so the clumping factor needed to reproduce the H α core is only f cl = 2 . 0, whereas themodel having 70% of solar abundance fits the observed spectra using a larger clumping factor f cl = 4 . 5. We find that the solar-like model is the best of the two, because it reproduced theemission in He II λ .2. RESULTS T eff (kK) 35.0 35.0log g R ∗ /R (cid:12) M ∗ /M (cid:12) . ∗ / FeG (cid:12) ] 1.0 –( k, α, δ ) (0 . , . , . M ( M (cid:12) yr − ) 1 . × − . × − v ∞ (km s − ) 2 970 1 200 f cl v rot (km s − ) 10 10 v turb (km s − ) 25 – v macro (km s − ) 20 62Table 5.2: Summary of stellar and wind parameters used to fit HD 57682 (Fig. 5.2). Values for v rot and v macro taken from Grunhut et al. (2017). Parameters HD 195592 hd195592fe10 hd195592fe07 de Becker et al. (2010) T eff (kK) 29.5 29.5 28.4log g R ∗ /R (cid:12) ∗ / FeG (cid:12) ] 1.0 0.7 –( k, α, δ ) (0 . , . , . . , . , . M ( M (cid:12) yr − ) 3 . × − . × − – v ∞ (km s − ) 1 060 1 050 – f cl v rot (km s − ) 60 70 60Table 5.3: Summary of stellar and wind parameters used to fit HD 195592. main conclusion we draw is the fact that we can use models with lower/higher iron groupabundance values for constraining the other parameters, especially the mass-loss rate.A remarkable difference compared with our previous fit for HD 57682 is that in this casewe are capable to reproduce the core of H α but not its wings because they are too wide.This problem is also present for the next stars. We therefore concentrate on fitting only the00 CHAPTER 5. SPECTRA WITH SELF-CONSISTENT SOLUTIONS λ ( Ang ) N o r m a li se dF l u x H α λ ( Ang ) N o r m a li se dF l u x H β λ ( Ang ) N o r m a li se dF l u x H δ λ ( Ang ) N o r m a li se dF l u x He II 4200 λ ( Ang ) N o r m a li se dF l u x He II 4541 λ ( Ang ) N o r m a li se dF l u x He II 4686 λ ( Ang ) N o r m a li se dF l u x He I 4471 λ ( Ang ) N o r m a li se dF l u x He I 4922 λ ( Ang ) N o r m a li se dF l u x He I 4713 Figure 5.3: Plot of the models with FeG/FeG (cid:12) = 1 . (cid:12) = 0 . H α core same as for the other lines. The consequences of this disagreement are a matter ofdiscussion. T eff = 34, log g = 3 . 36 and R ∗ = 20 . R (cid:12) , which are also the parameters givenin the VizieR Online Data Catalog (Nebot G´omez-Moran & Oskinova, 2018). Alternatively,Martins et al. (2015) found T eff = 33, log g = 3 . 35. The fitted spectrum and parameters arepresented on Fig. 5.4 and Table 5.4 respectively.Similar to HD 195592, the wings of H α are too wide even when we fit the emission intensity,although the model with the same ˙ M but without clumping fits the wings better but not thecore. Besides, He II λ g has improved with up to decimals. http://vizier.u-strasbg.fr/viz-bin/VizieR?-source=J/A+A/620/A89 .2. RESULTS λ ( Ang ) N o r m a li se dF l u x H α λ ( Ang ) N o r m a li se dF l u x H β λ ( Ang ) N o r m a li se dF l u x H δ λ ( Ang ) N o r m a li se dF l u x He II 4200 λ ( Ang ) N o r m a li se dF l u x He II 4541 λ ( Ang ) N o r m a li se dF l u x He II 4686 λ ( Ang ) N o r m a li se dF l u x He I 4471 λ ( Ang ) N o r m a li se dF l u x He I 4713 λ ( Ang ) N o r m a li se dF l u x He I 4922 Figure 5.4: Best fit for the star 9 Sge, using v rad = 30 km s − . Here we present our analysis using self-consistent solution for the supergiant O 7.5 I, HD192639, observed with Hermes in 2013. The initial stellar parameters ( T eff = 32 kK andlog g = 3 . 36) were taken from standard calibrations done by Martins et al. (2005) using theirtheoretical values for effective temperature. The initial rotational velocity is set to v rot = 80km s − .The stellar and wind parameters obtained for the fit to HD 192639 are given in Table 5.5,whereas the fit is shown in Fig. 5.5. This time, we have included all hydrogen and helium linesavailable on FASTWIND. The effective temperature is higher than the initial value, whereasthe surface gravity is ∼ 10% smaller. Concerning the spectral fitting, only for this star itwas possible to get a good fit to both the wings and the Hα emission peak. We suppose itis because HD 192639 is a standard star (Bouret et al., 2012; Martins et al., 2015), thereforeits spectrum does not show signs of any unusual behaviour such as magnetic fields or fastrotational velocity.02 CHAPTER 5. SPECTRA WITH SELF-CONSISTENT SOLUTIONS Parameters 9 SgeThis work Nebot G´omez-Moran & Oskinova (2018) T eff (kK) 34.5 34.0log g R ∗ /R (cid:12) M ∗ /M (cid:12) . ∗ / FeG (cid:12) ] 1.0 –( k, α, δ ) (0 . , . , . M ( M (cid:12) yr − ) 4 . × − . × − v ∞ (km s − ) 1 700 – f cl v rot (km s − ) 80 60 v turb (km s − ) 20 – v macro (km s − ) 30 46Table 5.4: Summary of stellar and wind parameters used to fit 9 Sge (Fig. 5.4). λ ( Ang ) N o r m a li se dF l u x H α λ ( Ang ) N o r m a li se dF l u x H β λ ( Ang ) N o r m a li se dF l u x H γ λ ( Ang ) N o r m a li se dF l u x H δ λ ( Ang ) N o r m a li se dF l u x H ϵ λ ( Ang ) N o r m a li se dF l u x He II 4200 λ ( Ang ) N o r m a li se dF l u x He II 4541 λ ( Ang ) N o r m a li se dF l u x He II 4686 λ ( Ang ) N o r m a li se dF l u x He II 6527 λ ( Ang ) N o r m a li se dF l u x He II 6683 λ ( Ang ) N o r m a li se dF l u x He I 4387 λ ( Ang ) N o r m a li se dF l u x He I 4471 λ ( Ang ) N o r m a li se dF l u x He I 4713 λ ( Ang ) N o r m a li se dF l u x He I 4922 λ ( Ang ) N o r m a li se dF l u x He I 6678 Figure 5.5: Best fit for the star HD 192639, using v rad = − 10 km s − . .3. SUMMARY AND FUTURE WORK T eff (kK) 33.5 33.5log g R ∗ /R (cid:12) M ∗ /M (cid:12) . ∗ / FeG (cid:12) ] 1.0 –( k, α, δ ) (0 . , . , . M ( M (cid:12) yr − ) 4 . × − . × − v ∞ (km s − ) 1 380 1 900 f cl v rot (km s − ) 100 90 v micro (km s − ) 25 – v macro (km s − ) 30 43Table 5.5: Summary of stellar and wind parameters used to fit HD 192639 (Fig. 5.5). We have presented spectral fits combined with our self-consistent m-CAK procedure (Gormaz-Matamala et al., 2019) calculating synthetic spectra with the NLTE radiative transfer codeFASTWIND. Different from previous studies, the parameters we obtain for the wind (mass-loss rate and terminal velocity) are self-consistent with the line-acceleration and the hydro-dynamics of the wind. For that reason, in almost all the cases they differ for the stellarparameters determined in previous studies. Our work presents new best parameters for char-acterising massive stars.The methodology presented here has demonstrated be able to accurately fit the heliumlines, with the exception of He II λ α ): H α yields a good fit on the line cores but not the wings. It is well known that H α is very importantfor the determination of mass-loss rates but also its shape can reveal other phenomena suchas accretion disks or magnetic fields (which are beyond the scope of the FASTWIND andself-consistent procedure presented here).Another important new result, is that we have been able to fit the spectra with stellarparameters below the threshold stated by Gormaz-Matamala et al. (2019). This is an indicatorthe the lower limits mentioned there should be revisited.As for future work, we mention the necessity of including all the hydrogen and helium04 CHAPTER 5. SPECTRA WITH SELF-CONSISTENT SOLUTIONS lines available on FASTWIND for the first fitted stars. Besides, we have not included yetmodifications of CNO abundances for 9 Sge and HD 192639, which are detailed in Martinset al. (2015). Finally, an appropriate prescription for constraining stellar parameters withreliable error bars is necessary in order to use the new parameters in subsequent studies suchas the determination of stellar masses of massive binary systems. hapter 6 Evolution of Massive Stars withself-consistent HydrodynamicModels The main goal of the present thesis has is development of a methodology capable to coupleboth, the calculation of line-acceleration and the hydrodynamics of the wind (i.e., a self-consistent solution). We have performed a self-consistent prescription based on the m-CAKtheory, which has demonstrated to give confident values for wind parameters (mass-loss rateand terminal velocity, see Chapter 3). The contrast between these m-CAK solutions, undera quasi-NLTE scenario for their atomic populations, and a self-consistent solution under afull NLTE scenario by CMFGEN was presented and discussed on Chapter 4. Even thoughsome differences, self-consistent m-CAK prescription has been capable to provide us accuratespectral fits for massive stars, as it is shown in Chapter 5. Thus, self-consistent mass-lossrates can be implemented for future studies, such as the evolution of massive stars.Mass-loss rate, as we already briefly mentioned on Chapter 1, plays a key role in theevolution of massive stars. The pioneer work from Meynet et al. (1994) summarised us thateven changes on ˙ M by a factor of two may dramatically affect the fate of a star. More recentstudies (Ekstr¨om et al., 2012; Groh et al., 2019) have confirmed the same trend. Particularly,recent studies as the mentioned above calculated their evolutionary tracks for massive starsusing mass-loss values set by Vink’s formula (Vink et al., 2001) which, as we pointed out onChapter 3, is not self-consistent. The study of the evolution of massive stars that includeself-consistent solutions for their stellar winds is still not fully performed.For that reason, on this chapter we present new evolutionary tracks for a set of massivestars, assuming mass-loss rate given by the self-consistent procedure instead of Vink’s formula.We use the Geneva evolutive code Genec (Maeder, 1983, 1987; Maeder, & Meynet, 1987) toperform the evolutionary tracks, editing the code to include our prescription. The foundations10506 CHAPTER 6. EVOLUTION OF SELF-CONSISTENT SOLUTIONS of the work presented in this Chapter was performed during a four months internship atTrinity College Dublin in Ireland, under the supervision of Professor Dr. Jose Groh. Followingthe outline presented on Section 3.3, this analysis will be done for a case with solar metallicity( Z (cid:12) = 0 . Z = 0 . Previously, when we presented the results for wind parameters under self-consistent m-CAKprescription, we introduced Eq. 3.13 and Eq. 3.14 as linear relationships to calculate self-consistent mass-loss rate from stellar parameters, in an analogous way to the formulae pre-sented by Vink et al. (2001). But, even when the coefficient of determination was extremelyconfident, these expressions were limited to the range of validity given in Gormaz-Matamalaet al. (2019), namely T eff ≥ 32 kK and log g ≥ . 4. Spectral fits shown in Chapter 5 havesuggested that these lower thresholds could be decreased in order to be valid for more stars,specially O supergiants. For that reason, we proceed to calculate more self-consistent solu-tions for a new set of standard stellar parameters in order to improve these relationships,complementing then the results presented on Table 3.3.For this time however, standard parameters will be taken for those obtained after running Genec . We select the evolutive tracks for non-rotating stars with 25, 40, 70 and 120 solarmasses, picking then four representative points for each track (only three for the case of 25 M (cid:12) ). The points were selected in order to be the most representative on temperature andalmost equitable, as it is shown in Fig. 6.1. The search of new line-force parameters is notlimited only for the variations on effective temperature, surface gravity and stellar radius,but also for the abundances. Even when metallicity can be considered as constant throughthe entire evolutionary track, changes in the He to H ratio or in the individual abundance ofmetal elements affects the resulting line-acceleration and therefore affecting the final mass-loss rate. These effects produced by abundances are also studied, in order to incorporatethem at the final expressions for ˙ M . Given Genec does not provide us surface abundancesfor important metals such as the iron-group (in order to include a mathematical expressionanalogous to the relationship shown in Eq. 5.2), our analysis of individual metal elements isconstrained to CNO elements only.The new line-force parameters, with their respective new self-consistent mass-loss rate,are presented in Table 6.1. We have also included columns to display the stellar mass and theluminosity at each stage. It is important to notice that changes on surface abundances for Genec provides us the change of abundances both on the core of the star and its surface. However, forthe calculation of line-acceleration we are interested only on the modification of abundances in the surface ofthe star. .1. METHODOLOGY M sun M sun M sun M sun T eff [ kK ] l og L * / L s un Figure 6.1: Evolutionary tracks for stars with 120, 70, 40 and 25 M (cid:12) at solar metallicity ( Z = 0 . CNO elements and the He to H ratio is just presented for the case of 120 M (cid:12) , because is theonly one evolutive track that exhibits this behaviour . In these cases, a self-consistent solutionwith the abundances unchanged is presented first, followed by self-consistent solutions afterthe single modification of He to H ratio and the single CNO elements, to finish showing thelast self-consistent solution with all the abundance modifications included.Same procedure must be done for the case of metallicity Z = 0 . This is a direct consequence of being working with non-rotating models: reactions in the core abruptlymodify the abundance structure at some point for only extreme high mass cases such as 120 M (cid:12) . Evolutivetracks that consider rotation exhibit a gradual modification of surface abundances even for stars with massesbelow 25 M (cid:12) (Ekstr¨om et al., 2012). CHAPTER 6. EVOLUTION OF SELF-CONSISTENT SOLUTIONS T eff log g R ∗ M ∗ log L ∗ k α δ v ∞ ˙ M SC ˙ M SC / ˙ M Vink [K] [ R (cid:12) ] [ M (cid:12) ] [ L (cid:12) ] [km s − ] [ M (cid:12) yr − ]54 000 4.15 15.0 116.0 6.24 0.126 0.660 0.014 5 148 8 . × − . × − . × − [He/H] = 0 . . × − [C/C (cid:12) ] = 0 . . × − [O/O (cid:12) ] = 0 . . × − [He/H] = 0 . 11, mod. in CNO abun. . × − . × − [He/H] = 0 . . × − [C/C (cid:12) ] = 0 . . × − [N/N (cid:12) ] = 6 . . × − [O/O (cid:12) ] = 0 . . × − [He/H] = 0 . 14, mod. in CNO abun. . × − . × − . × − . × − . × − . × − . × − . × − . × − . × − . × − . × − Table 6.1: Analogous to Table 3.3, self-consistent line-force parameters ( k, α, δ ) for adopted standardstellar parameters, together with the resulting terminal velocities and mass-loss rates ( ˙ M SC ). Ratiosbetween self-consistent mass-loss rates and Vink’s recipe values (re-scaled to match metallicity fromAsplund et al., 2009) using v ∞ /v esc = 2 . .2. EVOLUTIONARY TRACKS AT SOLAR METALLICITY M sun M sun M sun M sun T eff [ kK ] l og L * / L s un Figure 6.2: Evolutionary tracks for stars with 120, 70, 40 and 25 M (cid:12) at low metallicity ( Z = 0 . Combining solutions presented on both Table 3.3 and Table 6.1, we perform a new fit toobtain the self-consistent formula for mass-loss rate as function of the stellar parameters.log ˙ M Z = Z (cid:12) = 42 . × (cid:20) log (cid:18) T eff (cid:19)(cid:21) − − . × (cid:20) log (cid:18) T eff (cid:19)(cid:21) − − × (log g ) − + 178 . × (log g ) − − . × ( R ∗ /R (cid:12) ) − − . − . × [log(He/H) − log(He/H) i ] − . × log(C/C i ) − . × log(N/N i ) . (6.1)Unlike Eq. 3.13, the dependence of Eq. 6.1 on the stellar parameters is not linear. More-over, we included the already mentioned dependence on abundance modifications. This time,we consider this formula (Eq. 6.1) valid while log g ≥ . 2, whereas for surface gravities belowthat value, the mass-loss rate will be determined by Vink’s formula (Vink et al., 2001) .Resulting new evolutionary tracks are presented on Fig. 6.3. Our first comment, it is It is necessary to remind that Genec uses more than one recipe to calculate the mass-loss rate, dependingon the evolutive stage of the star (Ekstr¨om et al., 2012). For the particular case of stars with M ∗ ≥ M (cid:12) (which is the focus of this chapter), Vink’s formula is used from ZAMS to the end of the Main Sequence,whereas for the Wolf-Rayet stage it is used the recipe given by Gr¨afener, & Hamann (2008). CHAPTER 6. EVOLUTION OF SELF-CONSISTENT SOLUTIONS T eff log g R ∗ M ∗ log L ∗ k α δ v ∞ ˙ M SC ˙ M SC / ˙ M Vink [K] [ R (cid:12) ] [ M (cid:12) ] [ L (cid:12) ] [km s − ] [ M (cid:12) yr − ]58 000 4.3 12.8 120.0 6.23 0.114 0.537 0.010 3 790 1 . × − . × − [N/N (cid:12) ] = 0 . . × − . × − [N/N (cid:12) ] = 2 . . × − . × − [N/N (cid:12) ] = 2 . . × − . × − . × − . × − . × − . × − . × − . × − . × − . × − . × − . × − . × − Table 6.2: Summary of standard stars extracted from the initial evolutionary paths for low metal-licity ( Z = 0 . (cid:12) = 0 . clearly seen that the cases using self-consistent mass-loss rates exhibit slightly more luminouspaths. Additionally, we show the temporal evolution of the surface gravity on Fig. 6.4. Foreach track, we have highlighted three points, each one representing a specific stage: • Zero Age Main Sequence: the beginning of our evolutive tracks. • Abundances X He = X H = 0 . • Surface gravity log g = 3 . 3: the ending point of our prescription, entering then a lineartransition until return to Vink’s formula for log g = 3 . M issmaller that the ˙ M given by Vink’s recipe, the resulting evolutionary paths remain with alower mass-loss rate until the end of the self-consistent prescription and thus, the stellar massfunction for stars decreases slower than before . Result presented on left-panel of Fig. 6.5 seems to be in contradiction with the last column of Table 6.1,where it is shown that ˙ M SC > ˙ M Vink . However, ratios on the table represent a comparison between self- .2. EVOLUTIONARY TRACKS AT SOLAR METALLICITY M sun M sun M sun M sun T eff [ kK ] l og L * / L s un Figure 6.3: Comparison of evolutionary tracks for stars with 120, 70, 40 and 25 M (cid:12) at solarmetallicity ( Z = 0 . M SC (solid lines) and using Vink’sformula ˙ M Vink (dashed lines). M sun M sun M sun M sun [ Myr ] l ogg Figure 6.4: Evolution of log g of evolutionary tracks for stars with 120, 70, 40 and 25 M (cid:12) atsolar metallicity ( Z = 0 . M SC (solid lines) and using Vink’sformula ˙ M Vink (dashed lines). consistent and Vink’s mass-loss rates for the old evolutionary track using Vink’s recipe , whereas Fig. 6.5represents the evolution of the mass-loss rate for each different evolutionary tracks, so the ˙ M SC plotted are CHAPTER 6. EVOLUTION OF SELF-CONSISTENT SOLUTIONS M sun M sun M sun M sun - - - - - [ Myr ] l og ( M ) M sun M sun M sun M sun [ Myr ] M s t a r [ M s un ] Figure 6.5: Evolution of ˙ M (left panel) and M ∗ (right panel) of evolutionary tracks for stars with120, 70, 40 and 25 M (cid:12) at solar metallicity ( Z = 0 . M SC (solid lines) and using Vink’s formula ˙ M Vink (dashed lines). T eff log g R ∗ M ∗ log L ∗ k α δ v ∞ ˙ M SC [K] [ R (cid:12) ] [ M (cid:12) ] [ L (cid:12) ] [km s − ] [ M (cid:12) yr − ]45 000 4.25 7.8 39.8 5.35 0.191 0.574 0.023 3 540 4 . × − 41 000 3.99 10.5 39.6 5.46 0.148 0.605 0.024 3 300 7 . × − 30 000 3.30 22.9 40.0 5.35 0.061 0.671 0.142 1 750 8 . × − Table 6.3: Stellar and line-force parameters for the set of black points highlighted on Fig. 6.4 forthe case with 40 M (cid:12) . The most noticeable new result, is that the usage of self-consistent mass-loss rates producemore luminous paths, increased approximately a ∼ g = 3 . M SC . Even when these differences are on the edge of thesignificance, their consequences still need to be understood. In order to have an idea about how stellar spectra are changing through the evolutionarytracks, we decide to run self-consistent solutions for each one of the three points previouslymentioned. Motivation is, to perform an analysis on the evolution of spectra similar to theworks previously developed (Groh et al., 2014; Liermann, 2015), but using the self-consistent valid for the new path and then they do not belong to the parameters tabulated on Table 6.1. .2. EVOLUTIONARY TRACKS AT SOLAR METALLICITY λ ( Ang ) N o r m a li se dF l u x H α λ ( Ang ) N o r m a li se dF l u x H β λ ( Ang ) N o r m a li se dF l u x H δ λ ( Ang ) N o r m a li se dF l u x He II 4200 λ ( Ang ) N o r m a li se dF l u x He II 4541 λ ( Ang ) N o r m a li se dF l u x He II 4686 λ ( Ang ) N o r m a li se dF l u x He I 4026 λ ( Ang ) N o r m a li se dF l u x He I 4471 λ ( Ang ) N o r m a li se dF l u x He I 4922 Figure 6.6: Synthetic spectra created by FASTWIND for the three highlighted stars tabulated onTable 6.3: T eff = 45 kK (red), T eff = 41 kK (blue) and T eff = 30 kK (green). All models were executedwithout convolution by rotation and with a clumping factor f cl = 5 . hydrodynamics and FASTWIND. We will present here the evolution of the track for 40 M (cid:12) .Self-consistent solutions for the three stages highlighted on the 40 M (cid:12) evolutionary trackare shown in Table 6.3, whereas Fig. 6.6 presents the synthetic spectra. It is interesting toobserve that, for the first two stages there is no significant differences on the profiles of thehydrogen lines; the most remarkable difference is the absorption intensity of the He I lines,presumably because the lower ionisation at lower temperature (He II/He I ratio becomessmaller as effective temperature decreases). This change on the ionisation is visible for thethird stage too, specially with the weakness of the He II lines, but besides we observe emissioncomponents on H α . This emission is an effect produced by the reduction of the surface gravitymore than an effect produced by the clumping: from the fits previously performed for starsHD 164794 (Fig. 3.10, log g = 3 . 92) and HD 163758 (Fig. 3.9, log g = 3 . α even for homogeneous models,whereas increase on clumping does not produce emission for large values for log g . At thesame time, this emission component for H α may be produced by the change on the line-force parameters, specially the increase on the values for α and δ (see Fig. 2 from Araya &14 CHAPTER 6. EVOLUTION OF SELF-CONSISTENT SOLUTIONS M sun M sun M sun M sun T eff [ kK ] l og L * / L s un Figure 6.7: Comparison of evolutionary tracks for stars with 120, 70, 40 and 25 M (cid:12) at solarmetallicity ( Z = 0 . M SC (solid lines) and using Vink’sformula ˙ M Vink (dashed lines). Cur´e, 2017). However, a more complete analysis considering the evolutionary tracks for allthe studied stellar masses and with a more extended analysis at different clumping values isrequired in order to confirm or discard this preliminary trends. Repeating the procedure shown in previous section, we combine the self-consistent solutionsfrom Table 3.3 with those presented on Table 6.2 in order to create a specific formula for thecase Z/Z (cid:12) = 0 . M Z =0 . Z (cid:12) = − . × (cid:20) log (cid:18) T eff (cid:19)(cid:21) − − . × (cid:20) log (cid:18) T eff (cid:19)(cid:21) − + 33 . × (log g ) − − . × (log g ) − − . × ( R ∗ /R (cid:12) ) − − . − . × log(N/N i ) . (6.2)The new evolutionary tracks are presented on Fig. 6.7, whereas the evolution of surfacegravities is shown in Fig. 6.8. In this low-metallicity case, even though new paths are slightlymore luminous, differences with Vink’s tracks is negligible. .3. EVOLUTIONARY TRACKS AT LOW METALLICITY M sun M sun M sun M sun [ Myr ] l ogg Figure 6.8: Evolution of log g of evolutionary tracks for stars with 120, 70, 40 and 25 M (cid:12) atsolar metallicity ( Z = 0 . M SC (solid lines) and using Vink’sformula ˙ M Vink (dashed lines). M sun M sun M sun M sun - - - - - [ Myr ] l og ( M ) M sun M sun M sun M sun [ Myr ] M s t a r [ M s un ] Figure 6.9: Evolution of ˙ M (left panel) and M ∗ (right panel) of evolutionary tracks for stars with120, 70, 40 and 25 M (cid:12) at low metallicity ( Z = 0 . M SC (solidlines) and using Vink’s formula ˙ M Vink (dashed lines). The behaviour of the masses however, shows that mass-loss rates are comparatively farbelow previous values, in contrast with the case at solar metallicity. In spite of this, theeffects on the reduction of stellar mass is lower because the value of the mass-loss rate at lowmetallicity is smaller than the solar metallicity mass-loss rate.16 CHAPTER 6. EVOLUTION OF SELF-CONSISTENT SOLUTIONS Given the new self-consistent solutions for stellar winds on massive stars performed on thisThesis (Chapter 3, Gormaz-Matamala et al., 2019), we have started to use the resulting mass-loss rate to calculate new evolutionary tracks as an alternative to those performed underVink’s recipe. Implementing this in Genec code, these new paths over the Hertzsprung-Russell diagram for a set of massive stars at solar and low metallicities are shown in Fig. 6.3and Fig. 6.7. The analysis of these results and their consequences are still a work in progress.By the obtention of new self-consistent solutions for a new set of stellar parameters, formerquick formulae to calculate ˙ M SC (Gormaz-Matamala et al., 2019, Eq. 3.13 and Eq. 3.14) havebeen improved into Eq. 6.1 and Eq. 6.2, where it has been included the dependence on themodification of CNO abundances and the He to H ratio. In order to ensure our confident onthese new formulae, it would be necessary to run more self-consistent solutions to improvethe statistics.As we mention at the beginning of this chapter, evolutionary tracks for massive stars usingself-consistent hydrodynamics is still a not fully performed topic, and therefore it has a lotof potential to research. However, as a first stage we will concentrate on the non-rotatingcase only. Self-consistent hydrodynamics under m-CAK prescription are implemented toinclude the effects of rotation (Cur´e, 2004; Araya et al., 2018), but a good connection withthe correction factor over mass-loss rates described by Maeder, & Meynet (2000) is neededin order to perform a coherent methodology, a work that is currently beyond the resultspresented here. hapter 7 Summary and Conclusions The important role that the massive stars play on the field of Stellar Astrophysics, makes theminteresting objects to study. In particular, the main link that massive stars have with othertopics of Astronomy such as Galactic Astrophysics or Cosmology, is the chemical enrichmentand energy output produced by their powerful stellar winds.The main feature of the massive stars is their strong mass-loss rates produced by line-driven stellar winds, which largely affects their evolution and therefore their contribution tothe galactic chemical enrichment. For that reason it is necessary to properly understand thephysics behind the stellar wind in order to better constrain values for mass-loss rates. Formassive stars, line driven winds theory has provided a quite complete theoretical frameworkcapable to predict values for mass-loss rate in the order of the observed values, with smalluncertainties (Puls et al., 2008). In the frame of line-driven theory, m-CAK prescription(Castor et al., 1975; Pauldrach et al., 1986) have also demonstrated to be versatile andaccurate enough to describe the hydrodynamics of the wind and also predict mass-loss rates.Under the framework of this line-driven wind theory, through this thesis we have per-formed prescriptions to calculate the line-acceleration in an iterative way combined with thehydrodynamics of the wind, in other words, a self-consistent solution for the stellar windon massive stars. Solutions presented here not only satisfy the concordance between theline-acceleration and the hydrodynamics, but also have been calculated using new numericalprocedures, new mathematical tools and new atomic information not considered by previousstudies, having been self-consistent or not. Besides, from the provided solutions is easy to cal-culate the respective synthetic spectra in order to evaluate their stellar and wind parameterswith observations.The most important part of the work along this thesis has been performed under the m-CAK prescription (see Chapter 3). We have developed a methodology capable to parametrisethe line-acceleration calculating the force multiplier M ( t ) and their line-force parameters k , α and δ , using as input an updated atomic database and hydrodynamic profiles fromthe HydWind code (Cur´e, 2004), iterating until a convergence criterion is attained. This11718 CHAPTER 7. SUMMARY & CONCLUSIONS prescription has the big advantage to converge in a short timescale ( ∼ 10 minutes), but itincludes some assumptions and approximations such as an approximated NLTE treatmentfor atomic populations and the assumption that line-force parameters are constant throughthe wind. Moreover, the consideration of density inhomogeneities on the wind is not formallyincluded in the iterative calculation of the self-consistent solution and it is only introduced forthe execution of the synthetic spectra with FASTWIND. Even with these disadvantages, m-CAK prescription has proven to be successful in the calculation of a line-acceleration capableto lead into accurate mass-loss rates, and the synthetic spectra performed by FASTWINDquickly have delivered a good fit in the neighbourhood of the ”real” solution compared withobserved spectra of massive stars.It has been found that self-consistent line-force parameters k , α and δ strongly depends oninitial stellar parameters determined for the model. Dependence in temperature is known andexpected because of the force multiplier (Eq. 2.46) is implicitly dependent in temperature.However, k , α and δ have demonstrated to have a dependence on surface gravity too, despitethe fact that there is no term related with mass for M ( t ). This new dependence comes fromthe iterative procedure, because this time the optical depth t is properly calculated fromthe hydrodynamics via HydWind (instead using a standard grid as Abbott, 1982, or a β -law). More surprisingly, dependence on log g seems to be deeper than dependence on T eff forline-force multipliers, specially for surface gravities in the range of log g ∼ . − . g ≥ . ∼ v ∞ and the line-force parameter α , where we found that self-consistent solutions satisfy Eq. 3.15for low metallicities only. Solar metallicities presents a more spread distribution, where thelinear relationship is steeper and depends again on the initial surface gravity. The big influ-ence from the surface gravity (or implicitly influence from stellar mass) is a matter for futurediscussion and analysis, specially at the moment to evaluate evolutionary tracks at differentinitial masses.Concerning to mass-loss rate (the most interesting case for us), self-consistent valuesare higher than those observationally determined, but still below the predicted values by theVink’s formula (Vink et al., 2001). These discrepancies between the theoretical self-consistentmass-loss rate ˙ M SC and the observed one ˙ M obs might partially be explained by the usage of β -law on the spectral fitting, but also the consideration of inhomogeneities in the wind (i.e.,clumping) plays a role: self-consistent values for ζ -Puppis (HD 66811) and HD 163758 arein the order of 2 times higher than the clumped values found by Bouret et al. (2012), but˙ M SC for the unclumped stars analysed by Markova et al. (2018) differs only in the order of19 ∼ − 30% with their observed ˙ M obs . However, even for the clumped cases it has been foundthat self-consistent mass-loss rates reproduce accurate synthetic spectra with FASTWINDadjusting clumping factor (see Figures from 3.7 to 3.10). Therefore, it is possible to concludethat self-consistent solutions under the m-CAK procedure provide accurate mass-loss ratetheoretical values, that can be used for future studies. Indeed, m-CAK prescription has beenused as a basis to calculate new stellar and wind parameters for a set of spectra of massivestars by means of the spectral fitting with FASTWIND.Simultaneously, we have analysed the case of self-consistent solutions beyond the m-CAKprescription; that means, with a line-acceleration calculated without the assumptions im-plemented for the previous case such as Sobolev approximation or quasi-NLTE atomic pop-ulations. Besides, flux field is directly determined by solving the NLTE radiative transferequation. This line-acceleration is obtained from the output of the radiative transfer codeCMFGEN, whereas the new hydrodynamics are calculated with the help of the mathematicaltool called Lambert W -function. Iterative combination of the execution of CMFGEN and theLambert W -Function to recalculate v ( r ) is called Lambert-procedure. This full NLTE pre-scription also provides us an accurate synthetic spectra to fit our standard star ( ζ -Puppis), butthe self-consistent mass-loss rate obtained this time is ∼ 40% lower than under the m-CAKprescription, both cases using the same stellar parameters (see Table 4.5). This discrepancymay be attributed to the differences on both methods to determine the self-consistent line-acceleration, such as the treatment for atomic populations, the inclusion of clumping, theflux field used or the rotational effects. Concerning this last point, we have discussed the factthat differences on mass-loss rate would be only around a ∼ 25% if rotation had not beenincluded for the m-CAK self-consistent solution. However, it is important to remark also thatwe are not expecting to reproduce the same self-consistent wind parameters for m-CAK andLambert-procedure because both prescriptions are using different radiative transfer codes tocalculate their synthetic spectra, and the differences between them were already outlined byMassey et al. (2013). The main conclusion is that both self-consistent methodologies predicthigher values for mass-loss rates but with a lower clumping factor than previously calculatedby former studies. Indeed, both self-consistent values for f cl (see Table 4.5) are closer to thevalue calculated by Puls et al. (2006, f cl = 5 . 0) than to more recent extremely high valuesfrom Bouret et al. (2012, f cl = 20) or Sander et al. (2017, f cl = 10). Nevertheless, it isnecessary a more extended analysis in order to understand the physics about the implicationsof these differences on clumping.In spite of the differences on the m-CAK and Lambert-procedures, the first one has thegreat advantage of the time and CPU saving, so we employ it to obtain theoretical valueson wind parameters for subsequent studies. Besides the spectral fitting done over a set ofHERMES spectra, we use the self-consistent mass-loss rates to perform new evolutionarytracks using the Genec code. For main-sequence stars (log g ≥ . M SC arebelow the classical ones obtained from the Vink’s recipe, originating then evolutive tracks20 CHAPTER 7. SUMMARY & CONCLUSIONS with a higher retention of the stellar mass and more luminous. It is interesting to notice that,even when ˙ M SC > ˙ M Vink for log g ≤ . 5, the opposite initial trend produced at the beginningkeeps unchanged over the new resulting evolutionary track. Ergo, main difference on this newevolutionary tracks with the classical ones lies in the fact that self-consistent mass-loss ratesare below the Vink’s values for higher surface gravities. Both studies, even when their scopesare part of the forthcoming work beyond the thesis presented here, are examples of a novelresearch in the field of massive stars derived from the self-consistent solutions for the stellarwinds presented in this thesis. We aim to exploit their potential in the ongoing years. ibliography Abbott, D. C., 1982 ApJ, 259, 282 AAbbott, D. C., & Lucy, L. B. 1985, ApJ, 288, 679Araya, I., Cur´e, M., & Cidale, L. S. 2014, ApJ, 795, 81Araya, I., & Cur´e, M. 2017, The Lives and Death-throes of Massive Stars, 383Araya, I., Cur´e, M., ud-Doula, A., et al. 2018, MNRAS, 477, 755Asplund, M., Grevesse, N., Sauval, A. J., & Scott, P. 2009, ARAA, 47, 481Bouret, J.-C., Lanz, T., & Hillier, D. J. 2005, A&A, 438, 301Bouret, J.-C., Hillier, D. J., Lanz, T., & Fullerton, A. W. 2012, A&A, 544, A67Cardona, O., Mart´ınez-Arroyo, M., & L´opez-Castillo, M. A. 2010, ApJ, 711, 239Carroll, B. W., & Ostlie, D. A. 1996, Institute for Mathematics and Its Applications,Castor, J. I., Abbott, D. C. & Klein, R. I., 1975 ApJ, 195, 157C (CAK)Corless, R. M., Gonnet, G. H., Hare, D. E. G., & Jeffrey, D. J. 1993, The Maple Technical Newsletter,9, 12Corless, R. M., Gonnet, G. H., Hare, D. E. G., Jeffrey, D. J., & Knuth, D. E. 1996, Advances inComputational Mathematics, 5, 329Crowther, P. A., 2007, ARAA, 45, 177.Cur´e, M., 2004, ApJ, 614, 929Cur´e, M., & Rial, D. F. 2004, A&A, 428, 545Cur´e, M., & Rial, D. F. 2007, Astronomische Nachrichten, 328, 513Cur´e, Cidale, L., & Granada, A. 2011, ApJ, 737, 18Cur´e, M., Cidale, L., & Rial, D. F. 2012, ApJ, 757, 142de Becker, M., Linder, N., & Rauw, G. 2010, NA, 15, 76Ekstr¨om, S., Georgy, C., Eggenberger, P., et al. 2012, A&A, 537, A146 BIBLIOGRAPHY Fierro-Santill´an, C. R., Zsarg´o, J., Klapp, J., et al. 2018, ApJS, 236, 38Friend, D. B., & Abbott, D. C. 1986, ApJ, 311, 701Gayley, K. G. 1995, ApJ, 454, 410Georgy, C., Ekstr¨om, S., Meynet, G., et al. 2012, A&A, 542, A29Georgy, C., Ekstr¨om, S., Eggenberger, P., et al. 2013, A&A 558, A103Nebot G´omez-Moran, A., & Oskinova, L. M. 2018, VizieR Online Data Catalog, J/A+A/620/A89Gormaz-Matamala, A. C., Herv´e, A., Chen´e, A.-N., et al. 2015, New Windows on Massive Stars, 100Gormaz-Matamala, A. C., Cur´e, M., Cidale, L. S., & Venero, R. O. J. 2019, ApJ, 873, 131Gray, R. O., & Corbally, C. 2009, Stellar Spectral Classification by Richard O. Gray and ChristopherJ. Corbally. Princeton University PressGr¨afener, G., Koesterke, L., & Hamann, W.-R. 2002, A&A, 387, 244Gr¨afener, G., & Hamann, W.-R. 2008, A&A, 482, 945Groh, J. H., Meynet, G., Ekstr¨om, S., et al. 2014, A&A, 564, A30Groh, J. H., Ekstr¨om, S., Georgy, C., et al. 2019, A&A, 627, A24Grunhut, J. H., Wade, G. A., Sundqvist, J. O., et al. 2012, MNRAS, 426, 2208Grunhut, J. H., Wade, G. A., Neiner, C., et al. 2017, MNRAS, 465, 2432Hamann, W.-R., & Gr¨afener, G. 2003, A&A, 410, 993Hillier, D. J. 1990, A&A, 231, 111Hillier, D. J. 1990, A&A, 231, 116Hillier, D. J., & Miller, D. L. 1998, ApJ, 496, 407Hillier, D. J., & Lanz, T. 2001, in Astronomical Society of the Pacific Conference Series, Vol. 247,Spectroscopic Challenges of Photoionized Plasmas, ed. G. Ferland & D. W. Savin, 343Hubeny, I., & Lanz, T. 1995, ApJ, 439, 875Krtiˇcka, J., & Kub´at, J. 2010, A&A, 519, A50Krtiˇcka, J., Kub´at, J., & Krtiˇckov´a, I. 2015, A&A, 579, A111Krtiˇcka, J., & Kub´at, J. 2017, A&A, 606, A31Kroupa, P., 2001, MNRAS 322, 231K.Kudritzki, R. P., Pauldrach, A., Puls, J., & Abbott, D. C. 1989, A&A, 219, 205Kudritzki, R.-P., & Puls, J. 2000, Annual Review of Astronomy and Astrophysics, 38, 613 IBLIOGRAPHY Kudritzki, R. P. 2002, ApJ, 577, 389Kurucz, R. L. 1979, ApJS, 40, 1Lamers, H. & Cassinelli, J., Introduction to Stellar Winds , Cambridge University Press, 1999.Lanz, T., 2000, EAA Book E, 2102 L,