Coupling hydrodynamics with comoving frame radiative transfer
I. A unified approach for OB and WR stars
Institut für Physik und Astronomie, Universität Potsdam, KarlLiebknechtStr. 24/25, 14476 Potsdam, Germany
email: ansander@astro.physik.unipotsdam.de
Received: 16 February 2017
Accepted: 15 April 2017
Context. For more than two decades, stellar atmosphere codes have been used to derive the stellar and wind parameters of massive stars. Although they have become a powerful tool and sufficiently reproduce the observed spectral appearance, they can hardly be used for more than measuring parameters. One major obstacle is their inconsistency between the calculated radiation field and the wind stratification due to the usage of prescribed massloss rates and windvelocity fields.
Aims. We present the concepts for a new generation of hydrodynamically consistent nonlocal thermodynamical equilibrium (nonLTE) stellar atmosphere models that allow for detailed studies of radiationdriven stellar winds. As a first demonstration, this new kind of model is applied to a massive O star.
Methods. Based on earlier works, the PoWR code has been extended with the option to consistently solve the hydrodynamic equation together with the statistical equations and the radiative transfer in order to obtain a hydrodynamically consistent atmosphere stratification. In these models, the whole velocity field is iteratively updated together with an adjustment of the massloss rate.
Results. The concepts for obtaining hydrodynamically consistent models using a comovingframe radiative transfer are outlined. To provide a useful benchmark, we present a demonstration model, which was motivated to describe the wellstudied O4 supergiant ζPup. The obtained stellar and wind parameters are within the current range of literature values.
Conclusions. For the first time, the PoWR code has been used to obtain a hydrodynamically consistent model for a massive O star. This has been achieved by a profound revision of earlier concepts used for WolfRayet stars. The velocity field is shaped by various elements contributing to the radiative acceleration, especially in the outer wind. The results further indicate that for more dense winds deviations from a standard βlaw occur.
Key words: stars: massloss / stars: winds, outflows / stars: earlytype / stars: atmospheres / stars: fundamental parameters / stars: massive
© ESO, 2017
1. Introduction
In order to understand massive stars and their winds, stellar atmosphere models have become a powerful and widely used instrument. Typically applied for spectroscopic analysis, these models yield quantitative information on the stellar wind properties together with fundamental stellar parameters. The special conditions in stellar winds lead to the development of sophisticated model atmosphere codes performing several complex calculations. The outer layers are not even close to local thermodynamical equilibrium (LTE), requiring the population numbers to be calculated from a set of statistical equations. For a sufficient treatment, large model atoms with hundreds of levels in total have to be considered. Furthermore, the linedriven winds of hot stars demand a proper description of the radiative transfer in an expanding atmosphere.
Basically three different approaches exist to tackle the radiative transfer problem: the first are analytical descriptions, usually based on the concept of CAK theory (named after pioneering work of Castor et al. 1975), which allow for rapid calculation of the radiative acceleration at the cost of several approximations. While the theory has undergone several extensions relaxing the original assumptions (see, e.g., Friend & Abbott 1986; Pauldrach et al. 1986; Kudritzki et al. 1989; Gayley 1995; Puls et al. 2000; Kudritzki 2002), it has opened up the whole area of timedependent and even multidimensional calculations (e.g., Owocki et al. 1988; Feldmeier 1995; Owocki & Puls 1999; Dessart & Owocki 2005; Sundqvist et al. 2010). CAKlike concepts are therefore not only used in stellar atmosphere analyses, but also in most cases where a more detailed radiative transfer would be too costly from a computational standpoint; for example, detailed multidimensional hydrodynamical simulations of a stellar cluster or the interaction with a companion (e.g., Blondin et al. 1990; Manousakis et al. 2012; van Marle et al. 2012; Čechura & Hadrava 2015).
An alternative to (semi)analytical approximations is to calculate the radiative force with the help of Monte Carlo (MC) methods. Motivated already by the work of Lucy & Solomon (1970), this approach was first applied by Abbott & Lucy (1985) and later used for a variety of massloss studies (e.g., de Koter et al. 1993, 1997; Vink et al. 1999, 2000, 2001). More recently the application has been widely extended, including velocity field and clumping studies as well as multidimensional calculations (e.g., Müller & Vink 2008, 2014; Muijres et al. 2011, 2012; Šurlan et al. 2012; Noebauer & Sim 2015). Using the MC approach allows one to include effects such as multiple line scattering, but is computationally much more expensive than CAKlike calculations, especially in multidimensional approaches. It is therefore mostly used for fundamental studies and rarely applied when analyzing a particular stellar spectrum.
A third method to tackle the radiative transfer is the calculation in the comoving frame (CMF). Built on the conceptual work of Mihalas et al. (1975), it is essentially a bruteforce integration over the frequency range, using the advantage that opacity κ and emissivity η are isotropic in the comoving frame. Various studies using a CMF approach have been performed since the 1980s (e.g., Hamann 1980, 1981; Hillier 1987; Pauldrach et al. 1986; Sellmaier et al. 1993; Baron et al. 1996), culminating in the development of a handful of stellar atmosphere codes using the CMF radiative transfer either partially or exclusively, including phoenix (Hauschildt 1992; Hauschildt & Baron 1999), wmbasic (Pauldrach et al. 1994, 2001), fastwind (SantolayaRey et al. 1997; Puls et al. 2005), cmfgen (Hillier 1990b,a; Hillier & Miller 1998), and PoWR (Hamann 1985, 1986; Gräfener et al. 2002; Hamann & Gräfener 2003). For studying stars with denser winds and especially WolfRayet stars, almost all spectral analyses are performed with CMFbased atmosphere codes (e.g., Hillier & Miller 1999; Crowther et al. 2002; Hamann et al. 2006; Sander et al. 2012; Hainich et al. 2014). Due to the complexity of the CMF calculation, these codes typically assume spherical symmetry, allowing for a 1 or 1.5dimensional treatment, and a stationary wind situation.
Given that on top of the two major tasks, that is, solving the statistical equations and the radiative transfer, several further challenges exist, such as ironline blanketing or the need for a consistent calculation of the temperature stratification in an expanding, nonLTE environment, it is not that surprising that only a few codes exist that can adequately model a stellar atmosphere for a hot star with a dense, linedriven wind. So far, these codes are typically used for either measuring stellar and wind parameters or predicting fluxes and related quantities for a given set of parameters. Yet, only a few examples exist where they are used to actually predict the wind parameters. This lack of examples results from the fact that – at least in the wind part – most stellar atmosphere models use a prescribed velocity field instead of consistently calculating the wind stratification. While this approach is mostly sufficient for the current use of the atmosphere models, it also cuts off a variety of potential applications. In order to open up this perspective, we present a new approach for hydrodynamically consistent stellar atmosphere models using the Potsdam WolfRayet (PoWR) model atmosphere code. Originally starting from earlier efforts for WolfRayet (WR) stars, we have developed a brand new scheme to update the massloss rate and the velocity stratification that can finally be applied to both WR and OB models. For the first time, we will present a hydrodynamically consistent PoWR model for an O supergiant, closely reproducing most of the spectral features for ζ Pup.
In Sect. 2 of this work, we briefly discuss the underlying stationary wind hydrodynamics and introduce a special notation that is helpful in analyzing the status of PoWR models with respect to hydrodynamical consistency. A special emphasis is given in the following Sect. 3 on the meaning of the critical point. Afterwards, the basic concepts of the PoWR code and its set of input parameters are outlined in Sect. 4. Section 5 then deals with all the techniques used for obtaining a hydrodynamically consistent model before showing and discussing the results for an example model in Sect. 6. Finally, conclusions are drawn in Sect. 7.
2. Stationary wind hydrodynamics
In a hot stellar wind, which we here describe as a onedimensional, stationary outflow, the accelerations due to radiation and gas pressure have to balance gravity g(r) = GM_{∗}r^{2} and inertia . The corresponding equation of motion is therefore (1)with a_{rad} representing the total radiative acceleration, that is, using the same notation as in Sander et al. (2015). The term a_{press} describes the gas (and potentially turbulence) pressure, that is, (4)To replace the pressure P with the density, we use the equation of state for an ideal gas P(r) = ρ(r)·a^{2}(r) and define (5)with μ being the mean particle mass (including electrons) in units of the hydrogen atom mass m_{H}, T being the electron temperature, and ν_{mic} the microturbulence velocity, which is a free parameter in our models. For vanishing turbulence, a(r) is identical to the isothermal sound speed.
With the equation of state together with the equation of continuity (6)the a_{press}term in Eq. (1) can be rewritten in three terms that remove all explicit ρdependencies in favor of terms containing only ν, r, and a (see, e.g., Sander et al. 2015, for an explicit calculation). As a consequence, the hydrodynamic equation for a sphericallysymmetric wind reads with Γ_{rad}(r): = a_{rad}(r) /g(r). For an easier reading, we have dropped the explicit notation of the radius dependencies in Eqs. (7) and (8) and continue to do so unless absolutely necessary for the context. Using the two definitions one can write Eq. (7) in a more compact way: The quantities and are dimensionless and therefore ideal for visualizations. We note that in the subsonic regime, that is, ν≪ a, we obtain and thus Eq. (13) reduces to (14)which is a form of the hydrostatic equation discussed in a previous paper (Sander et al. 2015).
3. The critical point
The crucial difference between Eq. (14) and the full hydrodynamic equation (Eq. (13)) is the denominator . In contrast to the quasihydrostatic case, the hydrodynamic equation has a critical point at ν= a, that is, where the denominator becomes zero. In order to allow for a finite solution for the velocity gradient at the corresponding radius r_{c}, the nominator must also vanish at exactly this point. This leads to the constraint (15)As will be discussed below, this constraint implies a fixing of the massloss rate Ṁ, even though this quantity does not appear explicitly in the hydrodynamic Eq. (7). Since the hydrostatic equation does not have this constraint, it can provide a solution for ν(r) for any nonvanishing value of Ṁ.
Fig. 1 Dimensionless depthdependent quantities and are shown for a converged hydrodynamically consistent atmosphere model. The three terms adding up to are also indicated by the orange curves. At the critical point, both and become zero. 

