New Hydrodynamic Solutions for Line-driven Winds of Hot Massive Stars Using the Lambert W-function

Taking into account the differentiation between radiation and line-acceleration outlined in the previous paragraphs, we rewrite the equation of momentum as

$$\begin{eqnarray}&&{\unicode{x0028B}}\displaystyle \frac{{d}{\unicode{x0028B}}}{{dr}}=-\displaystyle \frac{1}{\rho }\displaystyle \frac{{dp}}{{dr}}-\displaystyle \frac{{{GM}}_{\mathrm{eff}}}{{r}^{2}}+{g}_{\mathrm{line}},\end{eqnarray} \tag{ 7 }$$

where p is the gas pressure, ʋ is the wind velocity, and GM_eff/r² is the effective gravitational acceleration introduced in Equation (5).

Following Müller & Vink (2008) and Araya et al. (2014), we express Equation (7) in dimensionless form by introducing the dimensionless variables

$$\begin{eqnarray}&&\hat{r}=\displaystyle \frac{r}{{R}_{* }},\ \ \ \hat{{\unicode{x0028B}}}=\displaystyle \frac{{\unicode{x0028B}}}{a},\end{eqnarray} \tag{ 8 }$$

where R_* is the stellar radius and a is the isothermal sound speed given by

$$\begin{eqnarray}&&{a}^{2}=\displaystyle \frac{{k}_{B}{T}_{\mathrm{eff}}}{\mu \,{m}_{{\rm{H}}}}+\displaystyle \frac{1}{2}{{\unicode{x0028B}}}_{\mathrm{turb}}^{2}.\end{eqnarray} \tag{ 9 }$$

Unlike Sander et al. (2017), we need to assume a constant sound speed. However, the assumption of an isothermal wind near the sonic point is generally an excellent approximation for O stars, and above the sonic point, the sound speed quickly becomes irrelevant for the dynamics. The turbulent velocity ʋ_turb is included to solve the e.o.m. and is commonly assumed to be on the order of ∼10–15 km s⁻¹ (Villamariz & Herrero 2000; Puls et al. 2005). In principle, the turbulent velocity can be derived by spectral fitting.

Defining the dimensionless line-acceleration

$$\begin{eqnarray}&&{\hat{g}}_{\mathrm{line}}=\displaystyle \frac{{R}_{* }}{{a}^{2}}{g}_{\mathrm{line}},\end{eqnarray} \tag{ 10 }$$

the e.o.m. reads

$$\begin{eqnarray}&&\hat{{\unicode{x0028B}}}\displaystyle \frac{d\hat{{\unicode{x0028B}}}}{d\hat{r}}=-\displaystyle \frac{1}{\rho }\displaystyle \frac{{dp}}{d\hat{r}}-\displaystyle \frac{{\hat{{\unicode{x0028B}}}}_{\mathrm{crit}}^{2}}{{\hat{r}}^{2}}+{\hat{g}}_{\mathrm{line}},\end{eqnarray} \tag{ 11 }$$

with ʋ_crit being the rotational breakup velocity in the case of a corresponding rotating star (Araya et al. 2014)

$$\begin{eqnarray}&&{\hat{{\unicode{x0028B}}}}_{\mathrm{crit}}=\displaystyle \frac{{{\unicode{x0028B}}}_{\mathrm{esc}}}{a\sqrt{2}}=\displaystyle \frac{1}{a}\sqrt{\displaystyle \frac{{{GM}}_{\mathrm{eff}}}{{R}_{* }}}.\end{eqnarray} \tag{ 12 }$$

With the use of the equation of state for an ideal gas (p = a² ρ), Equation (11) becomes

$$\begin{eqnarray}&&\left(\hat{{\unicode{x0028B}}}-\displaystyle \frac{1}{\hat{{\unicode{x0028B}}}}\right)\displaystyle \frac{d\hat{{\unicode{x0028B}}}}{d\hat{r}}=-\displaystyle \frac{{\hat{{\unicode{x0028B}}}}_{\mathrm{crit}}^{2}}{{\hat{r}}^{2}}+\displaystyle \frac{2}{\hat{r}}+{\hat{g}}_{\mathrm{line}}.\end{eqnarray} \tag{ 13 }$$

and is formally independent of density (and therefore independent of mass-loss rate).

This expression has the advantage that if the line-acceleration is expressed as a function of the radius only ${g}_{\mathrm{line}}(\hat{r})$, the velocity and radial terms are separated on either side of the equation. Therefore, by means of the Lambert W-function (see the mathematical definition and its main properties in Appendix A), Equation (13) can be solved by integrating both sides along the atmosphere and expressing $\hat{{\unicode{x0028B}}}(\hat{r})$ in terms of the defined W(z) (see Equation (A2)).

In Equation (13) the explicit dependence on the density ρ(r) has vanished and as a consequence the Lambert W-function cannot provide a new value for the mass-loss rate. Thus, $\dot{M}$ can be considered as a free parameter of the e.o.m. Our solution to the e.o.m. will be consistent with the assumed line force but this will (in general) be inconsistent with the new line force computed by a subsequent CMFGEN iteration. As discussed in more detail in Section 3.2, a fully consistent solution of the mass-loss rate will require an iterative adjustment for the mass-loss rate. However, given the uncertainties and assumptions, it can be desirable to leave the mass-loss rate fixed and simply iterate until our solution for the velocity law stabilizes (Section 4, Appendix C).

2.1.1. Subsonic versus Supersonic Region

