2 Optically thick wind models

In this section we will derive the equations that describe the flow of the radiation and gas in the optically thick part of the wind. The wind has to be accelerated from subsonic velocities deep down in the atmosphere to supersonic velocities in the wind. This implies that the wind must pass smoothly through the sonic point where the flowspeed is equal to the local sound speed. We first describe the general equations for radiation driven winds. We then apply these to the sonic point. The requirement of a smooth transonic flow sets specific conditions to the values and gradients of several parameters, including the mass-loss rate. We will derive these conditions.

   
2.1 General equations for radiation driven winds

We assume that the winds of WR-stars can be described by a steady, spherically-symmetric flow. The fluid flow equations for radiatively driven winds are:

$$\frac {\rm d}{{\rm d}r} (\rho v r^{2}) = 0,$$ (1)

$$v \frac{{\rm d}v}{{\rm d}r} + \frac{1}{\rho} \frac{{\rm d}P_{\rm g}}{{\rm d}r} = -\frac{G M} {r^{2}} + f_{\rm rad},$$ (2)

$$v \frac{{\rm d}U_{\rm g}}{{\rm d}r} + P_{\rm g} v \frac{\rm ... ...}^{\infty} 4 \pi \chi_{\nu} ( J_{\nu} - S_{\nu} ) {\rm d} \nu.$$ (3)

These equations express conservation of respectively mass, momentum and energy of the gas. The symbols are the mass density $\rho$, the gas pressure $P_{\rm g}$, the velocity v, the (angle-) mean intensity of radiation $J_{\nu}$, the mass extinction (absortion + scattering) coefficient $\chi_{\nu}$, the radiative source function $S_{\nu}$, the internal energy $U_{\rm g}$ of the gas per unit mass, the stellar mass M and the force $f_{\rm rad}$ produced by radiation pressure. The mass-loss rate of the star is

$$\dot{M}=4 \pi \rho v r^{2}.$$ (4)

The gas pressure $P_{\rm g}$ depends on temperature T and density $\rho$as

$$P_{\rm g} = \frac{(\gamma+1)}{\mu m_{\rm u}} \rho k T,$$ (5)

where $\gamma$ is the mean number of free electrons per atom, $\mu m_{\rm u}$ is the mean atomic weight and $m_{\rm u}$ is the atomic mass unit. The internal energy $U_{\rm g}$ is related to $P_{\rm g}$ by the rule

$$U_{\rm g} = \frac{3}{2} \frac{P_{\rm g}}{\rho}\cdot$$ (6)

The isothermal sound speed is given by

$$v^{2}_{\rm s}={\rm d}P_{\rm g}/{\rm d}\rho=(\gamma+1)k T/(\mu m_{\rm u})~.$$ (7)

2.1.1 The momentum equation

Using the formulae (1), (5) and (7) we can transform the momentum Eq. (2) into the form