Open with DEXTER 
An illustration of and for a converged hydrodynamically consistent model is shown in Fig. 1. Unlike in the CAKtype approaches, the critical point in Eq. (13) and Fig. 1 is identical to the sonic point, potentially corrected for a turbulence contribution. This is a direct consequence of our approach where we do not assume any semianalytical expression for the radiative acceleration a_{rad}. Instead, we simply treat a_{rad}(r) as a given quantity that can be expressed as a function of radius. Of course a_{rad} also reacts on changes of the velocity field and the massloss rate, but instead of trying to parametrize this in analytical form, we use an iterative approach and recalculate a_{rad}(r) after any adjustment of ν(r) or Ṁ until the thereby calculated acceleration does not enforce further velocity or massloss rate updates.
The PoWR code has already been used in the past for hydrodynamically consistent calculations of WR stars by Gräfener & Hamann (2005, 2008). In their approach, which turned out to be suitable only for thickwind WR stars, the critical point was not identical to the sonic point due to a semianalytical approach where the radiative acceleration was described with the help of an effective force multiplier parameter α(r). In the present work, we do not use such a parametrization and therefore the critical point in our hydrodynamic equation is identical to the sonic point. The implementation of the new method is also conceptually different from the one described in Gräfener & Hamann (2005) and will be further outlined in Sect. 5.
4. PoWR
4.1. Basic concepts
For the Potsdam WolfRayet (PoWR) model atmospheres we assume a spherically symmetric atmosphere with a stationary mass outflow. In order to properly describe the situation of an expanding atmosphere without the LTE approximation, the equations of statistical equilibrium and radiative transfer have to be solved iteratively until a consistent solution for the radiation field and the population numbers are obtained. In addition, the temperature stratification is updated iteratively to ensure energy conservation in the expanding atmosphere. This is performed using the improved UnsöldLucy method described in Hamann & Gräfener (2003) or alternatively via the electron thermal balance which has recently been added to the PoWR code (see Sander et al. 2015, and references therein).
If then all changes to the population numbers are smaller than a defined threshold, the atmosphere model is considered to be converged and the synthetic spectrum is calculated using a formal integration in the observer’s frame. The iron group elements are treated in a superlevel approach: the levels are grouped into energy bands, which are then represented by superlevels. While we assume LTE for the relative occupations of the individual levels inside a superlevel, the superlevels themselves are treated in full nonLTE. The detailed crosssections for the superlevel transtions have been prepared on a sufficiently fine frequency grid prior to the model iteration and contain all the individual transitions to ensure that radiative transfer treats each of these transitions at their proper frequency (see Gräfener et al. 2002, for more details.) PoWR is furthermore able to account for wind inhomogeneities in the socalled “microclumping” approximation (see Hamann & Koesterke 1998). In the calculation of the formal integral, PoWR can also account for optically thick clumps in an approximate way, see Oskinova et al. (2007) for details.
The radiative transfer is calculated in the CMF, thereby implicitly accounting for multiple scattering and avoiding all simplifications, which are done in the faster but more approximate concepts used, for example, in timedependent calculations. In particular, the solution is obtained by solving the moment equations via a differencing scheme based on the concepts of Mihalas et al. (1976a,b). The variable Eddington factors required in this scheme are obtained from a socalled “raybyray solution” where an angledependent radiation transfer is performed using shortcharacteristics integration (see Koesterke et al. 2002, for details). The CMF radiative transfer requires a strictly monotonic velocity field. Even in a stationary wind, this might not always be the case and therefore in some cases the solution for ν(r) resulting from the hydrodynamic equation cannot be applied. This will be discussed in more detail in a future paper and does not apply to the models presented in this work.
In the case of stars with low or moderate log g it is sufficient to assume that the intrinsic line profiles are Gaussians with a constant Doppler broadening velocity ν_{dop} during the CMF calculations while in the formal integral, where the emergent spectrum is eventually obtained, detailed thermal, microturbulent, and pressure broadening are accounted for with their depth dependence.
The necessary atomic data required for our calculations are taken from a variety of sources, combining Wiese et al. (1966), Eissner et al. (1974), Bashkin & Stoner (1975), the opacity project (OP, Cunto & Mendoza 1992), the Kurucz atomic database^{1}, the NIST atomic database^{2}, private communication with K. Butler, and several minor sources listed in Hamann et al. (1992). For argon, which turns out to be one of the more important elements driving the outer wind of our demonstration model, we use opacity project data combined with level energies from the NIST atomic database. The iron group elements, which are treated as one generic element with the help of the superlevel approach described in detail in Gräfener et al. (2002), are modeled using Kurucz data if available (usually up to ionization stage X), while opacity project data are applied to also cover the higher ions. The collisional crosssections are described with different formulae depending on the element and ion, most notably from Jefferies (1968, Eqs. (6.24), (6.25)), K. Butler (priv. comm.), and van Regemorter (1962, Eq. (22)). The latter is also used for all collisional transitions of argon and the iron group superlevels if the corresponding radiative transition is allowed. The collisional crosssections of forbidden transitions are mostly approximated by Mendoza (1983, Appendix 4), though for very few ions a more specialized treatment is applied (e.g., Berrington et al. 1982, for He i). For the boundfree transitions, we use Jefferies (1968, Eq. (6.39)) for all collisional ionizations while we branch between OP data fits, Mihalas (1967) and Seaton (1960), depending on the ion for the photoionization crosssections. The hydrogenic approximation (e.g., Cowan 1981) is widely used as fallback, which also applies for Ar and the iron group.
4.2. Model parameters
PoWR model atmospheres can be specified by a set of fundamental parameters. These are:

The chemical abundances of all considered elements, typicallygiven as mass fractions X_{i}.

Two out of the three quantities connected by StefanBoltzmann’s law, namely:

the stellar radius R_{∗}, defined at a specified Rosseland continuum optical depth τ_{∗} (default: τ_{∗} = 20);

the effective temperature T_{∗} at the radius R_{∗};

the luminosity .


The stellar mass M_{∗}, either given explicitly via input of M_{∗} or log g, or calculated from the luminosity if the stellar mass is not given otherwise. In the latter case, depending on the stellar type, the massluminosityrelations from Langer (1989) or Gräfener et al. (2011) are used.

The massloss rate Ṁ or an implying quantity (see below).

The terminal wind velocity ν_{∞} and the wind velocity law ν(r), directly implying the density stratification via the continuity Eq. (6).

The clump density contrast (cf. Hamann & Koesterke 1998, for a detailed description).
As an alternative to the massloss rate Ṁ, one can also specify a line emission measure in the form of either the transformed radius (16)(Schmutz et al. 1989; Hamann & Koesterke 1998, for the current form) or the wind strength parameter (17)(Puls et al. 1996, 2008, for the current form). Since all other quantities in Eqs. (16) and (17) have to be specified anyhow, these quantities imply a certain value of Ṁ. Using R_{t} or log Q_{ws} can be helpful when calculating model grids or searching for models with a similar emission line strength in their normalized spectra.
In hydrodynamically consistent models, Ṁ and ν(r) are adjusted in order to ensure that the hydrodynamic equation is fulfilled throughout the atmosphere. However, as they also define the density stratification, starting values are still required as an input for the calculations. Depending on the starting model, the resulting Ṁ and ν(r) of a converged hydrodynamically consistent model can differ significantly from their initial specifications.
4.3. Iteration scheme
With the main concepts and parameters given in the previous paragraphs, the overall iteration scheme for a PoWR model can now be summarized by the following steps:

1.
Model start: setup of radius and frequency grids, first velocitystratification, start approximation for the population numbersn(r) and the radiation field J_{ν}(r).

2.
Main iteration:

(a)
solution of the radiative transfer in the comoving frame;

(b)
temperature corrections (if necessary);

(c)
solution of the statistical equations;

(d)
(optional:) solution for the hydrostatic or hydrodynamic equation to update the velocity/density stratification.

(a)

3.
Formal integration: calculation of the emergent spectrum in the observer’s frame.
For efficiency, we use the method of variable Eddington factors, where the Eddington factors f_{ν} and g_{ν}, which have to be obtained from the more costly raybyray radiative transfer, are only updated every few iterations while otherwise only the momentum equations are solved in order to obtain the new radiation field and the radiative acceleration. For convenience, we schedule stratification updates immediately before the next renewal of the Eddington factors.
A sketch of this iteration scheme is given in Fig. 2. The stratification update, which can be either restricted to the quasihydrostatic part as outlined in Sander et al. (2015) or the full hydrodynamical update described in this work, is fully integrated into the main iteration. The particular details of the hydrodynamic stratification update are discussed in Sect. 5.
Fig. 2 Iteration scheme for the calculation of a PoWR atmosphere model. The inner cycle without the raybyray radiative transfer is typically applied for a few (typically 5–6) iterations before the Eddington factors must be renewed. If the criteria for a stratification update are fulfilled – see Sect. 5.4 for the case of a HD update – this is performed before the next raybyray transfer would be scheduled. On the right side of the figure, indicated by gray arrows, the most important quantities obtained by the current step are highlighted. 

Open with DEXTER 
5. Hydrodynamically consistent models
5.1. Start approximation
In order to integrate the hydrodynamic equation in the form of Eq. (13), the quantities a(r) and a_{rad}(r) have to be specified as a function of radius. This cannot be done from scratch and therefore a starting approximation for the stellar atmosphere has to be given, including a velocity stratification. Unless the new model is only a small variation of an already existing hydrodynamically consistent model, where one could employ the old velocity field, usually a model with a βlaw connected to a consistent hydrostatic solution (see Sander et al. 2015) is adopted as a starting approximation. For the massloss rate, it has turned out to be helpful, if at least the global energy budget is close to consistency. To obtain this budget, the hydrodynamic equation is written in the form (18)and then integrated over r and multiplied with Ṁ: (19)This last equation (Eq. (19)) describes the balance between the modeled wind luminosity L_{wind} and the provided power W_{wind}. Dividing Eq. (19) by L_{wind} yields the socalled work ratio (20)While stellar atmosphere models with Q < 1 do not provide a radiative acceleration that is sufficient to drive the wind, models with Q > 1 exhibit a radiation acceleration that could actually drive a stronger wind. Models with Q = 1 exactly supply the power that is required to drive the wind. The corresponding models are therefore consistent on a “global” scale, as they fulfill the integrated form of the hydrodynamic equation. Although they usually do not fulfill this equation locally, models with Q ≈ 1 are usually well suited as a starting model for the full hydrodynamic calculations. A similar approach to identify proper start approximations has been used by Gräfener & Hamann (2005) in their calculation of a hydrodynamically consistent WC atmosphere.
5.2. Obtaining a consistent velocity field
With a given starting model, all terms in the hydrodynamic equation are known, including the radiative acceleration a_{rad}(r) as a function of radius, and the terms and can be calculated. Since the ratio of the latter two essentially defines the right hand side of Eq. (13), both terms have to vanish at exactly the same radius r_{c}, defining the critical point of the equation. For the nonconsistent starting model, this is generally not the case and thus the radius will differ from the radius . In some cases, can become zero at more than one point, thus indicating a nonmonotonic solution for ν(r). While this does not prevent the integration of the hydrodynamic equation, a nonmonotonic ν(r) cannot be used in the CMF radiative transfer and thus such cases are discarded at the moment.
Fig. 3 As in Fig. 1, but now for a standard PoWR model where only the quasihydrostatic regime is treated selfconsistently. The dimensionless, depthdependent quantities and now cross zero at different radii and , which is typical for nonconsistent models and indicates a necessary adjustment of the massloss rate. 