Assuming a monotonic behavior of the velocity field throughout the atmosphere ($d\hat{{\unicode{x0028B}}}/d\hat{r}\gt 0$), it is seen that the left-hand side (hereafter lhs) of Equation (13) becomes zero when ʋ = a because then $\hat{{\unicode{x0028B}}}=1$. Thus, at the sonic (or critical) point ${\hat{r}}_{c}$, the right-hand side (hereafter rhs) must also be zero, i.e.,

$$\begin{eqnarray}&&-\displaystyle \frac{{\hat{{\unicode{x0028B}}}}_{\mathrm{crit}}^{2}}{{\hat{r}}_{c}^{2}}+\displaystyle \frac{2}{{\hat{r}}_{c}}+{\hat{g}}_{\mathrm{line}}({\hat{r}}_{c})=0.\end{eqnarray} \tag{ 14 }$$

Because we are assuming that the line force is only a function of r, the formal critical point of the e.o.m. is the sonic point, which was also assumed by Lucy & Solomon (1970). However, using the Sobolev theory, the line force depends on dv/dr, and this moves the critical point above the sonic point. The requirement that the wind moves smoothly through the critical point sets the velocity gradient at the critical point and (potentially) yields a unique solution for the mass-loss rate. In more recent work, Lucy (2007) assumed that radiative acceleration in the neighborhood of the sonic point was a function of ʋ, and this also leads to the critical point being at the sonic point.

Due to the monotonic behavior of ʋ(r) (and consequently $\hat{{\unicode{x0028B}}}$), the sonic point becomes a boundary between two regions. The first of them is for $\hat{\unicode{x0028B}}\lt 1$, which makes both sides of Equation (13) negative. It is called the subsonic region. The second, for $\hat{{\unicode{x0028B}}}\gt 1$, is then called the supersonic region. The difference between them is not only due to mathematical analysis but also to physical reasons: in a one-dimensional fluid moving at velocities below sound speed, perturbations are propagated both inwards and outwards, whereas perturbations occurring in a supersonic fluid are only propagated downstream.

This means that we actually are searching solutions for the equation of momentum in two different branches that merge at the critical point: the analytical solution for the velocity profile implies the integration from the sonic point both outwards and inwards.

In terms of the Lambert function, the general solution is

$$\begin{eqnarray}\begin{array}{rcl}{\hat{{\unicode{x0028B}}}}^{2}{e}^{-{\hat{{\unicode{x0028B}}}}^{2}} & = & {\left(\displaystyle \frac{{\hat{r}}_{c}}{\hat{r}}\right)}^{4}\exp \left[-1-2\,{\hat{{\unicode{x0028B}}}}_{\mathrm{crit}}^{2}\left(\displaystyle \frac{1}{\hat{r}}-\displaystyle \frac{1}{{\hat{r}}_{c}}\right)\right]\\ & & \times \exp \left[-2{\displaystyle \int }_{{\hat{r}}_{c}}^{\hat{r}}{\hat{g}}_{\mathrm{line}}\,d\hat{r}\right],\end{array}\end{eqnarray} \tag{ 15 }$$

See Appendix B for a detailed step-by-step of the integration approach.

The radiative acceleration, as a function of wind location, is obtained from a converged CMFGEN model. It is evaluated using the formula:

$$\begin{eqnarray}&&{g}_{\mathrm{rad}}(r)=\displaystyle \frac{{\chi }_{\nu }{L}_{\nu }}{4\pi c\rho {r}^{2}}\ \ =\displaystyle \frac{{\chi }_{{\rm{f}}}{L}_{* }}{4\pi c\rho {r}^{2}},\end{eqnarray} \tag{ 16 }$$

where χ_f is the flux mean opacity

CMFGEN provides the radiative acceleration as a function of radius, and this can be used directly with our Lambert-procedure to obtain a new velocity law. However, it is widely known that the radiative acceleration is also a function of the velocity and velocity gradient, and hence, in general, the radiative acceleration in CMFGEN will change when we use the new velocity law. Thus, an iterative procedure is adopted.

The dimensionless equation of momentum (Equation (13)) is not explicitly dependent on density because ρ(r) terms are mathematically canceled. As a consequence, the Lambert W-function cannot calculate a new mass-loss rate and hence $\dot{M}$ here is considered as a free parameter. However, the dependence on density, and therefore on mass-loss rate, is implicitly included in Equation (13) via the g_line term. From Equation (16), the line-acceleration appears to be inversely proportional to the mass-loss rate $\dot{M}$, and hence also inversely proportional to the density ρ(r). However, there is also a further hidden dependence arising from line driving. For O-type stellar winds, Sobolev theory shows that g_line scales as ∼${\dot{M}}^{-\alpha }$, where α is the CAK parameter and is around 2/3 for O-star winds (Puls et al. 2008). Because the line force is density dependent, different mass-loss rates will lead to different final Lambert hydrodynamics. Higher mass-loss rates produce slower line accelerations (see Figure 1), which yield velocity profiles with lower terminal velocities.

Zoom In Zoom Out Reset image size

Line-acceleration is also affected by small-scale inhomogeneities (i.e., clumping) present throughout the wind. Clumping alters the line force and modifies spectral features produced in the wind. Clumping is implemented in CMFGEN in terms of the volume-filling factor f. CMFGEN assumes a void interclump medium and that the clumps are much smaller than the mean free path of a photon. The density inside the clumps, ρ_cl, satisfies ρ_cl = ρ₀/f, where ρ₀ is the homogeneous (smooth) wind density at the same radius (Bouret et al. 2005). The volume-filling factor is defined in terms of the velocity field as