$$\left( v-\frac{v^{2}_{\rm s}}{v} \right) \frac{{\rm d}v}{{\r... ...ac{v_{\rm s}^{2}}{T} \frac{{\rm d}T}{{\rm d}r} + f_{\rm rad}.$$ (8)

The radiation pressure force is

$$f_{\rm rad} = \frac{4 \pi} {c} \int_{0}^{\infty} \chi_{\nu} H_{\nu} {\rm d} \nu~=~ \frac{\chi L(r) } {4 \pi r^{2} c},$$ (9)

where $H_{\nu}$ is the first moment of the radiation intensity (the Eddington flux), L(r) is the radiative luminosity in the frame comoving with the wind and $\chi$ is the flux-mean opacity (the flux-mean extinction (absorption + scattering) coefficient per unit mass).

2.1.2 The energy equation

The integrated overall energy conservation equation for both the gas and the radiation is

$$L(r) + L_{\rm adv}(r) + \mbox{$\dot{M}$ }\left\{ \frac{v^{2}... ...}{r} + \frac{P_{\rm g}}{\rho} + U_{\rm g} \right\} = \dot{E},$$ (10)

where $\dot{E}$ is the energy-loss constant (the total energy transported out per unit time across any spherical surface). The first term is the luminosity L(r) in the comoving frame. The second term is the radiative energy that is carried out (advected) by the gas flow in the comoving frame. Cassinelli & Castor (1973) have shown that the advective luminosity $L_{\rm adv}$ is

$$L_{\rm adv} = 16 {\pi}^{2} r^{2} \frac{v}{c} (J + K).$$ (11)

This term is independent of $\rho$ because the radiation density is independent of $\rho$. At large distance from the star the moments of the radiation approach $J=K=H=L(r=\infty)/(16 \pi^2 r^2)$ and so

$$\mbox{$L_{\rm adv}$ }(\infty)=\frac{2 \mbox{$v_\infty$ }L(\infty)}{c}\cdot$$ (12)

The sum of L(r) and $\mbox{$L_{\rm adv}$ }$ is the luminosity in the stationary frame (see Cassinelli & Castor 1973; Pistinner & Eichler 1995). The third term in Eq. (10) is the total gas energy (kinetic, thermal and potential) that is advected by the flow.

The energy-loss constant $\dot{E}$ can be found from general considerations. The energy rate $L_{\rm core}$ generated by the nuclear fusion must be equal to the luminosity at infinity plus the energy lost by accelerating the wind and lifting it out of the potential well (cf. Heger & Langer 1996):

$$L_{\rm core} \simeq L(\infty) \left( 1 + \frac{2 \mbox{$v_\in... ...left( \frac{v_{\infty}^{2}}{2} + \frac{GM}{R_{\rm hc}}\right),$$ (13)

where $R_{\rm hc}$ is the radius of the hydrostatic core. We see from Eqs. (10) and (12) that

$$L(\infty) \left( 1 + \frac{2 \mbox{$v_\infty$ }} {c}\right) + \dot{M} \frac{v_{\infty}^{2}}{2} \approx \dot{E}.$$ (14)

From these two equations we obtain an expression for the energy loss constant

$$\dot{E} = L_{\rm core} - \dot{M} \frac{GM}{R_{\rm hc}}\cdot$$ (15)

   
2.2 The optically thick part of the wind

The momentum Eq. (8) and the energy Eq. (10) together describe the radiation driven wind. The solution of these equations requires knowledge about the opacity and the temperature structure. For the optically thick part of the wind we can make some simplifying approximations, because the transfer of radiation can be described by diffusion. In the case of radiative diffusion the temperature gradient is

$$\frac{{\rm d} T} {{\rm d} r} = \frac{- 3 \chi_{\rm R} \rho L(r)} {64 \pi r^{2} \sigma T^{3} },$$ (16)

and the radiation-pressure force is

$$f_{\rm rad} = \frac{\chi_{\rm R} L(r) } {4 \pi r^{2} c},$$ (17)

where $\chi_{\rm R}$ is the Rosseland mean opacity. In the optically thick part of the wind the mean intensity of the radiation is $J=\sigma T^4/ \pi$ and K=J/3 (Cassinelli & Castor 1973) and so the advective luminosity (Eq. (11)) is

$$L_{\rm adv} = 16 {\pi}^{2} r^{2} \frac{v}{c} \frac{4 \sigma T^{4}} {3 \pi}\cdot$$ (18)

Defining the constants $c_{1}=64 \pi \sigma / (3 c)= 1.267 \times 10^{-13}$ and $a_{1}=k (\gamma+1) / (\mu m_{\rm u})= 8.314 \times 10^7 (\gamma+1)/\mu$ in cgs-units, where $\gamma$ is the mean number of free electrons per ion, we can write the energy equation in the optically thick part of the wind as

$$L(r) = \dot{E} - c_{1} r^{2} v T^{4}-\dot{M} \left( \frac{v^{2}}{2} - \frac{GM}{r} + \frac{5a_{1} T}{2} \right)\cdot$$ (19)

This is the form of the energy equation that we will use to find the radiative luminosity L(r) in the comoving frame that is needed for the calculation of the radiative force  $f_{\rm rad}$.

The diffusion approximation becomes valid at large optical depths (i.e., many photon mean-free-paths below the surface) where $T \gg T_{\rm eff}$ (Mihalas 1978, p. 49). For an estimate, let us adopt that if the photon mean-free-path is more than 10 times smaller than the sonic point radius then the diffusion approximation is valid. The photon mean-free-path ( $\overline{l_{\nu}}$) is defined as $1/(\rho \chi)$. Using the sonic point parameters from the paper (see Table 7), we find that in the case of the WN6 star WR136 $\overline{l_{\nu}} \approx 2.38\times 10^{9}$ cm ( $\overline{l_{\nu}}/R_{\rm s} \approx 0.0074$) and in the case of the WN5 star WR 139 $\overline{l_{\nu}} \approx 1.82 \times 10^{9}$ cm ( $\overline{l_{\nu}}/R_{\rm s} \approx 0.013$). Thus, the diffusion approximation is valid for the regions near the sonic points of WR-star winds. Note that this is not the case for the winds of O-stars, which generally have $\overline{l_{\nu}}/R_{\rm s} > 0.1$.

   
2.3 The conditions at the sonic point

The momentum Eq. (8) has the characteristic form of an equation with a critical point. The left hand side shows that when $v=v_{\rm s}$, the velocity gradient ${\rm d}v/{\rm d}r$ becomes infinite (positive or negative), unless the right hand side is zero at this sonic point. So the right hand side must vanish exactly at the critical point. The smooth passage of the flow through the sonic point also requires that the derivatives of the temperature and velocity are equal on both sides of the sonic point. These conditions are expressed by the de l'Hopital rule which describes the relations for the radial velocity gradients at the critical point (see e.g. Lamers & Cassinelli 1999, p. 421).

The momentum Eq. (8) can be written as

$$f_1~=~f_2 \frac{{\rm d}v}{{\rm d}r}$$ (20)

with

$$f_{1}= \frac{2 a_{1} T}{r}-\frac{GM}{r^{2}}-a_{1} \frac{{\rm d}T}{{\rm d}r} + \frac {\chi L(r)} {4 \pi c r^{2}}$$ (21)

and

$$f_{2}= v-\frac{a_{1} T}{v}\cdot$$ (22)

The conditions at the sonic point, $R_{\bf s}$, where the flow velocity is

$$v_{\rm s} \equiv v(R_{\rm s})=\sqrt{a_{1} T_{\rm s}}$$ (23)

with $T_{\rm s} \equiv T(R_{\rm s})$ are

$$f_1(R_{\rm s})=0~~~~~~~~{\rm and}~~~~~f_2(R_{\bf s})=0.$$ (24)

The first condition implies that

$$\chi_{\rm s} L_{\rm s} = 4 \pi c R_{\rm s}^{2} \left(-\frac{... ...+a_{1} \left(\frac{{\rm d}T}{{\rm d}r}\right)_{\rm s} \right),$$ (25)

where the subscript $\rm s$ denotes the values at the sonic point. The velocity gradient at the sonic point is given by de l'Hopital's rule

$$\left(\frac{{\rm d}v}{{\rm d}r}\right)_{\rm s}= \left(\frac{... ...rm s} / \left(\frac{{\rm d} f_2}{{\rm d}r}\right)_{\rm s}\cdot$$ (26)

Since f₂ increases from negative below the sonic point to positive above the sonic point, the gradient ${\rm d}f_1/{\rm d}r$ should be positive at the sonic point. In physical terms this means that either energy or momentum must be added to the wind as it passes through the sonic point (cf. Lamers & Cassinelli 1999, p. 100). The requirement that the velocity gradient should be continuous through the sonic point (i.e. the regularity condition) implies that

$$\left(\frac{{\rm d}^{2}v}{{\rm d}r^{2}}\right)_{\rm s}= \fra... ...rm s}} { \left(\frac{{\rm d} f_{2}}{{\rm d}r}\right)_{\rm s} }$$ (27)

with the derivatives ${\rm d}f_1/{\rm d}r$, ${\rm d}f_2/{\rm d}r$, ${\rm d}^2f_1/{\rm d}r^2$and ${\rm d}^2 f_2/{\rm d}r^2$ following from the definitions of f₁ and f₂.

We point out that the CAK-type line forces (actually the line force amplification due to Doppler shifts - Castor et al. 1975) can be neglected near the sonic point of WR-star winds. The CAK-type line force is usually expressed in relation to the radiation force by electron scattering as follows

$$f_{\rm L}=\frac{\sigma_{\rm e}^{\rm ref} M(t) L } {4 \pi c r^{2}},$$ (28)

where $\sigma_{\rm e}^{\rm {ref}}$ is the reference value for the electron scattering opacity ( $\sigma_{\rm e}^{\rm {ref}} = 0.325$ cm²g^-1, Lamers & Cassinelli 1999, p. 217). The force multiplier can be expressed as

$$M(t)=k t^{-\alpha} C,$$ (29)

where k and $\alpha$ are the force multiplier parameters, C is a correction term which is very close to unity and t is the dimensionless optical depth parameter defined as

$$t=\frac{\sigma_{\rm e}^{\rm ref} v_{\rm th} \rho}{{\rm d}v/{\rm d}r},$$ (30)

where $v_{\rm th}$ is the mean thermal velocity of protons. Our optically thick wind models show that ${\rm d}v/{\rm d}r \approx v_{\rm s}/R_{\rm s}$ near the sonic point. The value of t at the sonic point can be found from the relationship

$$t=\frac{ \sigma_{\rm e}^{\rm {ref}} \sqrt{(2 \mu/(\gamma +1))} \dot{M}} {4 \pi v_{\rm s} R_{\rm s}}\cdot$$ (31)

This shows that near the sonic point of the WN6 star WR136 $t \approx 140$ ( $R_{\rm s} \approx 4.6~ R_{\odot}$) and near the sonic point of WN5 star WR139 $t \approx 52$ ( $R_{\rm s} \approx 2.0 ~R_{\odot}$). With such large values of t and with the expected values of $\alpha$ ($\approx$0.5) and k ($\approx$0.6-0.7), the force multiplier is indeed very small, M(t)<0.05 to 0.08. This implies that the CAK-type line force due to Doppler shifts is negligible compared to electron scattering and other radiation forces near the sonic points of WR-stars. One might counterargue that the velocity gradient near the sonic point might be of order ${\rm d}v/{\rm d}r \approx v_{\infty}/R_{\rm s}$, as in optically thin (in the continuum) CAK-type O-star wind models, rather than ${\rm d}v/{\rm d}r \approx v_{\rm s}/R_{\rm s}$. If that were the case, the values of t near the sonic points of WR-stars would be smaller than given above ( $t \approx 3$ for WR136 and $t\approx 1$ for WR139), but the force multiplier is still very small 0.3 < M(t) <0.6 compared to $t \approx 10^{-2}$ and M(t)>10¹ for O-stars. In this estimate we neglected continuum absorption within the Sobolev length. This implies that the true CAK-forces near the sonic points of WR-star winds are even smaller than estimated above. So, we can safely conclude that near the sonic points of WR-star winds the CAK-type forces are not important, but they become the dominate driving force starting from some distance above the sonic point. The influence of these forces is not ignored in our models - they are taken into account indirectly by using the momentum equation in the supersonic part of the wind by assuming that the winds are radiatively driven.

   
2.4 The temperature structure

The conditions at the sonic point, described above, contain the derivatives of the temperature T(r), which depends on the velocity v(r) structure. In this section we describe these two functions. From the analytical solution of the spherically extended atmosphere in radiative equilibrium in the generalized Eddington approximation derived by Lucy (1971) we find that

$$T^{4}(r) \approx \frac{3}{4}~ \mbox{$T_{\rm eff}$ }^4(R_{\rm s})~ \Big(\tau^{\prime} + \frac{4}{3}W(r) \Bigr).$$ (32)

In this expression $\mbox{$T_{\rm eff}$ }(R_{\rm s}) \equiv (L / (4 \pi \sigma R_{\rm s}^2))^{0.25}$ where $R_{\rm s}$ is the radius of the sonic point of the wind, W(r) is the coefficient of geometrical dilution, and the effective optical depth $\tau^{\prime}$ is defined by

$$\tau^{\prime} = \int_{r}^{\infty} \chi \rho \frac {R_{\rm s}^{2}} {r^{2}} {\rm d}r,$$ (33)

where $\chi$ is the flux-mean extinction coefficient. We point out that Heger & Langer (1996) and Lucy & Abbott (1993) used a slightly different version of the original Lucy (1971) formula. They used the photospheric radius ( $R_{\rm ph}$ where $\tau^{\prime}(R_{\rm ph})=2/3$) as a reference level, whereas below the photospheric radius the geometrical dilution was neglected. Test computations have shown that our results (described below) would have been qualitatively the same if we had adopted this same reference level $R_{\rm ph}$ instead of the value of $R_{\rm s}$ that we adopted. The important argument against the use of $R_{\rm ph}$ as reference level is the large difference between $R_{\rm ph}$ and $R_{\rm s}$ for WR-stars. Therefore the effects of geometrical dilution should not be neglected in the region below $R_{\rm ph}$, where J/K is not constant. Our approach is almost in line with the arguments of the optically thick wind modeling study of the $1.4 ~M_{\odot}$ stars driven by super-Eddington luminosities by Quinn & Paczynski (1985). They stated (p. 635 in their paper) that "Models of radiatively driven winds with a critical point at a large optical depth have very diffuse photospheres with the density scale height approximately equal to the local radius. Therefore, the variation of radiation energy density with radius is due not only to diffusion through opaque matter, but also to geometrical dilution''.

We adopt the temperature structure of a spherically symmetric stationary wind in radiative equilibrium, given by Eq. (32) with the effective optical depth given by Eq. (33). This approximation is not valid at distances larger than $r \sim R(\tau \approx 2/3)$, because it assumes that $J \approx B$. For example Eq. (32) predicts that at large distance from the star, where $\mbox{$\tau^{\prime}$ }$ becomes very small, the temperature will go to zero. In reality the observations and empirical wind models show that the temperatures of WR-winds approach a constant finite value of order 10⁴ K at distances of a few tens of stellar radii. This is due to the radiative heating by photoionization. Nugis et al. (1998) have determined the asymptotic temperatures $T_{\rm asymp}$ of WR-winds as a function of WR-subtypes. We adopt these values. So the temperature structure of our models is described by Eq. (32) if $T(r) \ge T_{\rm asymp}$ and $T(r)=T_{\rm asymp}$ where Eq. (32) would predict a smaller value. We should point out that the choice of $T_{\rm asymp}$ plays only a minor role in our modeling, as most of the important physical effects occur in the optically thick part of the wind.

The temperature at the sonic point is

$$T_{\rm s}^{4} = \frac{3}{4} \mbox{$T_{\rm eff}$ }^{4}(R_{\rm s}) \Bigl( \tau_{\rm s}^{\prime} + 2 / 3 \Bigr),$$ (34)

with an effective optical depth $\tau_{\rm s}^{\prime}$ at the sonic point

$$\tau_{\rm s}^{\prime} = \int_{R_{\rm s}}^{\infty} \chi \rho \frac {R_{\rm s}^{2}} {r^{2}}{\rm d}r.$$ (35)

Notice that $T_{\rm s}$ depends on the effective optical depth at the sonic point and hence on the run of $\chi \rho$ at all layers above $R_{\rm s}$. Therefore we can only find the solution of the momentum and energy equations at the sonic point if we know the functions $\chi(r)$ and $\rho(r)$ at $r>R_{\rm s}$.

For determining the distribution of the flux-mean opacity in the wind, we use the approach of Lucy & Abbott (1993) who derived $\chi(r)$ from the demand of the momentum conservation in the wind for an adopted velocity law. This requirement follows from Eqs. (8) and (16):

$$\chi = \frac{4 \pi c r^{2}}{L(r)} \left\{ \left( v-\frac{a_{... ...ac{GM}{r^{2}} + a_{1} \frac{{\rm d}T}{{\rm d}r} \right\}\cdot$$ (36)

So for a given v(r), T(r) and L(r) the function $\chi(r)$ is known. Note that $\chi$ and v are related so that $\chi$ increases with the increase of v. The term consisting of ${\rm d}T/{\rm d}r$ is very small and is found by using the diffusion approximation for $r \leq 1.01 R_{\rm s}$ and by differentiating the temperature formula for $r > 1.01 R_{\rm s}$.

The value of L(r) is very close to $L_{\rm core}$ throughout the optically thick and optically thin part of the wind. This is because the other terms in the energy Eq. (19), i.e. $L_{\rm adv}$ and the flow of potential, kinetic energy and enthalpy, are all very small compared to $L_{\rm core}$ for WR-stars. So in principle we could have used $L(r)=L_{\rm core}$ for the calculation of $\chi(r)$ with Eq. (36). In the computations we took the other terms into account properly for the calculation of L(r). For the value of $L_{\rm adv}$ in the calculation of L(r) we used the following scheme: for the optically thick part we calculated $L_{\rm adv}$ from Eq. (18) and for the optically thin part we used $L_{\rm adv}(r)\simeq 2 v(r) L(r)/c$, which is valid if $J \simeq H\simeq K$ (see Eq. (12)). In any case, the exact value of $L_{\rm adv}$ makes a difference of less than 1 percent in the calculation of $\chi(r)$.

   
2.5 The velocity structure

The density in the supersonic part of the wind, which is needed for the calculation of $\chi$ and T(r), can be expressed in terms of the mass-loss rate and the velocity law. The velocity law in the supersonic part of the wind is approximately a $\beta$-law:

$$v(r)~\simeq~ \mbox{$v_\infty$ }~\left( 1- \frac{R_{\rm s}}{r}\right)^{\beta}$$ (37)

with the value of $\mbox{$v_\infty$ }$ taken from the observations and the value of $\beta$ still to be determined. However, near the sonic point the velocity must deviate from this $\beta$-law, because it would imply that $v(R_{\rm s})=0$, which contradicts the definition of the sonic point. Therefore we adopt a velocity law in the region of $r \ge R_{\rm s}$ of the type

 

$$v(r) = \left\{ w_1\left(\frac{R_{\rm s}}{r}\right) + w_2\left(\fr... ...\right)^3 \right\} + v_{\infty} \left( 1 - \frac{R_{\rm s}}{r} \right)^{\beta},$$ (38)

with $w_1+w_2+w_3=v_{\rm s}$. The $\beta$-law is only valid for $r>R_{\rm s}$. This expression for v(r) is practical for our purpose because (i) the first three terms are small and quickly vanish for $r>R_{\rm s}$ and (ii) the critial point conditions of Eqs. (26) and (27) directly describe the values of w₁ and w₂ with w₃ adjustable to give $v(R_{\rm s})=v_{\rm s}$.

We stress that the second part of the expression for v(r), i.e. the $\beta$-law, only enters into our analysis for the calculation of the velocity and density structure in the supersonic part of the wind. For the analysis of the sonic point conditions only the first part is important. This is equivalent to the statement that the $\beta$-law does not start exactly at $R_{\rm s}$ but at some slightly larger distance where the power-law part slowly merges with the $\beta$-law. (In principle we could have chosen a more elaborate expression for the velocity law that ensures a smooth transition from the subsonic to the supersonic part, than the one adopted here, but that would not have improved the analysis of the sonic point conditions.)

   
2.6 The opacity gradient near the sonic point

If the velocity law and the temperature structure are known, as described above, then the product $\chi(r)L(r)$ can be found from Eq. (36). The solution of the critical point equations, which will give the mass-loss rate, requires knowledge of the gradients of the optical depth and of the radiative luminosity near the sonic point. In normal modeling of a stellar wind one would assume that $\chi=f(\rho,T)$ is known and then solve the structure of the wind and find the mass-loss rate. However, as argued in the introduction, the flux-mean opacity in the atmospheres of WR-stars is not known. Therefore, the purpose of our paper is just the opposite: we try to find what the conditions for $\chi$ are that could explain the observed mass-loss rates in terms of optically thick radiation driven winds. A consequence of this approach is that we do not know $\chi(r)$, nor its gradient near the sonic point. Instead, we derive it from the requirement that the model must produce the observed mass-loss rate.

If the gradient of L(r) is known we can derive the gradient of $\chi(r)$ from the derivative of Eq. (36) and vice versa. Therefore we adopt two methods to estimate these gradients: variants A and B. In the first one we make an assumption about L(r) and find the gradient of $\chi(r)$. In variant B we make an assumption about the variation of $\chi(r)$ and find the gradient of L(r). In our modeling we will use both variants and show that the results are quite similar. We describe these two variants in more detail in the next section.