Open with DEXTER 
Assuming that there is only one radius at which we have , this point acts as the current candidate for the critical point. Starting at with , Eq. (13) is then integrated inwards and outwards to obtain the new velocity field. By setting the wind velocity to the current value of a at , the critical point condition is automatically fulfilled and we achieve a velocity field that smoothly passes through the critical point. In principle it would also be possible to start the integration at or any point inbetween and , but this would require a modification of during the hydrodynamic iteration in order to fulfill the critical point condition. Several approaches have been tested during the development phase and none of them turned out to be favorable. Essentially, the approach by Gräfener & Hamann (2005, 2008) made use of modifying when changing the massloss rate. While this worked fine for some WR models, it turned out to fail for OB models. Furthermore, their approach required the calculation of a force multiplier parameter α(r), which can only be obtained by a modified radiative transfer calculation, thereby essentially doubling the calculation times for the radiative transfer before each hydro stratification update. On the other hand, modifications of the Γ_{rad}term in without any prediction on how the radiative acceleration will react on changes of ν(r) or Ṁ turn out not to be precise enough. Thus the direct integration from proved to be the best method, both in terms of stability and performance.
5.3. Calculation of the massloss rate
By starting the integration of the velocity field outwards from the critical point of the hydrodynamic equation, one automatically obtains the terminal velocity ν_{∞} when reaching the outer boundary, since ν_{∞} ≈ ν(R_{max}) as long as R_{max} is chosen to be sufficiently large. This method so far provides a new velocity field that fulfills the hydrodynamic equation, but does not perform any update of the massloss rate Ṁ. This is already a powerful tool, but the major issue with such a solution is their inconsistency with some of the initial stellar parameters; since the integration starts from the critical point, also the inner boundary value ν_{min}: = ν(R_{∗}) is not fixed, but instead obtained from integrating the hydrodynamic equation. As long as and are not identical before the integration, the total optical depth at R_{∗} will change after the velocity update. This especially means that T_{∗} and the corresponding R_{∗} in a converged model would refer to a different optical depth than in the starting model. As long as the total optical depth is larger than the old one, one could infer the values for the original optical depth. However, the consequence would be that the obtained hydrodynamic model refers to a different temperature and radius than the starting model, thereby being inconsistent with the radiative transfer calculation.
This problem can be solved by introducing another constraint, namely the conservation of the total optical depth. More precisely, for practical reasons we demand the conservation of the total Rosseland continuum optical depth, which we denote as τ_{Ross}(R_{∗}) throughout this work. As a consequence, the massloss rate Ṁ needs to be updated, but the main model parameters T_{∗} and R_{∗} now keep their intended reference. Unfortunately, finding a proper update method for Ṁ is not a trivial task. A simple approach would be to use the definition of the optical depth and replace the density via Eq. (6) to obtain an expression that explicitly contains Ṁ: However, even though Eq. (23) seems like a straightforward approach to extract the density and thus the massloss rate, one must keep in mind that ϰ_{Ross}(r) is not generally depthindependent, as some of the contributions (e.g., the boundfree opacities) do not just have a linear dependence on ρ(r). In fact, using this expression leads to large changes in Ṁ, often overpredicting the required changes by orders of magnitude. Since the critical point tends to change significantly even for moderate Ṁ updates, this method can only be successfully applied in very few cases and thus cannot be considered as a standard approach. The sensitivity of the critical point was already found by Pauldrach et al. (1986) when they used an early form of this concept with the Thomson opacity instead of the Rosseland continuum opacity as they did not account for the freefree and boundfree continuum in their radiative force. Similar problems occur when employing other quantities with analog descriptions, such as the integrated density without the mass absorption coefficient (24)or the total atmosphere mass (25)In order to obtain a more stable method, this work utilizes a completely different approach, where we approximate the response of the radiative acceleration a_{rad} to a change of the massloss rate by a factor f. A typical example for an Ostar model is shown in Fig. 4, where the ratio of unmodified to modified acceleration is shown for different values of f^{1}. As the radiative acceleration is defined as there is a leading dependence with Ṁ^{1}, therefore making f^{1} the more interesting quantity for the plots. The lower panel of Fig. 4 also illustrates that the particular value of f changes the amplitude, but not the general behavior of the response, since one can scale all the results almost perfectly to the same curve resp(r) using the relation (28)With the help of Eq. (28) it would therefore be possible to implement the detailed response of a_{rad}(r) to suggested changes of Ṁ in order to improve the calculation of the massloss rate in a hydro iteration. However, similar to what has been discussed for the α(r)approach from Gräfener & Hamann (2005), this would double the CMF calculation time before each update of the velocity field. Furthermore, the relation (28) cannot account for the typically more complex changes of ν(r) and thus even a massloss rate obtained in this detailed way does not lead to a better model convergence. In fact, such methods have been tested and have turned out not to be better than the more approximate way described below.
Fig. 4 Upper panel: response of the radiative acceleration a_{rad}(r) to a change of the massloss rate Ṁ is illustrated by plotting the ratio of the modified to unmodified acceleration for different modification factors f. Lower panel: results for different factors f can be scaled to an almost unique curve. 

Open with DEXTER 
Fig. 5 Upper panel: response of the radiative line (solid curves) and continuum (dashed curves) acceleration for two different change factors f of the massloss rate Ṁ. Similar to Fig. 4, the ratio of the modified to unmodified acceleration is plotted. Lower panel: zoomin of the curves showing the continuum response. 

Open with DEXTER 
Apart from calculating the total response of the radiative acceleration in our test calculations, we also calculated the isolated response of the line and the total continuum term. In Fig. 5 the results for two different massloss modification factors f are shown. As we can see, the continuum shows only a very small response that can be neglected, while the line response is stronger and roughly of the order f^{1}. In fact f^{1} is never completely reached, especially not in the outer wind, but since we want to use our approximation only for the calculation of the massloss rate, this is not a problem. By slightly overpredicting the effect of the change in Ṁ, we avoid potential “overshooting” of the correction. We therefore assume for our calculations of the massloss rate update, that the radiative acceleration a_{rad} changes for a massloss rate modified by the factor f such that (29)with a_{rad}(Ṁ) = a_{lines}(Ṁ) + a_{cont}(Ṁ) denoting the acceleration with the original massloss rate. Using this assumption, we calculate the factor f which would be necessary to obtain at the current . For a converged model, this is automatically fulfilled and we obtain f = 1, that is, no more change in Ṁ. In the general case, the condition leads to (30)with Γ_{lines} = a_{lines}/g and Γ_{cont} = a_{cont}/g. To avoid overly large corrections, which could significantly disturb the model convergence, only 50% of the calculated change for Ṁ is usually applied in one iteration.
The described update of the massloss rate now leads to a convergence of and . However, there is so far no guarantee that the total optical depth of the converged model will be identical to that of the starting model. To ensure also the latter, the whole atmosphere stratification is radially adjusted before integrating the hydrodynamic equation. This adjustment is performed already with the new massloss rate and ensures the conservation of the total Rosseland continuum optical depth.
5.4. Scheme of the hydrodynamic stratification update
The implementation concept of the hydrodynamical stratification update in regards to the overall model iteration is similar to the one described in Sander et al. (2015) for models with a consistent quasihydrostatic part. In fact, the full hydrodynamic treatment and quasihydrostatic update are alternative branches for the step (2d) in the overall iteration scheme outlined in Sect. 4.3. In order to ensure that the overall model calculations are not vastly disrupted by a stratification update, such updates are not performed during each iteration, but only immediately before the Eddington factors in the following radiative transfer job are recalculated and if the overall corrections to the population numbers are below a certain prespecified level (in addition, it is also possible to require a certain level of flux consistency for a stratification update).
The hydrodynamic stratification update itself can be described by the following scheme:

1.
Check whether the hydrodynamic equation is fulfilled at alldepth points: if so, no update needs to be performed.

2.
The quantities Γ_{rad}(r) and a(r) are calculated based on the current model stratification.

3.
If , the mass loss rate Ṁ is updated by a factor f as described in Eq. (30).

4.
Iteration:

(a)
starting from , the new velocity field is obtained via integrating Eq. (13) with a fourthorder KungeKutta method using adaptive step sizes. The quantities ℱ and are calculated on the fly. To avoid numerical issues near the critical point, we make use of l’Hôpital’s rule. As the hydrodynamic integration requires a resolution below the regular depth grid, we use spline interpolation to obtain the interstice values for Γ_{rad} and a;

(b)
with the new ν(r) now given, we calculate the resulting new density stratification and total Rosseland continuum optical depth τ_{Ross}(R_{∗});

(c)
if the τ_{Ross}(R_{∗}) is conserved, the iteration ends. Otherwise Γ_{rad}(r) and a(r) are shifted radially and the next iteration cycle is started.

(a)

5.
If necessary, the depth grid spacing is updated to better reflect the new stratification. All necessary quantities are interpolated from the old to the new grid.
After the hydrodynamic stratification update is complete, the overall iteration cycle continues with the next radiative transfer calculation (we refer also to Sect. 4.3 for details of the overall iteration). Due to the inclusion in the overall iteration cycle, the values of Γ_{rad}(r) and a(r) used in the next hydrodynamic stratification update implicitly include all effects of the former stratification and massloss rate update. While this procedure can significantly increase the total number of overall iterations compared to nonHD models, this iterative approach allows us to refrain from any (semi)analytical assumptions for the radiative acceleration, making this method essentially applicable for the whole range of stars that can be described by PoWR atmosphere models. The atmosphere model is eventually considered to be converged if all of the following requirements are met:

The relative corrections to the population numbers are below acertain level (typically 10^{3}).

Flux consistency is achieved within a certain accuracy (typically setting: relative departures may not be larger than 10^{2}).

The hydrodynamical equation is fulfilled throughout the whole atmosphere within a specified accuracy (typically 5%, but we allow larger deviations at the inner and outer boundary).
With the exception of the last point of course, these criteria are the same as for nonHD models, which are used for reproducing observed spectra and empirically obtain stellar and wind parameters. Based on the converged atmosphere model, the emergent spectrum is subsequently calculated in the observer’s frame, allowing us to crosscheck the results with observed spectra.
6. Results
In order to test whether or not our new method is applicable to OB stars, where the approach from Gräfener & Hamann (2008) failed, we calculated a hydrodynamically consistent model for the wellstudied O supergiant ζ Pup/HD 66811. Starting from the parameters given in Bouret et al. (2012), we first calculated a standard PoWR model using a prescribed βlaw connected to a consistent quasihydrostatic part as described in Sander et al. (2015). While a model reproducing most spectral features already requires Ne, Mg, Si, P, and S to be considered, the work ratio of such a model is only Q = 0.74. By adding further elements, most notably Ar, the work ratio was close to unity and the model could be used as a starting approach for the hydrodynamic calculations.
In our first approach, we applied the same clumping stratification as Bouret et al. (2012), that is, depthdependent (micro)clumping with a maximum value of D_{∞} = 20 or and no interclump medium. This is a standard approach in stateoftheart atmosphere models for hot and massive stars (e.g., Hamann & Koesterke 1998; Hillier & Miller 1999) and allows one to calculate the population numbers for the clumped wind, which has a density increased by a factor D(r) compared to a smooth wind. For the radiative transfer, on the other hand, one can average between the clump and interclump medium as clumps are assumed to have a small size in comparison to the meanfree path of the photons. Furthermore, instead of the standard description of depthdependent clumping in PoWR (Gräfener & Hamann 2005), we employed the same parametrization as in the CMFGEN model from Bouret et al. (2012), namely (31)introduced in Martins et al. (2009), where the clumping “onset” is described by a velocity ν_{cl}. In their analysis for ζ Pup, Bouret et al. (2012) use ν_{cl} = 100 km s^{1}. We started with a similar stratification, but quickly realized that the value for ν_{cl} is not sufficient and leads to solutions with an overly high terminal velocity together with an overly low massloss rate. Stratifications where we set ν_{cl} = 0.5 ν_{sonic} lead to better results, but since the sonic point can change during the iterations, we decided to implement another clumping stratification with (32)that is, where we specify the clumping onset via an optical depth τ_{cl} instead of a particular velocity. This approach was eventually applied in the final model presented in this work. As we discuss later on, it was furthermore necessary to reduce the maximum clumping value to D_{∞} = 10 in our model. The complete set of input parameters for the final hydrodynamical model is compiled in Table 1.
Input parameters for the ζ Pup model.
Atomic data used in the hydrodynamic model.
In models with a predefined velocity law in the wind part, it is usually sufficient to include only those elements that can either be seen in the spectrum or contribute significantly to the blanketing. However, for the hydrodynamic models, it is essential to include all ions that have a significant contribution to the radiative force, even if they neither leave a noticable imprint in the spectrum nor significantly affect the blanketing. Yet, accounting for all elements and their ions from hydrogen up to the iron group would be numerically extremely expensive and thus practically impossible. Fortunately, various elements and ions are only important in a certain parameter regime and thus can be neglected outside of these. A list of the ions considered in the model presented in this work is given in Table 2.
The acceleration balance of the hydrodynamically consistent model is shown in Fig. 6. Throughout the atmosphere, an excellent agreement between the outward and inward forces is obtained. Up until the critical point, not only the line acceleration and the Thomson term are important, but there are also significant contributions from the gas pressure and the true continuum, that is, the continuum not produced by Thomson scattering, to the driving. In the wind part, both of the latter terms become negligible. However, it has to be noted that this is not necessarily the case for all kinds of hot stars. In the more dense WolfRayet winds, situations can occur where the true continuum is not negligible in the wind (see, e.g., the consistent model for WR 111 in Gräfener & Hamann 2005). The importance of the pressure term strongly depends on the assumptions for microturbulence. In this model, a constant value of ν_{mic} = 15 km s^{1} was used in the hydrodynamic calculations. When using larger values or depthdependent descriptions with ν_{mic} increasing outwards, the a_{press}term can become significant again in the outer wind.
The resulting velocity field for the converged hydrodynamic models is shown in Fig. 7, where it is compared to the stratification of two standard models using a prescribed velocity field in the form of socalled βlaws, that is, (33)When implementing the βlaw into a stellar atmosphere model, Eq. (33) is usually slightly modified due to several reasons, most prominently the necessity to connect the wind domain with a proper quasihydrostatic domain (see Sander et al. 2015, for more details) and the numerical issues that would occur for ν(R_{∗}) = 0. This modification can be done in more than one way and thus also differ between different stellar atmosphere codes. In PoWR, two ways of “fine parametrization” for the βlaw are available, namely (34)and (35)with their fine parameters p and R_{s} or f_{s}, respectively. While R_{s} or f_{s} are responsible for a proper connection of the wind and the quasihydrostatic domain, p ≈ ν_{∞} ensures that the specified terminal velocity is reached at the outer boundary. All fine parameters are automatically calculated depending on the choices of the velocity field, the connection criterion and whether the parametrization from Eq. (34) or from Eq. (35) should be used.
Fig. 6 Acceleration stratification for a hydrodynamically consistent model. The wind acceleration (thick red diamond line) is compared to the repulsive sum of inertia and gravitational acceleration g(r) (black line). The fact that these two curves are (almost) identical illustrates that the hydrodynamic equation is fulfilled throughout the stellar atmosphere. For a more convenient illustration, all terms are normalized to g(r). The input parameters of the model are compiled in Table 1 while the resulting quantities can be found in Table 3. 