$$\begin{eqnarray}&&f({\unicode{x0028B}}(r))={f}_{\infty }+(1-{f}_{\infty }){e}^{-{\unicode{x0028B}}(r)/{{\unicode{x0028B}}}_{\mathrm{cl}}}\end{eqnarray} \tag{ 17 }$$

although other options are also available. For some of the simulations presented in this paper, we used a slightly modified form that forced f to unity at the sonic point. f_∞ can be interpreted as the volume-filling factor at a larger distance, while ʋ_cl indicates the onset velocity of clumping and can be interpreted as the point in the velocity field where clumping effects start.

When clumping is included in stellar winds, inconsistencies will potentially arise in the solution of the steady-state hydrodynamic equations, which assume a smooth flow. With the inclusion of clumping in CMFGEN, the mass conservation equation becomes

$$\begin{eqnarray}&&{\rho }_{\mathrm{cl}}=\displaystyle \frac{\dot{M}}{4\pi {r}^{2}{\unicode{x0028B}}(r)f},\end{eqnarray} \tag{ 18 }$$

and the e.o.m. is

$$\begin{eqnarray}&&\left({\unicode{x0028B}}-\displaystyle \frac{{a}^{2}}{{\unicode{x0028B}}}\right)\displaystyle \frac{{d}{\unicode{x0028B}}}{{dr}}={g}_{\mathrm{rad}}-\displaystyle \frac{{{GM}}_{* }}{{r}^{2}}+2\displaystyle \frac{{a}^{2}}{r}+\displaystyle \frac{{a}^{2}}{f}\displaystyle \frac{{df}}{{dr}}.\end{eqnarray} \tag{ 19 }$$

This last expression is a more general articulation from Equation (2), using the clumped density for the equation of state of the ideal gas, p = a² f ρ_cl, and relaxing the isothermal condition previously imposed. ⁸ Similar to the case of Equation (13), we ignore the radial derivative over the sound speed term a²(r) (Equation (9)) because it has been empirically demonstrated that its value is almost constant over the wind (see Figure 1 from Sander et al. 2017).

If we consider a homogeneous wind, the last term of the rhs of Equation (19) vanishes and we can recover the original Equation (13). The term will also become unimportant above the sonic point where gas pressure forces become negligible. Unfortunately, studies of O-star spectra suggest that clumping effects extend down to the photosphere (Bouret et al. 2003; Hillier et al. 2003; Puls et al. 2006). However, a more important issue is that the hydrodynamic equation as written is inconsistent and takes no account of the interclump medium and clump formation.

Above the sonic point, a stronger clumping factor, expressed by a smaller filling factor f_∞, gives a stronger line-acceleration (see Figure 2), which, for a given mass-loss rate, will lead to higher terminal velocities. Clumping boosts line-acceleration because clumping enhances recombination, which enhances line-acceleration because the atomic levels producing the lines responsible for driving the flow will have a larger population. This effect can also be seen in m-CAK theory through the δ parameter, which indicates that the line force has an additional dependence proportional to ${({N}_{{\rm{e}}}/W(r))}^{\delta }$, where N_e is the electron density and W(r) the dilution factor ⁹ (Abbott 1982).

Zoom In Zoom Out Reset image size

As readily apparent, and which has been noted many times previously (e.g., Lucy & Solomon 1970), Equation (19) implies that the radiative force is equal to the gravitational force at the sonic point, i.e.,

$$\begin{eqnarray}&&{g}_{\mathrm{grav}}-{g}_{\mathrm{rad}}=0\ \ \ \mathrm{when}\ \ \ {\unicode{x0028B}}=a\end{eqnarray} \tag{ 20 }$$

(see also Lucy 2007). ¹⁰ This singularity condition is also introduced in the self-consistent studies of Sander et al. (2017) and Sundqvist et al. (2019) by means of the new factor Γ = g_rad/g_grav, which must be 1 at the sonic point. ¹¹ It becomes evident then that the initial value given in our CMFGEN models for mass-loss rate and for clumping will alter the equilibrium of Equation (19), especially around the sonic point singularity, and then the derived hydrodynamic solution.

Following this, the accuracy of the mass-loss rate could be checked using the singularity criterion (Equation (20)) and the hydrodynamic error determined from the CMFGEN modeling. Typically for CMFGEN models this is defined by

$$\begin{eqnarray}&&\mathrm{Err}(r)=\displaystyle \frac{200\times \left({\unicode{x0028B}}\tfrac{{d}{\unicode{x0028B}}}{{dr}}+\tfrac{1}{\rho }\tfrac{{dP}}{{dr}}+\tfrac{P}{\rho }\tfrac{{dT}}{{dr}}-{g}_{\mathrm{tot}}\right)}{\left|{\unicode{x0028B}}\tfrac{{d}{\unicode{x0028B}}}{{dr}}\right|+\left|\tfrac{1}{\rho }\tfrac{{dP}}{{dr}}\right|+\left|\tfrac{P}{\rho }\tfrac{{dT}}{{dr}}\right|+| {g}_{\mathrm{tot}}| }.\end{eqnarray} \tag{ 21 }$$

Other error prescriptions are also possible. In the above, for example, the use of ∣g_tot∣ rather than ∣g∣ + ∣g_r ∣ exacerbates the error at the sonic point (as ∣g_tot∣ = 0 at the sonic point). Another alternative is to express the error relative to the radiative force.

The two criteria discussed above can be used to constrain the mass-loss rate to be used in the Lambert-procedure. In principle we can adjust $\dot{M}$ to look for the value that minimizes the error. However, the hydrodynamic solution is influenced by several other (uncertain) free parameters besides the mass-loss rate (e.g., ʋ_turb, X-rays, and parameters associated with clumping), and the neglect of other factors that will influence the hydrodynamics (e.g., rotation). In practice we have adjusted some of these parameters to improve the hydrodynamics, but given the uncertainties we did not attempt to get perfect agreement with the hydrodynamics at all velocities. Importantly, we wanted to retain some consistency of the parameters, particularly the mass-loss rate, with values obtained from empirical model fitting.

Calculation of the new velocity profile is done by solving Equation (15), using for ${\hat{g}}_{\mathrm{line}}(\hat{r})$ the output of a CMFGEN model. The formal expression for ʋ(r) comes from Equation (B6) (see Appendix B).

From the two real branches of the Lambert W-function (see Appendix A), the one with asymptotic behavior is W₋₁. We use this branch to obtain the velocity profile in the supersonic region, i.e., from the sonic point outwards. As a consequence, each new step in the Lambert-procedure will have a new terminal velocity and a new exponential behavior different from the initial β-law.

To obtain the solution in the subsonic region, the W₀ branch of the Lambert W-function can be used. However, we decided not to do so, as large errors were produced for the acceleration and velocity profiles in the lowest part of the wind (in the hydrostatic region) when using the CMFGEN solution. Rather we used a modification of ʋ(r) below the sonic point by rescaling the velocity profile with respect to the initial model to ensure a smooth connection to the supersonic solution.

4.1.1. Junction of Subsonic and Supersonic Regions

After determining the critical point r_c from Equation (14) and calculating the new Lambert velocity profile, we spline both regions (new supersonic and old subsonic) at the critical point, where ʋ(r_c ) = a. This coupling must ensure a smooth transition between subsonic and supersonic zones. The conditions

$$\begin{eqnarray}&&\mathop{\mathrm{lim}}\limits_{r\to {r}_{{\rm{c}}}-}{\unicode{x0028B}}(r)=\mathop{\mathrm{lim}}\limits_{r\to {r}_{{\rm{c}}}+}{\unicode{x0028B}}(r)\end{eqnarray} \tag{ 22 }$$

and

$$\begin{eqnarray}&&\mathop{\mathrm{lim}}\limits_{r\to {r}_{{\rm{c}}}}\frac{{d}{\unicode{x0028B}}(r)}{{dr}}{\left|r\lt {r}_{{\rm{c}}}=\mathop{\mathrm{lim}}\limits_{r\to {r}_{{\rm{c}}}}\frac{{d}{\unicode{x0028B}}(r)}{{dr}}\right|}_{r\gt {r}_{{\rm{c}}}}\end{eqnarray} \tag{ 23 }$$

must be fulfilled.

To be in agreement with the CMFGEN units system (where the derivative of the velocity profile is expressed by the term $\sigma \equiv d\mathrm{ln}{\unicode{x0028B}}/d\mathrm{ln}r+1$), the second limit is easier to evaluate if we use logarithm derivatives, which are related by

$$\begin{eqnarray}&&\displaystyle \frac{d\mathrm{ln}{\unicode{x0028B}}}{d\mathrm{ln}r}=\displaystyle \frac{r}{{\unicode{x0028B}}(r)}\displaystyle \frac{{d}{\unicode{x0028B}}}{{dr}}.\end{eqnarray} \tag{ 24 }$$

Then,

$$\begin{eqnarray}&&\mathop{\mathrm{lim}}\limits_{r\to {r}_{{\rm{c}}}}\displaystyle \frac{d\mathrm{ln}{\unicode{x0028B}}(r)}{d\mathrm{ln}r}{\left|r\lt {r}_{{\rm{c}}}=\mathop{\mathrm{lim}}\limits_{r\to {r}_{{\rm{c}}}}\displaystyle \frac{d\mathrm{ln}{\unicode{x0028B}}(r)}{d\mathrm{ln}r}\right|}_{r\gt {r}_{{\rm{c}}}}\ .\end{eqnarray} \tag{ 25 }$$

As a consequence of these conditions, the subsonic region must be readapted, too, to have a smooth transition and avoid instabilities. Because our r_c is determined from the singularity of the e.o.m., this does not even coincide with the old velocity profile at ʋ = a (differences on these two sonic points were also previously referred by Sander et al. 2017, see their Figure 3). Figure 4 shows the disconnection of ʋ(r) and $d\mathrm{ln}{\unicode{x0028B}}/d\mathrm{ln}r$ around the sonic point for a velocity profile after solving the e.o.m. via the Lambert W-function. Therefore, we proceed to rescale the radius in the next model to match ʋ(r_c ) = a. The value of the rescaling depends on where in the old ʋ(r) we get ʋ = a from both sides as in Equation (23). Because it can happen that the differences in the velocity gradients between the subsonic and supersonic regions are too high, the spline can be also done below or over the critical point (up to 3 ND depth points depending of the case). In any case, the value of the rescaling is typically around 1.001 or 1.003 (i.e., the radius is never modified by more than 0.3%).

Zoom In Zoom Out Reset image size

Given the line-acceleration as an input, the Lambert W-function provides a new velocity profile, whereas CMFGEN provides a new line-acceleration given a prescribed velocity profile. Therefore, a self-consistent solution requires the combination of both these calculations in an iterative loop until reaching some convergence criteria.

The convergence of the Lambert iterations is checked by evaluating both the velocity field ʋ(r) and the line-acceleration g_line(r). To accept a Lambert model as well converged, we require that the conditions

$$\begin{eqnarray}&&\left|\mathrm{log}\displaystyle \frac{{\unicode{x0028B}}{(r)}_{n}}{{\unicode{x0028B}}{(r)}_{n-1}}\right|\leqslant 0.01,\end{eqnarray} \tag{ 26 }$$