Open with DEXTER 
Fig. 7 Normalized velocity field for the hydrodynamicallyconsistent model (red) versus a model using a βlaw connected to a consistent quasihydrostatic part. The upper panel shows the velocity in nonlogarithmic units, thereby highlighting the wind part, while the lower panel displays the normalized velocity in logarithmic units, thus focusing on the inner layers. 

Open with DEXTER 
For the hydrodynamically consistent model presented in this work, there is of course no prescribed wind velocity field, but for the comparison calculations we had to make a choice and used Eq. (34) since it is more widely used in modern PoWR models. Interestingly the velocity field in the lower part of the wind, just above the sonic point, can be approximated with a beta law using β = 2.4, but the outer part of the wind is best matched with β = 0.9. Due to the fact that the fine parametrization is not unique as described above, these deduced βlaw approximations can vary slightly (10 to 20%) when using different fine parametrizations. Inbetween the two parts there is a steep increase of the velocity, steeper than could be modeled by a βlaw connected to the quasihydrostatic part. The reasons for this kind of velocity field are revealed when looking at the particular contributions to the radiative acceleration plotted in Fig. 8. Around and shortly above the critical point, only the iron group elements and the electron scattering contribute significantly to the radiative acceleration. In contrast, further outwards, many more elements contribute, and for r> 1.6 R_{∗} N, O, and Ar start to exceed not only Γ_{e}, but also the contribution of the iron group elements. S, Cl, and Ne follow further out and for r ≳ 2 R_{∗} also the contributions of C and P are comparable to Γ_{e}, having even a bit more impact than the iron group at this distance.
The complex contribution to the radiative force from the various elements is directly imprinted in the resulting velocity field and thus “naturally” explains the deviations from a standard βlaw which has been derived from a taylored fit of the UV resonance line profiles already by Hamann (1980). While we do not see a noticeable velocity plateau in the final model, such a plateau occurred in several of the test models derived during the preparation of this work. The plateau becomes visible if the iron group contribution already exceeds Γ_{e} around or even below the sonic point, which can already happen for slightly higher massloss rates than derived for ζ Pup here.
Fig. 8 Absolute (upper panel) and relative (lower panel) contributions to the radiative acceleration from the different elements considered in the hydrodynamically consistent atmosphere model. The total radiative acceleration and the acceleration due to gas pressure are also shown in the upper panel for comparison. The lower panel shows the extent to which electron scattering (pink solid curve) and the various elements contribute to the total radiative acceleration. 

Open with DEXTER 
Fig. 9 Comparison of the hydrodynamically consistent Ostar model with observed spectra of the O4 supergiant ζ Pup. The uppermost panel shows the spectral energy distribution while the other panels compare the normalized line spectra between the model (red) and the observations (blue) in various wavelength ranges. 