and

$$\begin{eqnarray}&&\left|\mathrm{log}\displaystyle \frac{{g}_{\mathrm{line}}{(r)}_{n}}{{g}_{\mathrm{line}}{(r)}_{n-1}}\right|\leqslant 0.01,\end{eqnarray} \tag{ 27 }$$

are fulfilled for all r. Here, n denotes the nth Lambert iteration executed. ¹²

The definition of convergence we apply in this section only relates to the execution of the Lambert-procedure to solve the dimensionless e.o.m. We assume then that conditions Equations (26) and (27) imply that Equation (13) is satisfied thought the full wind. Unfortunately, we cannot extend this implication to the full e.o.m. introduced in Equation (19), because the Lambert W-function is applied over the simplified version. As a consequence, the convergence of a Lambert-procedure does not assure the fulfillment of the e.o.m. inside CMFGEN. Given the explicit independence on density, we can mathematically find a wide family of well-converged solutions for velocity profiles from different sets of mass-loss rates and clumping factors, independent of whether only one specific set of $\dot{M}$ and f_∞ would fully satisfy Equation (19).

The thresholds from Equations (26) and (27) are applied over all the range of the radial coordinate r, from the photosphere outwards and with special consideration around the sonic point. Because this is the zone where the solution coupling is done, it is also the zone with more instability at the moment of the convergence.

Besides purporting to be a solution for the velocity profile, these solutions under the Lambert W-function also must be proven to be unique. Wind hydrodynamics expressed under the β-law depends on two parameters: terminal velocity ʋ_∞ and the β exponent. If the Lambert-procedure is able to generate a well-converged hydrodynamic solution, it should be the same independently of any initial velocity profiles (i.e., initial ʋ_∞ and β parameters) chosen.

Typically, the Lambert-procedure converges after five or six Lambert iterations, each one of them corresponding to a CMFGEN model executed with about 12 inner iterations. In some cases the Lambert-procedure may not converge, and in such cases the mass-loss rate is adjusted to facilitate convergence.

Comparisons between the initial and converged self-consistent velocity fields are shown in Figure 5. Both velocity profile and line-acceleration were evaluated to satisfy the threshold conditions from Equations (26) and (27). Because the new model is self-consistent, both the lhs and rhs of Equation (2) are in agreement as shown in Figure 6 (see also Figure 20).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Unicity of the Lambert solution is verified by running CMFGEN models with different initial values for β and ʋ_∞. Figure 7 shows the unique convergence of all these parallel models: the left panel represents the initial velocity profiles with the β-law, whereas the right panel shows the final converged Lambert velocity profiles (where the overlap demonstrates all the initial models converging into the same solution).

Zoom In Zoom Out Reset image size

Because the condition Γ = 1 was already satisfied in the initial model, we did not modify the mass-loss rate nor the clumping factor. The run of Γ as a function of normalized velocity is shown in Figure 8.

Zoom In Zoom Out Reset image size

With the converged self-consistent hydrodynamics, the resulting terminal velocity has increased by a factor of ∼1.2 from 2300 to 2740 km s⁻¹. Besides, due to the rescaling in the subsonic region, there are no significant differences in the resulting velocity profile below the sonic point, just above ∼12–13 km s⁻¹.

Moreover, from Figure 9, we observe that it is possible to find a β value that can closely fit the new velocity field. This demonstrates that velocity profiles can be approximated by a β-law at large scale only, but around the critical point it would be necessary to use an alternative prescription to recover the shape of ʋ(r), the reason why we defined this new exponent as "quasi-β"). Several such procedures are implemented in CMFGEN, one of which is very similar to that proposed by Björklund et al. (2021).

Zoom In Zoom Out Reset image size

A good solution for HD 163758 was found by decreasing the mass-loss rate compared with the initial value provided by Bouret et al. (2012) and tabulated in Table 1. The run of Γ as a function of normalized velocity is shown in Figure 10. We observe, in this case, that our initial model for HD 163758 did not satisfy the criteria Γ = 1 at the sonic point, even after adjustments of the velocity law in the neighborhood of the sonic point. Thus, we reduced the mass-loss rate from 1.5 × 10⁻⁶ to 1.2 × 10⁻⁶ M_⊙ yr⁻¹.

Zoom In Zoom Out Reset image size

Checks for the unicity of the solution were the same as for ζ Puppis. Comparisons between the old and the self-consistent new velocity profile are shown in Figure 11. The increment in the terminal velocity for a Lambert model is 1.14 times, from 2100 to 2400 km s⁻¹, lower than that for ζ Puppis. However, we observe the "departure" of the Lambert velocity profile from the initial ʋ(r) at a much deeper zone, near ʋ ∼ 2 km s⁻¹. This occurs because the new Lambert velocity profile has a lower value for its quasi-β (0.85) compared with the initial introduced β parameter (1.1).

Zoom In Zoom Out Reset image size

Despite these encouraging result, the instability of the model could not be reduced to the same level as for ζ Puppis. As can be seen from Figure 12, the agreement for the e.o.m. this time is not satisfied as well as in the case for ζ Puppis.

Zoom In Zoom Out Reset image size

Particularly for this star, we tested the effects of the modification of the turbulent velocity ʋ_turb. We found that its reduction from 10 to 7 km s⁻¹ reduced the error of the CMFGEN model, without any significant variation over the resulting velocity profile.

Alternatively, we have found that a significant decrease of the mass-loss rate led to a more stable converged Lambert model and therefore a better agreement of Equation (2). This is observed in Figure 13, where both sides exhibit solid agreement. The derived $\dot{M}=8.3\times {10}^{-7}{M}_{\odot }\,{\mathrm{yr}}^{-1}$ is almost a factor of 2 lower than the initial model, while ʋ_∞ = 3100 km s⁻¹ is factor of 1.5 times higher. Such strong differences affect the determination of the Γ. Figure 14 shows that the condition Γ = 1 is not satisfied for the sonic point now, which can be assumed as a reinforcement that the alternative mass-loss rate is too low. This implies that minimization of the CMFGEN error may not lead to a smooth solution through the sonic point.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

A solution for α Cam was found for a higher mass-loss rate and an increased volume-filling factor (i.e., decreased clumping) compared to the initial values displayed in Table 1 that are similar to the ones provided by Najarro et al. (2011). The run of Γ as a function of normalized velocity is shown in Figure 15 where, as for HD 163758, Γ = 1 at the sonic point was not satisfied in the initial model.

Zoom In Zoom Out Reset image size

Comparisons between the old and the new self-consistent velocity profile are shown in Figure 16. This star presents the highest increment in the terminal velocity for a Lambert model among the stars studied, from 1550 to 3420 km s⁻¹, a factor of ∼2.2. The differences in the velocity profiles (for both the Lambert and initial models) go much deeper into the subsonic region (ʋ ≲ 1 km s⁻¹), which again can be attributed to the large departure of the resulting quasi-β (0.95) from the initial value of the β parameter (1.5).

Zoom In Zoom Out Reset image size

This great enhancement of the velocity profile also causes instabilities in the hydrodynamic structure of the model, making it impossible to obtain an accurate agreement for the equation of momentum, as illustrated in Figure 17.

Zoom In Zoom Out Reset image size

For α Cam, we tested the effects of the modification of the ʋ_cl (see Equation (17)). The modification from 37.5 to 100 km s⁻¹ reduced the error of the CMFGEN model but led to a Lambert solution with a terminal velocity increased by ∼500 km s⁻¹. However, the increase in mass-loss rate to 8.5 ×10⁻⁷ M_⊙ yr⁻¹ (a factor of ∼1.2 higher than the previous Lambert model) led to a more stable model and a lower terminal velocity (2650 km s⁻¹). Figure 18 shows a more precise agreement for both sides of Equation (2), and Figure 19 shows that the condition Γ = 1 for the sonic point is also satisfied. Thus, the alternative model satisfies better our two criteria and therefore provides a more reliable value for the mass-loss rate.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

An example of the discrepancy between observation and theory is presented in Section 6.5 of the study of Bouret et al. (2012), where they checked the consistency of the line force for ζ Puppis. A truly self-consistent solution must satisfy the stationary equation of momentum (Equation (2)). Both the lhs and rhs of this equation are plotted in Figure 20. This g_rad(r) was calculated internally by CMFGEN starting from the initial parameters outlined in Table 1. However, even when the choice of these parameters achieved a better agreement between the lhs and rhs for Equation (2) compared with those stellar and wind parameters determined by the spectral fit (see the discussion in Bouret et al. 2012), both sides are still far from being in agreement. Hence, this particular CMFGEN model is not hydrodynamically self-consistent.

The discrepancy shown in Figure 20 arises because the velocity field is calculated from a β-law rather than from the line acceleration. It could be argued that this discrepancy might be only in a specific star, but this was previously also noticed by Puebla et al. (2016). Hence, the lack of consistency between hydrodynamics and radiative acceleration for a CMFGEN model with a velocity profile fixed by a β-law seems to be a general rule, and we can infer the same inconsistency for other radiative transfer codes (e.g., FASTWIND, PoWR) if they are also using the β-law. On the contrary, Lambert solutions present a remarkably good agreement between both sides of Equation (2), especially for those Lambert models with the smallest error. In Figure 21, we compared the resulting Lambert model for ζ Puppis with our other two stars HD 163758 and α Cam. Figure 21 is in concordance with that shown in Figures 6, 12, 13, 17, and 18. The lower the error, the better the agreement between the lhs and rhs of the e.o.m. for the CMFGEN model. This provides evidence for the validity of the velocity profiles computed from the Lambert W-function (instead if using the β-law) and also that the e.o.m. was properly satisfied for the CMFGEN model (Equation (19)), despite the condition Γ = 1 not being fulfilled at ʋ = a. In such a case the critical point for the equation is not the sonic point but in its close neighborhood.

The search for a fully self-consistent solution (coupling line acceleration, hydrodynamics, and radiative transfer) for the wind of massive stars has been approached previously by other studies (Puls et al. 2000; Kudritzki 2002; Pauldrach et al. 2012; Krtička & Kubát 2017; Sander et al. 2017). In particular, we emphasize the work done by Sander et al. (2017), where another iterative procedure was implemented in order to obtain a self-consistent solution for the wind hydrodynamics of ζ Puppis under the non-LTE regime, using the radiative transfer code PoWR (Gräfener et al. 2002; Hamann & Gräfener 2003). In that study, radiative acceleration was also obtained from the output of the radiative transfer solution (PoWR for them, CMFGEN for us) and that acceleration was reused to obtain a new hydrodynamic profile.

Although their study shares a similar philosophy to ours, there are important differences. First, our new velocity profile is calculated by means of the Lambert W-function, whereas their ʋ(r) is recalculated by updating the stratification of the wind. The advantage of the Lambert-procedure is that it ensures the existence of a unique solution, i.e., the final converged Lambert model does not depend on the initial velocity profile. Besides, our Lambert-procedure considers iterative changes in the velocity profile only, letting the stellar parameters (temperature, radius, mass, etc.) free and the other wind parameters (mass-loss rate and clumping factor) can be simultaneously modified to reduce the error in the models. Such relaxation imposed by us allows the search by eye inspection of the best spectral fit, reducing the number of free parameters (hydrodynamic will depend now on the initial stellar parameters and the constrained $\dot{M}$ and f_∞). Hence, we gain the flexibility to constrain the mass-loss rate either by observations or by theoretical criteria, keeping in every case the self-consistency for the velocity profile.