Open with DEXTER 
The derived parameters of our hydrodynamically consistent model are compiled in Table 3, while the spectral energy distribution and important parts of the normalized spectrum are compared to observations in Fig. 9. The UV observation was obtained with the IUE satellite (SWP15296) while the optical spectrum stems from an earlier observation within our group (Hamann, priv. comm.). The photometric data used in the SED plot have been taken from Ducati (2002). While the displayed model is rather to demonstrate the new technique described in the previous section and therefore has not been finetuned to precisely reproduce the spectral features of ζ Pup, one can still compare the main parameters to spectral analyses. While the starting parameters were motivated by the nonhydrodynamical model results from Bouret et al. (2012), their high clumping factor of 20 would lead to an underprediction of the Hα electron scattering wings and thus was reduced to the more typical value of 10. The massloss rate of log Ṁ = −5.8 is slightly lower than in Bouret et al. (2012), but higher than the value obtained by Pauldrach et al. (2012).
The emergent spectrum of our hydrodynamical model shows all the typical features of an early Oftype star with a fast wind and strong emission in both N iii λ46344042 and He ii λ4686 (see, e.g., Sota et al. 2011, for a classification scheme). When comparing the detailed spectral appearance to the observation of ζ Pup in Fig. 9, one can conclude that the UV spectrum, including the iron forest, is well reproduced apart from the precise shape of the nitrogen profiles, which can be affected by various parameters including socalled “superionization” due to Xrays, which were not included in our model. The optical spectrum reveals that the massloss rate might be slightly too high to reproduce this observation, since emission in Hα and the other prominent emission lines appears slightly too strong. Also Hβ and Hγ seem to be filled up by wind emission. The overall appearance, however, as well as the spectral energy distribution, are nicely reproduced, illustrating that we have a realistic stellar atmosphere model that would require only minor parameter adjustments to allow for a more detailed discussion. Unfortunately even these minor changes can lead to a significant amount of calculation effort when constructing hydrodynamically consistent models, which is why we refrain from further efforts in this introductory paper.
7. Conclusions
In this work we constructed the first hydrodynamically consistent PoWR model for an O supergiant. A new method for the consistent solution of the hydrodynamic equation, together with the solution of the statistical equations, the temperature stratification, and the radiative transfer has been developed and successfully applied. This new technique enables us to construct a new generation of PoWR models where the velocity field and the massloss rate are calculated consistently. To obtain the velocity field, the hydrodynamic equation is integrated inwards and outwards from the critical point. Since we provide the radiative acceleration calculated in the comoving frame as a function of radius, the critical point in our hydrodynamic equation is identical to the sonic point. The uniqueness of the critical point also provides the necessary condition to obtain the massloss rate.
As we calculate the velocity field from the hydrodynamic equation, it is mandatory to include all ions that significantly contribute to the radiative acceleration somewhere in the atmosphere. In the case of our demonstration model, especially the inclusion of Ar was crucial as it provides a major contribution to the driving in the outer wind, comparable to N and O, although it does not leave detectable features in the spectral ranges typically observed.
In the region around the critical point, the most important line driving contribution stems from the iron group elements. Although these elements turn out to be important contributors throughout the whole atmosphere for our demonstration model, there are several other elements exceeding their input in the outer part, namely N, O, S, Ar, Cl, C, P and partly Ne. Further followup calculations for a wider parameter range will be necessary to shed light on details; namely, which ions are responsible, and how this picture will change when transitioning to different massloss or temperature regimes.
The obtained velocity field cannot be approximated by a βlaw. The resulting massloss rate of our hydrodynamic model is in the range of what has been determined by empirical analyses for the O4 supergiant ζ Pup, and the resulting spectrum resembles the observed line spectrum and the spectral energy distribution (cf. Fig. 9). Our calculations confirm that the value of the massloss rate crucially depends on the location of the critical point, which in turn reacts to several factors, such as the assumed microturbulence, the onset of clumping, and the Fe abundance. Followup research will therefore be required in order to study the precise influence of these and other parameters.
Acknowledgments
We would like to thank the anonymous referee for the fruitful suggestions that helped to improve this paper. We would also like to acknowledge helpful discussions with D. John Hillier. The first author of this work (A.S.) is supported by the Deutsche Forschungsgemeinschaft (DFG) under grant HA 1455/26. T.S. is grateful for financial support from the Leibniz Graduate School for Quantitative Spectroscopy in Astrophysics, a joint project of the Leibniz Institute for Astrophysics Potsdam (AIP) and the Institute of Physics and Astronomy of the University of Potsdam. This research made use of the SIMBAD and VizieR databases, operated at CDS, Strasbourg, France.
References
 Abbott, D. C., & Lucy, L. B. 1985, ApJ, 288, 679 [NASA ADS] [CrossRef] [Google Scholar]
 Asplund, M., Grevesse, N., Sauval, A. J., & Scott, P. 2009, ARA&A, 47, 481 [NASA ADS] [CrossRef] [Google Scholar]
 Baron, E., Hauschildt, P. H., & Mezzacappa, A. 1996, MNRAS, 278, 763 [NASA ADS] [CrossRef] [Google Scholar]
 Bashkin, S., & Stoner, J. O. 1975, Atomic energy levels and Grotrian Diagrams, Vol. 1, Hydrogen I, Phosphorus XV, Vol. 2, Sulfur I, Titanium XXII (Amsterdam: NorthHolland Publ. Co.) [Google Scholar]
 Berrington, K. A., Fon, W. C., & Kingston, A. E. 1982, MNRAS, 200, 347 [NASA ADS] [CrossRef] [Google Scholar]
 Blondin, J. M., Kallman, T. R., Fryxell, B. A., & Taam, R. E. 1990, ApJ, 356, 591 [NASA ADS] [CrossRef] [Google Scholar]
 Bouret, J.C., Hillier, D. J., Lanz, T., & Fullerton, A. W. 2012, A&A, 544, A67 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Castor, J. I., Abbott, D. C., & Klein, R. I. 1975, ApJ, 195, 157 [NASA ADS] [CrossRef] [Google Scholar]
 Čechura, J., & Hadrava, P. 2015, A&A, 575, A5 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Cowan, R. D. 1981, The theory of atomic structure and spectra (Berkeley: University of California Press) [Google Scholar]
 Crowther, P. A., Dessart, L., Hillier, D. J., Abbott, J. B., & Fullerton, A. W. 2002, A&A, 392, 653 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Cunto, W., & Mendoza, C. 1992, Rev. Mexicana Astron. Astrofis., 23 [Google Scholar]
 de Koter, A., Schmutz, W., & Lamers, H. J. G. L. M. 1993, A&A, 277, 561 [NASA ADS] [Google Scholar]
 de Koter, A., Heap, S. R., & Hubeny, I. 1997, ApJ, 477, 792 [NASA ADS] [CrossRef] [Google Scholar]
 Dessart, L., & Owocki, S. P. 2005, A&A, 437, 657 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Ducati, J. R. 2002, VizieR Online Data Catalog, 2237 [Google Scholar]
 Eissner, W., Jones, M., & Nussbaumer, H. 1974, Comput. Phys. Commun., 8, 270 [NASA ADS] [CrossRef] [Google Scholar]
 Feldmeier, A. 1995, A&A, 299, 523 [NASA ADS] [Google Scholar]
 Friend, D. B., & Abbott, D. C. 1986, ApJ, 311, 701 [NASA ADS] [CrossRef] [Google Scholar]
 Gayley, K. G. 1995, ApJ, 454, 410 [NASA ADS] [CrossRef] [Google Scholar]
 Gräfener, G., & Hamann, W.R. 2005, A&A, 432, 633 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Gräfener, G., & Hamann, W.R. 2008, A&A, 482, 945 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Gräfener, G., Koesterke, L., & Hamann, W.R. 2002, A&A, 387, 244 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Gräfener, G., Vink, J. S., de Koter, A., & Langer, N. 2011, A&A, 535, A56 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Hainich, R., Rühling, U., Todt, H., et al. 2014, A&A, 565, A27 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Hamann, W.R. 1980, A&A, 84, 342 [NASA ADS] [Google Scholar]
 Hamann, W.R. 1981, A&A, 93, 353 [NASA ADS] [Google Scholar]
 Hamann, W.R. 1985, A&A, 148, 364 [NASA ADS] [Google Scholar]
 Hamann, W.R. 1986, A&A, 160, 347 [NASA ADS] [Google Scholar]
 Hamann, W.R., & Gräfener, G. 2003, A&A, 410, 993 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Hamann, W., & Koesterke, L. 1998, A&A, 335, 1003 [NASA ADS] [Google Scholar]
 Hamann, W.R., Leuenhagen, U., Koesterke, L., & Wessolowski, U. 1992, A&A, 255, 200 [NASA ADS] [Google Scholar]
 Hamann, W., Gräfener, G., & Liermann, A. 2006, A&A, 457, 1015 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Hauschildt, P. H. 1992, J. Quant. Spec. Radiat. Transf., 47, 433 [NASA ADS] [CrossRef] [Google Scholar]
 Hauschildt, P. H., & Baron, E. 1999, J. Comput. Appl. Math., 109, 41 [NASA ADS] [CrossRef] [Google Scholar]
 Hillier, D. J. 1987, ApJS, 63, 947 [NASA ADS] [CrossRef] [Google Scholar]
 Hillier, D. J. 1990a, A&A, 231, 116 [NASA ADS] [Google Scholar]
 Hillier, D. J. 1990b, A&A, 231, 111 [NASA ADS] [Google Scholar]
 Hillier, D. J., & Miller, D. L. 1998, ApJ, 496, 407 [NASA ADS] [CrossRef] [Google Scholar]
 Hillier, D. J., & Miller, D. L. 1999, ApJ, 519, 354 [NASA ADS] [CrossRef] [Google Scholar]
 Jefferies, J. T. 1968, Spectral line formation (Waltham, Mass.: Blaisdell, 1968) [Google Scholar]
 Koesterke, L., Hamann, W.R., & Gräfener, G. 2002, A&A, 384, 562 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Kudritzki, R. P. 2002, ApJ, 577, 389 [NASA ADS] [CrossRef] [Google Scholar]
 Kudritzki, R.P., & Puls, J. 2000, ARA&A, 38, 613 [NASA ADS] [CrossRef] [Google Scholar]
 Kudritzki, R. P., Pauldrach, A., Puls, J., & Abbott, D. C. 1989, A&A, 219, 205 [NASA ADS] [Google Scholar]
 Langer, N. 1989, A&A, 210, 93 [NASA ADS] [Google Scholar]
 Lucy, L. B., & Solomon, P. M. 1970, ApJ, 159, 879 [NASA ADS] [CrossRef] [Google Scholar]
 Manousakis, A., Walter, R., & Blondin, J. M. 2012, A&A, 547, A20 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Martins, F., Hillier, D. J., Bouret, J. C., et al. 2009, A&A, 495, 257 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Mendoza, C. 1983, in Planetary Nebulae, ed. D. R. Flower, IAU Symp., 103, 143 [NASA ADS] [CrossRef] [Google Scholar]
 Mihalas, D. 1967, Methods in Computational Physics, 7, 1 [NASA ADS] [Google Scholar]
 Mihalas, D., Kunasz, P. B., & Hummer, D. G. 1975, ApJ, 202, 465 [NASA ADS] [CrossRef] [Google Scholar]
 Mihalas, D., Kunasz, P. B., & Hummer, D. G. 1976a, ApJ, 206, 515 [NASA ADS] [CrossRef] [Google Scholar]
 Mihalas, D., Kunasz, P. B., & Hummer, D. G. 1976b, ApJ, 210, 419 [NASA ADS] [CrossRef] [Google Scholar]
 Muijres, L. E., de Koter, A., Vink, J. S., et al. 2011, A&A, 526, A32 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Muijres, L. E., Vink, J. S., de Koter, A., Müller, P. E., & Langer, N. 2012, A&A, 537, A37 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Müller, P. E., & Vink, J. S. 2008, A&A, 492, 493 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Müller, P. E., & Vink, J. S. 2014, A&A, 564, A57 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Noebauer, U. M., & Sim, S. A. 2015, MNRAS, 453, 3120 [NASA ADS] [CrossRef] [Google Scholar]
 Oskinova, L. M., Hamann, W.R., & Feldmeier, A. 2007, A&A, 476, 1331 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Owocki, S. P., & Puls, J. 1999, ApJ, 510, 355 [NASA ADS] [CrossRef] [Google Scholar]
 Owocki, S. P., Castor, J. I., & Rybicki, G. B. 1988, ApJ, 335, 914 [NASA ADS] [CrossRef] [Google Scholar]
 Pauldrach, A., Puls, J., & Kudritzki, R. P. 1986, A&A, 164, 86 [NASA ADS] [Google Scholar]
 Pauldrach, A. W. A., Kudritzki, R. P., Puls, J., Butler, K., & Hunsinger, J. 1994, A&A, 283, 525 [NASA ADS] [Google Scholar]
 Pauldrach, A. W. A., Hoffmann, T. L., & Lennon, M. 2001, A&A, 375, 161 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Pauldrach, A. W. A., Vanbeveren, D., & Hoffmann, T. L. 2012, A&A, 538, A75 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Puls, J., Kudritzki, R.P., Herrero, A., et al. 1996, A&A, 305, 171 [NASA ADS] [Google Scholar]
 Puls, J., Springmann, U., & Lennon, M. 2000, A&AS, 141, 23 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Puls, J., Urbaneja, M. A., Venero, R., et al. 2005, A&A, 435, 669 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Puls, J., Vink, J. S., & Najarro, F. 2008, A&ARv, 16, 209 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Sander, A., Hamann, W.R., & Todt, H. 2012, A&A, 540, A144 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Sander, A., Shenar, T., Hainich, R., et al. 2015, A&A, 577, A13 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 SantolayaRey, A. E., Puls, J., & Herrero, A. 1997, A&A, 323, 488 [NASA ADS] [Google Scholar]
 Schmutz, W., Hamann, W.R., & Wessolowski, U. 1989, A&A, 210, 236 [NASA ADS] [Google Scholar]
 Seaton, M. J. 1960, Rep. Prog. Phys., 23, 313 [NASA ADS] [CrossRef] [Google Scholar]
 Sellmaier, F., Puls, J., Kudritzki, R. P., et al. 1993, A&A, 273, 533 [NASA ADS] [Google Scholar]
 Sota, A., Maíz Apellániz, J., Walborn, N. R., et al. 2011, ApJS, 193, 24 [NASA ADS] [CrossRef] [Google Scholar]
 Sundqvist, J. O., Puls, J., & Feldmeier, A. 2010, A&A, 510, A11 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Šurlan, B., Hamann, W.R., Kubát, J., Oskinova, L. M., & Feldmeier, A. 2012, A&A, 541, A37 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 van Marle, A. J., Meliani, Z., & Marcowith, A. 2012, A&A, 541, L8 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 van Regemorter, H. 1962, ApJ, 136, 906 [NASA ADS] [CrossRef] [Google Scholar]
 Vink, J. S., de Koter, A., & Lamers, H. J. G. L. M. 1999, A&A, 350, 181 [NASA ADS] [Google Scholar]
 Vink, J. S., de Koter, A., & Lamers, H. J. G. L. M. 2000, A&A, 362, 295 [NASA ADS] [Google Scholar]
 Vink, J. S., de Koter, A., & Lamers, H. J. G. L. M. 2001, A&A, 369, 574 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Wiese, W. L., Smith, M. W., & Glennon, B. M. 1966, Atomic transition probabilities. Vol.: Hydrogen through Neon. A critical data compilation (NSRDSNBS 4, Washington, D.C.: US Department of Commerce, National Buereau of Standards) [Google Scholar]