Because of the different numerical approaches adopted in the literature, we select these two previous studies to compare with the results given by our Lambert-procedure. The comparison among these three self-consistent solutions for ζ Puppis is summarized in Table 2. In addition, we could have included the study made by Pauldrach et al. (2012), who also analyzed extensively ζ Puppis by deriving self-consistent values for the stellar and wind parameters of the star. However, they assumed a different value for the stellar radius of ∼28 R_⊙, far from the ∼18 R_⊙ adopted by Sander et al. (2017), Gormaz-Matamala et al. (2019), and the present work, and hence we decided not to include it in the comparisons.

Table 2. Summary of Self-consistent Solutions Performed by Sander et al. (2017), Gormaz-Matamala et al. (2019), and This Present Study (Table 1)

	Sander et al. (2017)	Gormaz-Matamala et al. (2019)	Björklund et al. (2021)	This work
RT code	PoWR	FASTWIND	FASTWIND	CMFGEN
Hydro method	⋯	HydWind	⋯	Lambert W-function
Teff [kK]	40.7	40	40.7	41
$\mathrm{log}g$	3.63	3.64	3.65	3.6
R*/R⊙	15.9	18.7	18.91	17.9
v∞ [km s−1]	2046	2700	2302	2740
$\dot{M}$ [10−6 M⊙ yr−1]	1.6	5.2	2.5	2.7
f∞ [1/D∞]	0.1	0.2	⋯	0.10

Note. These three models assume the same abundances as Bouret et al. (2012). The radiative transfer code used to perform the respective synthetic spectra is also marked, together with the code/methodology used to calculate the hydrodynamics. Equivalence between filling factor and clumping factor D_∞ (used by PoWR) was made assuming the interclump medium is void, following Sundqvist & Puls (2018). Additionally, we include from the grid of Björklund et al. (2021), the closest parameters to ζ Puppis.

Download table as:  ASCIITypeset image

As an initial comment, we remark on the discrepancy between the wind parameters derived by Sander et al. (2017) and those derived by us, even when differences between stellar parameters are small. First, for the case of the terminal velocity, Sander's value of 2046 km s⁻¹ lies below typical values obtained by spectral fitting: 2250 km s⁻¹ by Puls et al. (2006) and 2300 km s⁻¹ by Bouret et al. (2012). On the contrary, the final terminal velocity provided by the Lambert-procedure is ∼20% higher than the "observed" ʋ_∞; this is seen, for example, in the fit for C iv λ λ 1548, 1551. However, their terminal velocity is close to that determined by Paper I's m-CAK prescription (see Table 5 from Gormaz-Matamala et al. 2019).

The difference in methodology between the current paper, which uses the Lambert W-function, and Paper I, which uses m-CAK theory, can be summarized as follows. Paper I:

• 1.  
  used a simplified "quasi-NLTE" scenario for the treatment of atomic populations, following formulations performed by Mazzali & Lucy (1993) and Puls et al. (2000) whereas in this work we solve the proper statistical equilibrium equations when CMFGEN is executed.
• 2.  
  used a flux field calculated by Tlusty, which uses the plane-parallel approximation whereas our flux field is calculated with CMFGEN. The radiation field computed in this case includes a full treatment of the line blanketing and the multiline scattering in spherical geometry (Puls 1987).
• 3.  
  did not consider effects from clumping upon the resulting g_line, whereas, from Figure 2, it is clear that line-acceleration changes when the mass-loss rate stays constant but the clumping factor is modified.
• 4.  
  creates hydrodynamic models starting from the photosphere of the star using HydWind, whereas with the Lambert-procedure velocity profiles are calculated from the sonic point outwards.

Therefore, these differences between the codes could easily explain the differences in wind parameters for ζ Puppis.

Further, it is well known that rotation enhances the values for mass-loss rate on the equator of the star. The hydrodynamics calculated by HydWind in Paper I used a value for the normalized stellar angular velocity ¹⁵ of Ω = 0.39 in order to reproduce the known value of $\unicode{x0028B}\sin i=210$ km s⁻¹ for ζ Puppis (Bouret et al. 2012), whereas the Lambert-procedure does not consider rotational effects because CMFGEN assumes spherical symmetry. According to Venero et al. (2016), mass-loss rates increase their values by a factor of ∼1.2–1.3 for Ω = 0.4, which leads us to think that the self-consistent value for $\dot{M}$ obtained in Paper I would decrease down to ∼4 × 10⁻⁶ M_⊙ yr⁻¹ if Ω = 0, closer to the $\dot{M}$ determined by the Lambert-procedure.

While both m-CAK and Lambert prescriptions obtain self-consistent solutions, they work under different philosophies. In Paper I, the m-CAK prescription calculates its own self-consistent mass-loss rate, which is later tested by spectral analysis in order to check how near or far it falls from the real solution. On the other hand, the Lambert-procedure presented in this work calculates a self-consistent solution for the wind hydrodynamics with the mass-loss rate as a free parameter, which is later (or can be) constrained by the spectral fitting.