All Tables
All Figures
Fig. 1 Dimensionless depthdependent quantities and are shown for a converged hydrodynamically consistent atmosphere model. The three terms adding up to are also indicated by the orange curves. At the critical point, both and become zero. 

Open with DEXTER  
In the text 
Fig. 2 Iteration scheme for the calculation of a PoWR atmosphere model. The inner cycle without the raybyray radiative transfer is typically applied for a few (typically 5–6) iterations before the Eddington factors must be renewed. If the criteria for a stratification update are fulfilled – see Sect. 5.4 for the case of a HD update – this is performed before the next raybyray transfer would be scheduled. On the right side of the figure, indicated by gray arrows, the most important quantities obtained by the current step are highlighted. 

Open with DEXTER  
In the text 
Fig. 3 As in Fig. 1, but now for a standard PoWR model where only the quasihydrostatic regime is treated selfconsistently. The dimensionless, depthdependent quantities and now cross zero at different radii and , which is typical for nonconsistent models and indicates a necessary adjustment of the massloss rate. 

Open with DEXTER  
In the text 
Fig. 4 Upper panel: response of the radiative acceleration a_{rad}(r) to a change of the massloss rate Ṁ is illustrated by plotting the ratio of the modified to unmodified acceleration for different modification factors f. Lower panel: results for different factors f can be scaled to an almost unique curve. 

Open with DEXTER  
In the text 
Fig. 5 Upper panel: response of the radiative line (solid curves) and continuum (dashed curves) acceleration for two different change factors f of the massloss rate Ṁ. Similar to Fig. 4, the ratio of the modified to unmodified acceleration is plotted. Lower panel: zoomin of the curves showing the continuum response. 

Open with DEXTER  
In the text 
Fig. 6 Acceleration stratification for a hydrodynamically consistent model. The wind acceleration (thick red diamond line) is compared to the repulsive sum of inertia and gravitational acceleration g(r) (black line). The fact that these two curves are (almost) identical illustrates that the hydrodynamic equation is fulfilled throughout the stellar atmosphere. For a more convenient illustration, all terms are normalized to g(r). The input parameters of the model are compiled in Table 1 while the resulting quantities can be found in Table 3. 

Open with DEXTER  
In the text 
Fig. 7 Normalized velocity field for the hydrodynamicallyconsistent model (red) versus a model using a βlaw connected to a consistent quasihydrostatic part. The upper panel shows the velocity in nonlogarithmic units, thereby highlighting the wind part, while the lower panel displays the normalized velocity in logarithmic units, thus focusing on the inner layers. 

Open with DEXTER  
In the text 
Fig. 8 Absolute (upper panel) and relative (lower panel) contributions to the radiative acceleration from the different elements considered in the hydrodynamically consistent atmosphere model. The total radiative acceleration and the acceleration due to gas pressure are also shown in the upper panel for comparison. The lower panel shows the extent to which electron scattering (pink solid curve) and the various elements contribute to the total radiative acceleration. 

Open with DEXTER  
In the text 
Fig. 9 Comparison of the hydrodynamically consistent Ostar model with observed spectra of the O4 supergiant ζ Pup. The uppermost panel shows the spectral energy distribution while the other panels compare the normalized line spectra between the model (red) and the observations (blue) in various wavelength ranges. 

Open with DEXTER  
In the text 