We have also included in Table 2 a comparison with the study performed by Björklund et al. (2021), who provide a grid of self-consistent solutions following the iterative procedure using FASTWIND introduced in Sundqvist et al. (2019). We selected the model with the stellar parameters closest to those assumed for ζ Puppis, and we observe that their predicted mass-loss rate is in agreement with the value set by us, although their value for terminal velocity is slightly lower than ours and better matches the "observational" value of 2300 km s⁻¹ (Bouret et al. 2012). The difference from our result may be due to their neglect of clumping. A summary of the best hydrodynamical solution for HD 162758 and α Cam found in this study, compared with the closest parameters from the grid of Björklund et al. (2021), is displayed in Tables 3 and 4, respectively.

Table 3. Summary of the Most Reliable Self-consistent Solutions for HD 163758 Performed by This Present Study (Section 5.2), Contrasted with the Closest Model from the Grid of Björklund et al. (2021)

	Björklund et al. (2021)	This Work
RT code	FASTWIND	CMFGEN
Hydro method	⋯	Lambert W-function
Teff [kK]	34.6	34.5
$\mathrm{log}g$	3.43	3.4
R*/R⊙	20.68	19.1
v∞ [km s−1]	2377	2400
$\dot{M}$ [10−6 M⊙ yr−1]	0.64	1.2
f∞ [1/D∞]	⋯	0.05

Download table as:  ASCIITypeset image

Table 4. Summary of the Most Reliable Self-consistent Solutions for α Cam Performed by This Present Study (Section 5.3), Contrasted with the Closest Model from the Grid of Björklund et al. (2021)

	Björklund et al. (2021)	This Work
RT code	FASTWIND	CMFGEN
Hydro method	⋯	Lambert W-function
Teff [kK]	28.4	28.2
$\mathrm{log}g$	3.18	2.975
R*/R⊙	23.11	30.3
v∞ [km s−1]	2972	2650
$\dot{M}$ [10−6 M⊙ yr−1]	0.18	0.85
f∞ [1/D∞]	⋯	0.05

Download table as:  ASCIITypeset image

The self-consistent velocity profiles introduced in this work have the advantage of having been analytically calculated from the e.o.m. of the wind using the Lambert-procedure. However, the e.o.m. introduced in Equation (13) assumes an isothermal and homogeneous wind. Such "simplifications" were necessary to allow the e.o.m. to be put in a form that can be solved using the Lambert W-function.

The adopted iterative procedure (Appendix C) allows a converged solution to be obtained for the wind velocity even though the full e.o.m. is not satisfied everywhere (particularly in the neighborhood of the sonic point). This was observed for HD 163758 and α Cam, where we could find two different new velocity profiles calculated by the Lambert W-function for each star and where the initial choice of mass-loss rate led to different terminal velocities.

As we have stated, different values for the mass-loss rate of each star produce different grades of errors associated with the CMFGEN model along all the velocity structures. The Lambert-procedure can provide a converged solution for the velocity profile (i.e., we can recover the previous ʋ(r) after running the Lambert W-function) for any potential set of $\dot{M}$ and f_∞ if the g_line is given, but just one set of values for the mass-loss rate and clumping factor will provide the most stable line-acceleration from the CMFGEN model. However, one of the limitations of that strategy is that the error reduction can be compared only with the intrinsic error of the previous model but does not guarantee achieving a standard threshold valid for any model (Figure 21). As a consequence, values for mass-loss rates presented in this work are those that lie close enough to the hypothetical true value, and thus we talk about constraining $\dot{M}$ instead of calculating it. Also, the CMFGEN models do not depend only on the mass-loss rate; this is particularly evident if we consider the extra features that affect the calculation of g_line.

Hydrodynamic solutions from CMFGEN models depend not only on the selected wind parameters but also on other physics, which are not commonly considered in the standard model of line-driven winds, for instance, the inclusion of X-rays, which modifies the ionization balances of the atomic populations and therefore potentially affects the lines contributing to the line acceleration. The turbulent velocity affects the line-acceleration near the sonic point because of the enhancement of the Doppler width. This influence was extensively discussed by Lucy (2007) and Sundqvist et al. (2019), where different choices of ʋ_turb produced different wind parameters. Adopting a turbulent velocity also modifies the sound speed (Equation (9)) and therefore the location of the critical point for the e.o.m. However, for HD 163758 a reduction of ʋ_turb did not produce any significant change in the final Lambert solution. In part this is expected as the biggest change will occur around the sonic point. A more detailed study of the influence of the turbulent velocity on the Lambert W-function solution is warranted.

Another important factor that influences both the resulting Lambert models and the hydrodynamic solution obtained with CMFGEN is clumping. The models presented in this paper generally characterize the clumping via two parameters, the filling factor at infinity, f_∞, and an onset velocity, ʋ_cl. Our functional form was adopted for simplicity—other forms may be more relevant, and in addition the clumping may not be a monotonic function (Puls et al. 2006; Najarro et al. 2011). The adopted value of f_∞ is crucial for the wind dynamics, whereas the influence of ʋ_cl is more uncertain. Alternative exponential functions using the Rosseland optical depth instead of the velocity field were proposed by Sander et al. (2017, see their Equation (32)), but still there is an extra free parameter to be determined. CMFGEN currently ignores porosity, but the inclusion of porosity, which is probably more important in models with a "low" mass-loss rate, could help to explain why the ʋ_∞ in hydrodynamic models is sometimes larger than that observed.

Given the numerous factors that influence the line-acceleration near the sonic point, it is extremely difficult to determine an accurate mass-loss rate from first principles. However, using the Lambert model to derive a velocity law does allow a range of mass-loss rates to be utilized in spectral fitting, and hence potentially provides tighter constraints on the derived mass-loss rates and wind properties.