A&A 377, 522-537 (2001)
DOI: 10.1051/0004-6361:20011115

Self-consistent coupling of radiative transfer and dynamics in dust-driven winds

S. Liberatore¹ - J.-P. J. Lafon¹ - N. Berruyer²

1 - Observatoire de Paris-Meudon, DASGAL/UMR 8633, 92195 Meudon, France
2 - Observatoire de la Côte d'Azur, Laboratoire G. D. Cassini/UMR 6529, 06304 Nice, France

Received 12 May 2000 / Accepted 3 August 2001

Abstract
Models of stationary dust-driven winds of late-type stars are investigated. The flow is described successively using three models which couple, in spherical symmetry, the grain-gas dynamics in a self-consistent way with radiative transfer. Complete radiative transfer including multiple scattering, absorption and thermal emission is taken into account to determine the temperature of dust grains which in turn governs their thermal emission. The medium is not necessarily optically thin. The first model is used to check one classical hypothesis, that where the gas and the grains expand at the same velocity, the flow is described by a one-fluid model. This model is then improved in a second model, to include complete momentum coupling between gas and grains by friction. Finally, a third model includes grains and gas coupled by friction as well as the effects of inertial force on grains. By means of a numerical iteration, dynamics and radiative transfer are coupled in order to achieve a self-consistent solution in all cases. Even for fairly low non-zero optical depths, coupling of radiation with dynamics is found to be important for wind models which are all highly sensitive to input data. In conclusion, approximations (position and momentum coupling) for the dynamics should be relaxed.

Key words: stars: mass loss - circumstellar matter - hydrodynamics - radiative transfer

1 Introduction

The existence of circumstellar dust grains and mass loss in cool stars is by now well established (Cassinelli 1979; Zuckerman 1980). A wealth of papers have been devoted to the investigation of the hydrodynamical and thermodynamical role of dust in the circumstellar environment (see for instance Winters et al. 2000 and the reviews of Lafon & Berruyer 1991 and Sedlmayr & Dominik 1995, and references therein).

Whatever the mechanisms causing the wind, the mass loss is controlled by the radiative pressure on dust grains. Photons emitted by the star are absorbed by grains and partially scattered in a more or less isotropic way. A part of the momentum carried by the absorbed photons is transferred to the grains. Friction between grains and atoms or molecules then transmits part of the grain's momentum to the gas.

Much effort was devoted to the study of grain formation and growth (Salpeter 1974a; Draine 1979; Deguchi 1980; Lafon 1991; Patzer et al. 1998), and the construction of mass outflow models (Salpeter 1974b; Kwok 1975; Lucy 1976; Menietti & Fix 1978; Phillips 1979; Deguchi 1980; Winters et al. 2000).

Many models assume (following Gilman 1972) that the whole momentum transferred to the dust grains from the stellar radiation is transferred to the gas by collisions. This "momentum coupling assumption'' provides a highly simplified method of calculating grain and gas velocities, but its influence on the overall structure of the wind models has never really been quantitatively checked.

Krüger et al. (1994) and Krüger & Sedlmayr (1997) built a two-fluid model for stationary dust-driven winds where the dust and the gas components are coupled through the condensation of dust from the gas phase and grain-gas collisions. The radiative transfer on dust is derived from a semi-analytical method developed by Lucy (1971). Moreover, the medium used in the model is optically thin, and the radiative transfer is not self-consistent with the dynamics.

Another work that is close to our purpose has been carried out by Simis et al. (2001). They have combined time-dependent hydrodynamics with a two-fluid model for dust-driven AGB winds. However, they estimate the gas-grain momentum transfer from each time step to the next before which the terminal velocity is assumed to be reached (as a function of some mean free path typical for this process). They also assumed a grey dust opacity.

Tielens (1983) or Netzer & Elitzur (1993) do not explicitly solve the equation of motion for the dust component. The dust is assumed to move with its equilibrium drift velocity relative to the gas.

In order to remove highly inaccurate approximations of dynamics and radiative transfer, we choose hereafter to limit the number of mechanisms to the minimum necessary. Radiative transfer is treated for all frequencies, self-consistently. Simplifying assumptions such as grey opacity, isotropic scattering or the Eddington approximation are relaxed. The question of dynamics is treated for three types of models based respectively on the position cloupling assumption, the momentum coupling assumption and finally on the full problem where the inertial force on grains is taken into account.

The equations for the three dynamical models and for radiative transfer are displayed in Sect. 2. Solutions of dynamical equations (in the three cases), radiative transfer and the coupling between the two problems are explained in Sect. 3. In Sect. 4 we discuss and compare self-consistent models of dust driven wind (Position Coupling, Momentum Coupling and Full Problem). The main conclusions are given in Sect. 5. They will be developed more extensively in future papers.

2 Wind models

2.1 Assumptions

The following simplifying assumptions are used: the flow is stationary, spherically symmetric and grains exist above some reference level r₀, which we refer to as "the base of the wind''. The central star is a blackbody at a given temperature $T_{\ast }$. The wind is neutral and isothermal. In fact, the scale height of the gas temperature is very much larger than any other scales in the considered envelope (Berruyer 1991). The balance between radiation absorbed and emitted by the grains determines the equilibrium grain temperature.

There is an energy equation for the dust alone. To determine the gas temperature profile, it would be necessary to solve the radiative transfer in the gas, which is in fact never optically thick enough to require this. The gas and dust temperatures are consistent in all processes. Variations in gas temperature, however, can induce variations of the nucleation and growth rates, and subsequently their optical properties. Finally, this creates feedback effects on radiative transfer and consistent dynamics. The variations in gas temperature in the considered range of radial distances are very small and hereafter all the effects will be ignored.

Berruyer (1991) studied, for a set of assumed kinetic temperature profiles, the stability of a static model of envelope in the presence of dust grains. Conditions under which the presence of dust grains may generate an instability in a static envelope and drive the motion of the gas component were determined: convectively unstable configurations exist, in spite of the mechanical equilibrium. Furthermore, Berruyer (1991) showed that for envelopes where temperature decreases less steeply than 1/r, the hydrostatic equilibrium is impossible in the outer layers, and the dust-driving is an extra mechanism which essentially modifies the velocity profile and the mass loss rate. Nevertheless, stricto sensu, the dust grains do not "create'' the wind.

Both time-dependent effects and thermal coupling between gas and grains will be neglected. The dust grains are spherical in shape and uniform in density. They scatter almost isotropically, i.e., only first order effects in anisotropic scattering are taken into account.

We can neglect electrostatic drag, which is about a factor 10⁴ smaller than the collisional drag for a cool star (Tabak et al. 1975). We have neglected the fragmentation and coalescence effects (Scalo 1977; Simpson et al. 1980; Biermann & Harwitt 1980) and have ignored grain pressure. Indeed, if the grains were treated as a perfect gas, with a kinetic temperature similar to the gas kinetic temperature, grain pressure for any reasonable wind model would be smaller than gas pressure by about a factor of 10^-14 and would therefore have no great effect.

In this paper, three models corresponding to three descriptions of the dynamics problem are presented (PC, MC and FP), and in each case dynamics are fully coupled with radiative transfer (described bellow). The results are compared in Sect. 2.3.

2.2 Dynamical equations

2.2.1 Full problem

The transfer of momentum gained by the grains through radiative force to the gas is reduced at least in the inner layers, since a part of the radiative force is used to overcome inertia and gravitational forces. the flow is described by a set of two differential equations, one for the grains and one for the gas, coupled by collisional drag.
The motion of the gas is described by the equation of mass conservation (stationary):

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

and momentum conservation:

$$\rho v \frac{{\rm d}v}{{\rm d}r} = -\frac{{\rm d}p}{{\rm d}r} - \rho \gamma_{\rm {grav}} + n_{\rm {gr}} F_{\rm {frict}}.$$ (2)

The state equation for the gas may be written as:

$$p=\rho C_{\rm s}^2 \qquad \textrm{with} \qquad C_{\rm s}^2=\frac{k T_0}{\mu m_{\rm H}},$$ (3)

where T₀ is the constant gas temperature, $\mu$ the number of atoms per free particle, k the Boltzmann constant, $m_{\rm H}$ the mass of the hydrogen atom and $C_{\rm s}$ the sound velocity.

The motion of the grains is submitted to the equations of mass conservation:

$$4 \pi r^2 \rho_{\rm {gr}} v_{\rm {gr}} = \dot{M}_{\rm {gr}}=\textrm{const.},$$ (4)

and momentum conservation:

$$\rho_{\rm {gr}} v_{\rm {gr}} \frac{{\rm d}v_{\rm {gr}}}{{\r... ...rm {gr}} \gamma_{\rm {grav}} - n_{\rm {gr}} F_{\rm {frict}},$$ (5)

where r is the radial distance, and $\rho$ ( $\rho_{\rm {gr}}$), v ( $v_{\rm {gr}}$), p, $\dot{M}$ ( $\dot{M_{\rm {gr}}}$) are the gas density (grain density), the gas velocity (grain velocity), the pressure and the mass loss rate of the gas (of dust grains) respectively. $n_{\rm {gr}}(r)$ is the number density of dust grains in particle per unit volume.

$F_{\rm {rad}}(r)$ is the radiative force on each grain.

$F_{\rm {frict}}(r, v, v_{\rm {gr}})$ is the frictional force which couples the gas to each grain.

$\gamma_{\rm {grav}} (r)$ is the gravitational acceleration.

The mass density of the grains $\rho_{\rm {gr}}$ is defined by:

$$\rho_{\rm gr}(r) = m_{\rm gr} n_{\rm gr}(r) \qquad \textrm{with} \qquad m_{\rm gr}=\frac{4}{3} \pi \rho_{\rm s} a^3,$$ (6)

where $m_{\rm gr}$ and a are the mass and the radius of a grain, respectively and  $\rho_{\rm s}$ the bulk density of the grain material. This full problem has been studied for $F_{\rm rad}(r) \propto r^{-2}$ by Berruyer & Frisch (1983), where we can find, among others, discussion about expressions of the frictional force $F_{\rm {frict}}$.

The fact that the size of the grains is much smaller than the mean free-path of the gas particles implies that the classical Stoke's law from Batchelor (1967) does not apply. The collisional drag for circumstellar grains has been carefully calculated by Baines et al. (1965). The transfer of momentum assuming specular reflexion of the gas particles on the grains ensures the maximum transfer.

The expression of the frictional force is:

$$F_{\rm frict}(r, v, v_{\rm {gr}}) = 2 \sqrt{\pi} C_{\rm s}^2 \rho a^2 F(V_{\rm d}),$$ (7)

with:

 

$${F(V_{\rm d}) = {\rm sgn} (V_{\rm d}) \left( \mid V_{\rm d}\mid +... ...V_{\rm d}^2 + 1 - \frac{1}{4 V_{\rm d}^2}\right) {\rm erf}(\mid V_{\rm d}\mid),$$ (8)

and:

$$V_{\rm d} = \frac{1}{\sqrt{2}}\bigg( \frac{v_{\rm gr} - v}{C_{\rm s}}\bigg).$$ (9)

$V_{\rm d}$ is essentially the drift velocity of a grain relative to the gas, in units of the sound velocity on the gas. Equations (7), (8) and (9) correspond to the formula of Schaaf (1963) for the particular case of elastic gas-grain collisions and large Knudsen numbers. Expression (8) is not handy for numerical calculations and Berruyer & Frisch (1983) constructed an accurate analytical approximation:

$$F(V_{\rm d}) = {\rm sgn} (V_{\rm d})\sqrt{\pi}\left( 1 + V_{... ... - \frac{1}{(aV_{\rm d}^2 + bV_{\rm d} + 1)^{\alpha}}\right).$$ (10)

Equation (10) enables us to retrieve the usual linear law for the subsonic region of the outflow, and the quadratic law for the supersonic part of the wind.

For our forthcoming models, it is convenient to write the frictional force solely in terms of velocities in order to limit the number of independent variables in the momentum conservation equations. With Eq. (1) we have thus:

$$F_{\rm frict}(r, v, v_{\rm {gr})} = \frac{{\cal A} F(V_{\rm d})}{r^2 v},$$ (11)

where:

$${\cal A} = \frac{2 C_{\rm s}^2 a^2 \dot{M}}{4 \sqrt{\pi}}\cdot$$ (12)

Two boundary conditions are required. Any reasonable stellar wind model should be in pressure balance with the interstellar medium at infinity. Since the latter is small compared to pressure anywhere inside the envelope, we shall assume the gas pressure p to be zero at infinity. It may be noted that at large radii, the wind will interact with the interstellar medium, giving rise to a breakdown in the present model assumptions (see for example Weymann 1960). The transition from stellar wind to interstellar medium is not studied here. This transition is usually considered to be a shock wave, for example see Mihalas (1978), p. 533. This shock wave is located at about $1000~{ R_{\ast}}$ and our models are calculated up to a distance about $160~{R_{\ast}}$ within which the gas remains isothermal.

For cool stars, grains are presumably generated inside the flow. Our model starts at the moment the grains formed. It therefore appears reasonable to assume that at the base of the flow the gas and the grains have the same velocity. This will be our second boundary condition.

A discussion about the behaviour of the flow when grains appear is given by Berruyer & Frisch (1983). Two different dynamical regimes exist, a boundary layer regime where the grains are submitted to a very strong acceleration (the gas keeping an almost constant velocity), and an outer regime where the two fluids are strongly coupled. Our numerical results scrutinize the velocity of grains just outside the boundary layer because this boundary layer is very thin and cannot be computed.

2.2.2 Position coupling

This approximation corresponds to a strongly local coupling where there are enough collisions between gas and grains for the two phases to move together. The grains are thus "frozen'' in the gas; there is only one fluid velocity:

$$v(r) = v_{\rm gr}(r),$$ (13)

and thus just one equation for momentum.

Addition of Eqs. (5) and (2) of momentum conservation, using Eq. (13) leads to:

$$(\rho + \rho_{\rm gr}) v \frac{{\rm d}v}{{\rm d}r} = -\fr... ...rho_{\rm gr}) \gamma_{\rm {grav}} + n_{\rm gr} F_{\rm {rad}}.$$ (14)

This is a one fluid model equation since $\dot{M}= 4 \pi r^2 \rho v$, $\dot{M}_{\rm gr}= 4 \pi r^2 \rho_{\rm gr} v$. However, to get solutions more comparable to the other cases (MCH and FPB) we prefer to consider the (realistic) case where:

$$\rho_{\rm gr} \ll \rho.$$ (15)

Then, Eq. (14) reduce to another one-fluid equation:

$$\rho v \frac{{\rm d}v}{{\rm d}r} = -\frac{{\rm d}p}{{\rm d}r} - \rho \gamma_{\rm {grav}} + n_{\rm gr} F_{\rm {rad}}.$$ (16)

Boundary conditions are the same as those for the gas in the full problem.

2.2.3 Momentum coupling

Gas and grains are coupled by friction. Under this approximation, the whole momentum given to grains by the radiative force is transferred to the gas by grain-gas collisions. In other words, the radiative force is set equal to the frictional force (the grain mass is taken to be equal to zero):

$$F_{\rm {rad}}(r) = F_{\rm {frict}}(r, v, v_{\rm {gr}}).$$ (17)

This approximation is also called the complete momentum coupling.

One can also use another approximation in which the gravity of grains is taken into account. we shall call it the "Partial Momentum Coupling'' because only part of the momentum given to the grains by radiative pressure is released to the gas by collisions. The other part counter-balance the gravitationnal force and:

$$F_{\rm rad}(r)= F_{\rm {frict}}(r, v, v_{\rm {gr}}) + m_{\rm gr} \gamma_{\rm grav}(r).$$ (18)

Comparing the results obtained using these approximations allows us to study the effect of gravity on the grains.

Boundary conditions are the same as in full problem except for those concerning initial velocities. In all cases except in the full problem, initial grain and gas velocities are not equal because the boundary layer is removed by degeneracy of the differential equation for the grains.

   
2.3 Formulation of the radiative transfer problem

The present work is based on the "Quasi-Diffusion Method'' developped by Leung (1975) for radiation transport in dense interstellar dust clouds. The boundary conditions have been modified accordingly to include a star with finite extent at the center of an extended dust shell. Hereafter we summerize the results of the conversion.

2.3.1 Equations of the model

The time-independent equation of radiative transfer for anisotropic scattering and isotropic source can be written as:

 

$${\mu \frac{\partial I_{\nu}(r, \mu)}{\partial r} + \frac{1-\mu^2}... ...frac{1}{2}\int_{-1}^{+1} I_{\nu}(r, \mu') p_{\nu}(\mu, \mu'){\rm d}\mu'\right],$$ (19)

where $I_{\nu}(r, \mu)$ is the specific intensity of radiation with frequency $\nu$ at some radial distance r and in the direction $\mu$, $B_{\nu}(r)$ is the Planck function, $\sigma_{\nu}^a(r)$ is the absorption coefficient per unit volume, $\sigma_{\nu}^s(r)$ is the scattering coefficient per unit volume, and $p_{\nu}(\mu, \mu')$ is the phase function which determines the the probability of scattering from direction $\mu'$ into the direction $\mu$, for photons with frequency $\nu$. $\mu$ is the cosine of the angle between the outwards directed radius and the direction of photon propagation at point r. We now define the quantities:

1.

$I^+_{\nu}$ for $0 \leq \mu \leq 1$ and $I^-_{\nu}$ for $-1 \leq \mu \leq 0$. Physically $I^+_{\nu}$ and $I^-_{\nu}$ refer to radiation propagating along the positive and negative directions of the line of sight.

2.

A quantity $u_{\nu}$ which has a mean-intensity-like character, and $v_{\nu}$ which has a flux-like character, where:

  

$$u_{\nu} \equiv \frac{1}{2}(I^+_{\nu} + I^-_{\nu})$$ (20)

$$v_{\nu} \equiv \frac{1}{2}(I^+_{\nu} - I^-_{\nu}).$$ (21)

3.

The moment integrals defined in terms of $u_{\nu}$ and $v_{\nu}$:

$$J_{\nu} \int_0^1 u_{\nu}(r, \mu){\rm d}\mu;$$ = (22)

$$H_{\nu} \int_0^1 v_{\nu}(r, \mu) \mu {\rm d}\mu;$$ = (23)

$$K_{\nu} \int_0^1 u_{\nu}(r, \mu) \mu^2 {\rm d}\mu.$$ = (24)

The thermal radiation of the dust grains has been assumed to be isotropic; but not the scattering, which occurs mainly in a forward direction. The net effect of the anisotropic scattering is a reduction in the degree of multiple scattering of a given photon through the shell. Indeed, a photon can be absorbed after it has been scattered.

The discussion on anisotropic scattering can be found in Leung's (1975) paper. According to Leung (1975), and using the first order asymmetry parameter $g_{\nu}$ (computed using the Mie scattering theory) which is equal to zero for isotropic scattering, to 1 for forward scattering and to -1 for backward scattering, the zero and first order moments of Eq. (19) with respect to $\mu$ are:

$$\frac{1}{r^2}\frac{\partial}{\partial r}\left[r^2 H_{\nu}(r) \right] = - \sigma_{\nu}^a(r)\left[J_{\nu}(r)-B_{\nu}(r)\right]$$ (25)

 

$${\frac{\partial K_{\nu}(r)}{\partial r} + \frac{3K_{\nu}(r)-J_{\n... ...igma_{\nu}^a(r) + \left[1- g_{\nu} \right]\sigma_{\nu}^s(r) \right\}H_{\nu}(r).$$ (26)

Using the variable Eddington factor:

$$f_{\nu}(r) \equiv \frac{K_{\nu}(r)}{J_{\nu}(r)}$$ (27)

introduced by Auer & Mihalas (1970) and the sphericity function:

$$\zeta_{\nu}(r) \equiv \textrm{exp} \left[ \int_{r_0}^r \left( 3 - \frac{1}{f_{\nu}(x)} \right)\frac{{\rm d}x}{x}\right],$$ (28)

suggested by Auer (1971), we combined Eqs. (25) and (26) to yield the combined moment equation:

1.

Substituting $f_{\nu}$ and $\zeta_{\nu}$ into Eq. (26), we get:

$$\frac{\partial}{\partial r}\left(f_{\nu}\zeta_{\nu}J_{\nu}\right) = - \kappa_{\nu} \zeta_{\nu} H_{\nu},$$ (29)

where $\kappa_{\nu}(r) = \sigma_{\nu}^a(r) + \left[1-g_{\nu} \right]\sigma_{\nu}^s(r)$.

2.

Eliminating $H_{\nu}$ from Eq. (25), we obtain the combined moment equation:

$$\frac{\partial}{\partial r}\left[ - \frac{r^2}{\kappa_{\nu}\... ...ft(\kappa^{\ast}_{\nu}J_{\nu} - \epsilon^{\ast}_{\nu} \right),$$ (30)

where $\kappa^{\ast}_{\nu} \equiv \sigma^a_{\nu}$ and $\epsilon^{\ast}_{\nu} \equiv \sigma^a_{\nu}B_{\nu}$.

If $Q_{\rm abs}$, $Q_{\rm sca}$ are the absorption and scattering efficiency factors in the Mie theory and $n_{\rm gr}$ denotes the number density of dust grains, then we have:

$$\kappa_{\nu}(r) n_{\rm gr}(r)\left[Q_{\rm abs}(\nu) \pi a^2+ (1-g_{\nu})Q_{\rm sca}(\nu)\pi a^2\right]$$ = (31)

$$\epsilon^{\ast}_{\nu}(r) n_{\rm gr}(r) Q_{\rm abs}(\nu) \pi a^2 B_{\nu}(r)$$ = (32)

$$\kappa^{\ast}_{\nu}(r) n_{\rm gr}(r) Q_{\rm abs}(\nu) \pi a^2.$$ = (33)

The variable Eddington factor and sphericity function are evaluated from an angle-dependent formal solution of the transfer equation, obtained by, integrating it along parallel directions tangential to spherical shells as described in Hummer & Rybicki (1971). Following the technique introduced by Feautrier (1964), they used a system of cylindrical coordinates $(\xi,z)$ where:

$$r\mu$$ z = (34)

$$\xi$$ = (r² - z²)^1/2. (35)

According to Leung (1975), with these coordinates we obtain two coupled first-order differential equations:

 

$$\frac{{\rm d} u_{\nu}}{{\rm d} z} - \left( \sigma^a_{\nu} + \sigma^d_{\nu}\right) v_{\nu};$$ = (36)

$$\frac{ {\rm d}v_{\nu}}{{\rm d} z} - \left( \sigma^a_{\nu} + \sigma^d_{\nu}\right)\left( u_{\nu} - S_{\nu}\right);$$ = (37)

from which we eliminate $v_{\nu}$ to obtain a simple second-order differential equation:

$$\frac{1}{\sigma^a_{\nu} + \sigma^d_{\nu}} \frac{{\rm d}}{{\r... ...frac{{\rm d} u_{\nu}}{{\rm d} z} \right) = u_{\nu} - S_{\nu},$$ (38)

where:

$$S_{\nu}(\mu) = \frac{1}{(\sigma^a_{\nu} + \sigma^d_{\nu})} \left[... ...1} I_{\nu}(\mu')p_{\nu}(\mu, \mu'){\rm d}\mu' + \sigma^a_{\nu} B_{\nu} \right].$$ (39)

With the above mentioned anisotropic scattering approximation using the $g_{\nu}$ parameter, according to Leung (1976), we can derive:

$$S_{\nu}(r, \mu) = \frac{\eta_{\nu}(r,\mu)}{\chi_{\nu}(r)},$$ (40)

with:

$$\chi_{\nu}(r) n_{\rm gr}(r)\left[Q_{\rm abs}(\nu) \pi a^2+ Q_{\rm sca}(\nu)\pi a^2\right],$$ (41)

$$\eta_{\nu}(r,\mu) = n_{\rm gr}(r) Q_{\rm abs}(\nu) \pi a^2B_{\nu}... ...(r) Q_{\rm sca}(\nu)\pi a^2 \left[ J_{\nu}(r) +3 g_{\nu} \mu H_{\nu}(r)\right].$$ (42)

The solution of Eq. (38) with boundary conditions at the inner and outer radius of the dust shell yields the angular dependence of $u_{\nu}$.

The evaluation of $f_{\nu}$ and $\zeta_{\nu}$ follows from Eqs. (27) and (28). The equation of radiative equilibrium which specifies the temperature $T_{\rm gr}(r)$, is given by (radiative equilibrium):

$$\int_0^{\infty}\left[ \epsilon_{\nu}^{\ast}(r)-\kappa_{\nu}^{\ast}(r)J_{\nu}(r)\right]{\rm d}\nu = 0.$$ (43)

2.3.2 Boundary conditions

The radiation coming from the interstellar medium was taken zero:

$$I^-_{\nu}(r_{\rm f}) = 0.$$ (44)

$\bullet$ For moment equations.

The following quantities are defined at each boundary $r=r_{\rm B}$:
Downward flux:

$$H^-_{\nu}(r_{\rm B}) \equiv \frac{1}{2}\int_0^1 I_{\nu}(r_{\rm B}, -\mu)\mu{\rm d}\mu;$$ (45)

Upward flux:

$$H^+_{\nu}(r_{\rm B}) \equiv \frac{1}{2}\int_0^1 I_{\nu}(r_{\rm B}, +\mu)\mu{\rm d}\mu;$$ (46)

Boundary factor:

$$f^{\rm B}_{\nu} \equiv \frac{1}{2} \frac{\displaystyle \int_0^1 \left[ I_{\nu}(r_{\rm B}... ..._{\nu}(r_{\rm B}, -\mu\right] \mu{\rm d}\mu }{\displaystyle J_{\nu}(r_{\rm B})}$$ (47)

$$\frac{\displaystyle \int_0^1u_{\nu}(r_{\rm B}, \mu)\mu{\rm d}\mu } {\displaystyle \int_0^1u_{\nu}(r_{\rm B}, \mu){\rm d}\mu}\cdot$$ = (48)

Usually $H_{\nu}^-$ is known at the outer boundary $r_{\rm f}$ while $H_{\nu}^+$ is partialy known at the inner boundary r₀. The net flux at each boundary is noted  $H_{\nu}^{\ast}$. It can be expressed in terms of the above quantities:

$$H_{\nu}^{\ast}(r_{\rm B}) = H_{\nu}^+(r_{\rm B}) - H_{\nu}^-(r_{\rm B}),$$ (49)

where $H_{\nu}^+(r_{\rm B})$ is calculated with Eq. (46) and  $H_{\nu}^-(r_{\rm B})$ with Eq. (45).

At the outer boundary, $r=r_{\rm f}$, Eq. (26) together with Eqs. (27) and (28) yields:

$$\frac{\partial}{\partial r}\left(f_{\nu} \zeta_{\nu} J_{\nu} ... ...^{\rm B}_{\nu} J_{\nu} + 2 \kappa_{\nu} \zeta_{\nu} H_{\nu}^-.$$ (50)

In our model no radiation comes from the outside. Thus $H_{\nu}^- = 0$ and the following equation can be taken as a boundary condition:

$$\frac{\partial}{\partial r}\left(f_{\nu} \zeta_{\nu} J_{\nu} \right) = - \kappa_{\nu} \zeta_{\nu} f^{\rm B}_{\nu} J_{\nu}.$$ (51)

The boundary factor $f^B_{\nu}$ can be evaluated from the solution of ray equations (if the radiation field is isotropic at the boundary, $f^B_{\nu} = \frac{1}{2}$). Note that to get $f_{\nu}$ and $f^B_{\nu}$ we need to evaluate only $u_{\nu}$, and not both $u_{\nu}$ and $v_{\nu}$.

At the inner boundary r=r₀, in the same way the boundary condition can be expressed as:

$$\frac{\partial}{\partial r}\left(f_{\nu} \zeta_{\nu} J_{\nu} \right) = - \kappa_{\nu} \zeta_{\nu} H_{\nu}^{\ast},$$ (52)

where $H_{\nu}^{\ast}$ is the net flux, at frequency $\nu$, at r=r₀. The radiative flux $H_{\nu}(r_0)$ at the "base of the wind'' r₀ is, by continuity, equal to $H_{\nu}^{\ast}(r_0)$ derived from Eq. (49).

There are various contributions to $H_{\nu}^{\ast}(r_0)$:
- The flux $H_{\nu}^+(r_0)=H_{\nu \ast}^+(r_0)+H_{\nu \rm shell}^+(r_0)$, where $H_{\nu \ast}^+(r_0)=(r_{\ast}^2/4 r_0^2) B_{\nu}(T_{\ast})$ is the flux of photons from the star, assumed to be that of a Planck function $B_{\nu}(T_{\ast})$ (blackbody of temperature $T_{\ast }$) with a dilution factor of $\frac{1}{4}(r_{\ast}/r)^2$, where $r_{\ast}$ is the stellar radius. $H_{\nu \rm shell}^+(r_0)$ corresponds to the re-emitted and scattered photons from the dust shell which cross the inner cavity ( $r_{\ast} < r < r_0$) and re-enter the dust shell. The contribution of $H_{\nu \rm shell}^+(r_0)$ is taken into account by Leung (1976) in the radiative transfer calculation.
- The flux $H_{\nu}^-(r_0)$ correspond to the re-emitted and scattered photons from the dust shell which cross the inner boundary in the reverse direction.

$\bullet$ For ray equations.

The boundary conditions associated with Eq. (38) are:
At the surface $r=r_{\rm f}$,

$$\frac{1}{\sigma^a_{\nu} + \sigma^d_{\nu}}\frac{{\rm d} u_{\nu}}{{\rm d} z} = - u_{\nu},$$ (53)

and at $r=\xi$, if the ray does not cross the core (i.e. lie in the circumstellar envelope) and because of the symmetry, implying $I^+_{\nu}(r, \xi) = I^-_{\nu}(r, \xi)$, we obtain from Eq. (36):

$$\frac{1}{\sigma^a_{\nu} + \sigma^d_{\nu}}\frac{{\rm d} u_{\nu}}{{\rm d} z} = 0.$$ (54)

However if for a given impact parameter a ray crosses the core defined by the inner radius of the dust shell, we then use from Eq. (36):

$$\frac{1}{\sigma^a_{\nu} + \sigma^d_{\nu}}\frac{{\rm d} u_{\nu}}{{\rm d} z} = - v_{\nu}^{\ast},$$ (55)

where $v_{\nu}$ is denoted by $v_{\nu}^{\ast}$ at each boundary.

At each step of the iteration scheme concerning the radiative transfer, $H_{\nu}^{\ast}$ is determined as:

$$H_{\nu}^{\ast}=\int_0^1 v_{\nu}^{\ast} \mu {\rm d}\mu.$$ (56)

Note that it follows directly from the definitions that the emergent flux and intensity are simply given by: $H^+(r_{\rm f})= f^B_{\nu} J_{\nu}(r_{\rm f})$ and $I^+_{\nu}(r_{\rm f}, \mu) = 2 u_{\nu}(r_{\rm f}, \mu)$.

We solve the radiative transfer problem as a two-point boundary value problem. Using two auxiliary functions - the Eddington factor $f_{\nu}$ and the sphericity function $\zeta_{\nu}$, which respectively characterize the anisotropy of the radiation field, the geometry of the problem (and depend on $f_{\nu}$) - we transform the equation of radiation transport into a quasi-diffusion equation. The "Quasi-Diffusion Method'' essentially consists in two parts: assuming the radiation field isotropic, i.e., setting $f_{\nu}$, $f^B_{\nu}$ and $\zeta_{\nu}$, we solve a set of moment and energy-balance equations to determine $J_{\nu}$, $K_{\nu}$ and then the source function. We then find the angular distribution of the radiation field using a ray-tracing method and hence update $f_{\nu}$, $f^B_{\nu}$ and $\zeta_{\nu}$. Even anisotropic scattering can easily be taken into account.

2.3.3 The radiative force

The dimensionless cross section $Q_{\rm pr}(\nu)$ of the radiative force can be written as:

$$Q_{\rm pr}(\nu)=Q_{\rm abs}(\nu) + (1-g_{\nu})Q_{\rm sca}(\nu)$$ (57)

and the radiative force $F_{\rm rad}$ then reads straight forwardly:

$$F_{\rm rad}(r)=\frac{4 \pi}{c}\int_0^{+\infty} \pi a^2 Q_{\rm pr}(\nu)H_{\nu}(r){\rm d}\nu.$$ (58)

Since a and $Q_{\rm pr}(\nu)$ are inputs of the model, $F_{\rm rad}$ is deduced from the solution of the radiative transfer.

In the case of an optically thin medium, the Eddington flux H is that of a blackbody:

$$H_{\nu}(r) =\frac{1}{4}(\frac{r_{\ast}}{r})^2 B_{\nu}(T_{\ast}),$$ (59)

where $r_{\ast}$ is the stellar radius. Equation (58) can be written as:

$$F_{\rm rad}(r)=(\frac{4 \pi}{c})(\pi a^2 {\overline Q_{\rm pr}})\int_0^{+\infty} H_{\nu}(r){\rm d}\nu,$$ (60)

where the averaged $Q_{\rm pr}$ is:

$${\overline Q_{\rm pr}} = \frac{\displaystyle \pi r^2_{\ast}}... ...} \int_0^{+\infty} Q_{\rm pr}(\nu)B_{\nu}(T_{\ast}){\rm d}\nu.$$ (61)

If we use the general expression of the Eddington flux H (for optically thin and thick mediums):

$$\int_0^{+\infty}H_{\nu}(r){\rm d}\nu=H(r) =\frac{L_{\ast}}{(4 \pi)^2 r^2},$$ (62)

where $L_{\ast }$ is the stellar luminosity. The radiative force becomes:

$$F_{\rm rad}(r) =\frac{\displaystyle a^2 L_{\ast}{\overline Q_... ...splaystyle 4 c} \frac{\displaystyle 1}{\displaystyle r^2}\cdot$$ (63)

In terms of the dynamical dimensionless radius (R=r/r₀), the radiative force reads:

$${\cal F}_{\rm rad}(R)=\frac{F_{\rm rad_0}}{R^2}$$ (64)

where $F_{\rm rad_0}$ denotes the radiative force at the inner boundary when the medium is optically thin.

Since attention is drawn mainly to dynamics we arbitrarily take $a=0.1~\mu$m which is a likely grain size. The values of $Q_{\rm abs}(\nu)$, $Q_{\rm sca}(\nu)$ and $g_{\nu}$ are those for graphite grains from Draine (1985). With these values at $T_{\ast}=3000\,{\rm K}$, ${\overline Q_{\rm pr}}=1.11$.

For optically thin media, iterations start with this expression of the radiative force. It will also be useful to compare the model based on this force profile to the full solution.

In these cases we compute the optical depth at some frequency of reference  $\nu_{\rm ref}$. This frequency has been choosen in the spectral range where the radiative force is most efficient ( $\nu _{\rm ref}=1.3\times 10^{15}\,{\rm Hz}$), that is to say in a spectral zone where $Q_{\rm pr}(\nu) H_{\nu}(r)$ is maximum. Hence, we will denote the optical depth at the frequency reference by  $\tau _{\nu _{\rm ref}}$.

3 Numerical solutions

Due to the coupling of dynamics with radiative transfer, the numerical solution in any of the three considered dynamical regimes requires iterations. After assuming a radiative pressure profile we first look for the velocity profile, then determine the radiative pressure profile from the radiative transfer in a given velocity/density field. The process is then repeated. At each step of iterations (indexed by some integer n), dynamics is governed by an equation with a critical point at $v=C_{\rm s}$, but at various radial distances  $r_{\rm c}^n$. The solution of the problem will be found if convergence of the iterative scheme is obtained. In this case, $r_{\rm c}^n$ tends to a limit say $r_{\rm c}$ (where the velocity is equal to the sound velocity).

3.1 Dynamics

3.1.1 Solution under the position coupling hypothesis (PCH)

 

Using Eq. (16) with Eqs. (1) and (3) we obtain a wind equation of the form:

$$\frac{1}{v} \frac{{\rm d}v}{{\rm d}r} = \frac{\displaystyle ... ...+ {\cal M} F_{\rm {rad}}(r)}{\displaystyle v^2 - C_{\rm s}^2},$$ (65)

where ${\cal M}=\frac{\displaystyle \dot{M}_{\rm gr}}{\displaystyle \dot{M} m_{\rm gr}} =\frac{\displaystyle n_{\rm gr}}{\displaystyle \rho}$. We introduce the dimensionless space and velocity variables:

$$R=\frac{r}{r_0}; \qquad V=\frac{v}{C_{\rm s}}\cdot$$ (66)

where r₀ denotes the radial distance at the base of the wind. Defining the constants $\beta_1= \frac{\displaystyle {\cal M} r_0}{\displaystyle C^2_{\rm s}}$ and $\beta_2 = \frac{\displaystyle G M_{\ast}}{\displaystyle r_0 C_{\rm s}^2}$, the equation can thus be written:

$$\frac{1}{V} \frac{{\rm d}V}{{\rm d}R} = \frac{1}{R^2} \left... ...1}\right) = \frac{\displaystyle f(R)}{\displaystyle V^2 - 1},$$ (67)

where ${\cal F}_{\rm rad}(R) \equiv F_{\rm rad}(r)$.

This is a wind equation similar to that of Parker (1958,1960a,b) for solar wind with a critical point at V=1, $R_{\rm c}=\beta_2/2$ which is a saddle point. The only solution with pressure reaching zero at infinity is the transsonic critical solution that crosses the sonic point (V=1) with an outwardly increasing velocity. Velocity is subsonic at the base of the wind and the gas pressure falls to zero at infinity (Mihalas 1978).

At each step of iteration the wind has the same topological structure as a coronal wind. The topology of the solutions ensures numerical stability, provided the integration is started slightly from the right and left of the critical point.

The critical (sonic) radius $R_{\rm c}$, where $V(R_{\rm c})=V_{\rm c}=1$, is the solution of the equation:

$$f(R_{\rm c}) = \frac{2}{R_{\rm c}} - \frac{\beta_2}{R^2_{\rm c}} + \beta_1{\cal F}_{\rm {rad}}(R_{\rm c}) = 0.$$ (68)

In comparison with the isothermal model of solar wind studied by Parker (1958,1960a,b), we can write Eq. (68) as:

$$R_{\rm c}= R_{\rm C_{Park}} - \frac{\beta_1}{2}{\cal F}_{\rm {rad}}(R_{\rm c}) R_{\rm c}^2,$$ (69)

where $R_{\rm C_{Park}}=\beta_2/2$ is the undimensioned radial distance of the critical point for Parker's isothermal wind model (with our notations).

It is obvious that under the position coupling approximation, the critical distance $R_{\rm c}$ is smaller than  $R_{\rm C_{Park}}$, which, due to the difference of temperatures (absence of corona), does not actually implies that $R_{\rm c}$ is smaller in the case of an evolved star, than for the solar wind.

If the medium is considered as optically thin, the radiative force varies as R^-2, and the position of the sonic point is: $R_{\rm C_{mince}} = R_{\rm C_{Park}} - \frac{\beta_1}{2} F_{\rm rad0}$. Using a radiative pressure derived from a radiative transfer code, we solve Eq. (69) using a classical numerical method. If Eq. (69) has no solution, it is impossible to find a sonic point in the circumstellar envelope. If there is more than one solution, there is a sonic point which exists as a critical point and also a point for which the gradient of the velocity is zero and beyond which the velocity begins to decrease before once again increasing. This case will be investigated in more detail in a forthcoming paper.

Now, we want to study only models with one solution for Eq. (69). The velocity gradient at the sonic point is given by the l'Hospital's rule:

$$\left(\frac{1}{V} \frac{{\rm d}V}{{\rm d}R}\right)_{R_{\rm c}... ... \lim_{R \to R_{\rm c}} \frac{{\rm d}}{{\rm d} R} (V^2 - 1)},$$ (70)

in other words:

$$\left( \frac{{\rm d}V}{{\rm d}R}\right)_{R_{\rm c}}^2 \frac{1}{2} \left(\frac{{\rm d}f}{{\rm d}R}\right)_{R_{\rm c}}$$ =

$$\frac{1}{R_{\rm c}^2} \left(\frac{\beta_2}{R_{\rm c}} - 1\right) ... ..._1}{2} \left(\frac{{\rm d}{\cal F}_{\rm rad}}{{\rm d}R}\right)_{R_{\rm c}}\cdot$$ = (71)

If, $\frac{\displaystyle 1}{\displaystyle R_{\rm c}^2} \left(\frac{\displaystyle \b... ... \frac{\displaystyle {\rm d}{\cal F}_{\rm rad}} {\displaystyle {\rm d}R} \geq 0$, we have a critical (saddle) point with the gradient velocity:

$$\left(\frac{{\rm d}V}{{\rm d}R}\right)_{R_{\rm c}} \sqrt{\frac{1}{2}\left(\frac{{\rm d}f}{{\rm d}R}\right)_{R_{\rm c}}}$$ = (72)

$$\sqrt{\frac{1}{R_{\rm c}^2} \left(\frac{\beta_2}{R_{\rm c}}- 1\right) + \frac{\beta_1}{2} \frac{{\rm d}{\cal F}_{\rm rad}}{{\rm d}R}}\cdot$$ =

With $R_{\rm c},V_{\rm c}$ and $\left( \frac{{\rm d}V}{{\rm d}R}\right)_{R_{\rm c}}$ we can find numerically the velocity profile V(R) using an extrapolation method developed by Bulirsch & Stoer (1966) with variable step sizes. This method gives more accurate results than the classical Runge-Kutta method.

3.1.2 Solution of models including drift velocity

Now we can write Eq. (2) as:

$$\frac{1}{v} \frac{{\rm d}v}{{\rm d}r} = \frac{\displaystyle ... ...ict}}(r, v, v_{\rm gr})}{\displaystyle v^2 - C_{\rm s}^2}\cdot$$ (73)

With the dimensionless variables:

$$R=\frac{r}{r_0}; \quad V=\frac{v}{C_{\rm s}}; \quad V_{\rm gr}=\frac{v}{C_{\rm s}},$$ (74)

expression (11) of the frictional force becomes:

$${\cal F}_{\rm frict}(R, V, V_{\rm gr}) = \frac{\displaystyl... ...rm s}} \frac{\displaystyle F(V_{\rm d})}{\displaystyle R^2 V},$$ (75)

with:

$${\cal F}_{\rm frict}(R, V, V_{\rm gr}) \equiv F_{\rm frict}(r, v, v_{\rm gr})$$ (76)

and:

$$V_{\rm d} = \frac{1}{\sqrt{2}}(V_{\rm gr} - V).$$ (77)

Equation (73) can be rewritten as:

$$\frac{1}{\displaystyle V}\frac{\displaystyle {\rm d}V }{\dis... ...aystyle V^2 - 1}\right) = \frac{f(R, V, V_{\rm gr})}{V^2 - 1}$$ (78)

where $\beta_3 = \frac{\displaystyle {\cal M A}}{\displaystyle r_0 C_{\rm s}^3}\cdot$

To find $R_{\rm c}$, where $V(R_{\rm c})=V_{\rm c}=1$, we need to solve the non-linear equation:

$$\left(f(R_{\rm c}, 1, V_{{\rm gr}_{\rm c}})\right) = \frac{... ...gr}_{\rm c}})}{\displaystyle V_{{\rm gr}_{\rm c}}}\right) = 0,$$ (79)

where $V_{{\rm gr}_{\rm c}} \equiv V_{\rm gr}(R_{\rm c})$ et $F_{\rm c}(V_{{\rm gr}_{\rm c}}) \equiv \left[F(V_{\rm d})\right]_{R_{\rm c}}$. If we define the function:

$$P_{\rm crit}(V, V_{\rm gr}) = R_{{\rm c}_{\rm Park}} - \frac{\beta_3}{2} \frac{F(V_{\rm d})}{V_{\rm gr}},$$ (80)

Eq. (79) can be written as:

$$R_{\rm c} = P_{\rm crit}(1, V_{{\rm gr}_{\rm c}}).$$ (81)

It is necessary to know the velocity $V_{{\rm gr}_{\rm c}}$ of grains at the critical point to have the position of the critical point.

Then we can distinguish two cases.
$\bullet$ Momentum coupling.

Here we analyse the case of momentum coupling using Eq. (17). A priori it would be more fitting to use Eq. (18) which includes gravity. In fact all calculations performed in this case reveal that the difference between the relevant results is very small whereas the equations and computation are much less lengthy. For this reason, we analyse only the case where grain mass is neglected.

Using Eq. (11), with dimensionless variables, the degenerated form of Eq. (17) for grains becomes:

$$\frac{F(V_{\rm d})}{V} = \Phi_1 R^2 {\cal F}_{\rm rad}(R) = g(R)$$ (82)

where:

$$\Phi_1 = \frac{r_0^2 C_{\rm s}}{\cal A} = \frac{2 \sqrt{\pi} r_0^2}{C_s a^2 \dot{M}} = \frac{\beta_1}{\beta_3}\cdot$$ (83)

Equation (81) becomes:

$$P_{\rm crit}(1, V_{{\rm gr}_{\rm c}}) = R_{{\rm c}_{\rm Par... ...gr}_{\rm c}}}\right) {\cal F}_{\rm rad}(R_{\rm c})R^2_{\rm c}.$$ (84)

To find both $R_{\rm c}$ and the grain velocity  $V_{\rm gr_c}$ we must solve the system of two non-linear equations:

$$F_{\rm c}(V_{{\rm gr}_{\rm c}}) = g(R_{\rm c}),$$ (85)

$$R_{\rm c} = P_{\rm crit}(V=1, V_{{\rm gr}_{\rm c}}).$$ (86)

The l'Hospital rule leads to a second degree equation for the gradient velocity of the gas:

$$\left(\frac{{\rm d}V}{{\rm d}R}\right)_{R_{\rm c}}^2 + B(R_{... ...\rm d}R}\right)_{R_{\rm c}} + D(R_{\rm c}, V_{\rm gr_c}) = 0,$$ (87)

with:

$$B(R_{\rm c}, V_{\rm gr_c}) = - \frac{\beta_3}{2}\frac{1}{R_{\rm c... ...r_c}}\right)\left( 1 + \frac{F_{\rm c}}{H_{\rm c}}\right) - H_{\rm c} \right] ,$$ (88)

and:

$$D(R_{\rm c},V_{\rm gr_c}) = -\frac{1}{R_{\rm c}^2}\left( \frac{\b... ...\rm d}g(R)}{{\rm d}R}\right)_{R_{\rm c}} -\frac{2F_{\rm c}}{R_{\rm c}}\right] ,$$ (89)

where:

$$F_{\rm c} \equiv \left(F(V_{\rm d})\right)_{R_{\rm c}},$$ (90)

and:

$$H_{\rm c} \equiv \left(H(V_{\rm d})\right)_{R_{\rm c}}= \fra... ...eft(\frac{{\rm d}F}{{\rm d} V_{\rm d}}\right)_{R_{\rm c}}\cdot$$ (91)

The positive gradient velocity of the gas is thus:

$$\left(\frac{{\rm d}V}{{\rm d}R}\right)_{R_{\rm c}} = \frac{1... ...V_{\rm gr_c}) + \sqrt{\Delta(R_{\rm c}, V_{\rm gr_c})}\right),$$ (92)

with:

$$\Delta(R_{\rm c}, V_{\rm gr_c}) =B^2(R_{\rm c}, V_{\rm gr_c}) - 4 D'R_{\rm c}, V_{\rm gr_c}).$$ (93)

Let us notice that the critical point has the topology of a saddle point only if $\Delta(R_{\rm c}, V_{\rm gr_c}) \ge 0$, and it is solely in this case that it is possible to find a transsonic solution.

Numerical integration can be performed in the case of Position Coupling. A test on the gradient velocity at each point is made to ascertain whether it is positive or negative. In the cases where this quantity is negative, an instability, characterised by an inverse gradient pressure, arises and develops at each iteration step leading to a divergence.

$\bullet$ Full problem.

When the momentum coupling approximation is relaxed, the problem described by Eq. (78) and Eq. (5) reads:

$$\frac{{\rm d}V_{\rm gr}}{{\rm d}R} = \frac{1}{V_{\rm gr}} \l... ... \frac{\Phi_4}{R^2}\left(\frac{F(V_{\rm d})}{V}\right)\right),$$ (94)

where:

$$\Phi_3 = \frac{r_0}{C_{\rm s}^2 m_{\rm gr}} \quad \textrm{and} \quad \Phi_4 = \frac{1}{\Phi_1 m_{\rm gr}}\cdot$$ (95)

The model now becomes much more complicated. We have now two independent variables V and  $V_{\rm gr}$. In the three dimensional space $(R,V,V_{\rm gr})$, there is a line of points solutions of Eq. (81) located in the plane V=1. On part of this line the critical points are saddle points. The solution to our problem will be sought among the critical wind-type solutions which pass through these saddle points. The boundary condition at the base of the wind ( $V=V_{\rm gr}$) will allow us to make a unique determination.

Using Eq. (92), Eq. (81) may be written as:

$$R_{\rm c} = R_{{\rm c}_{\rm Park}} - \frac{\beta_1}{2}\left(\frac... ... \left(\frac{{\rm d}V_{\rm gr}}{{\rm d}R}\right)_{R_{\rm c}}R^2_{\rm c}\right).$$ (96)

A comparison with the case of momentum coupling shows that for an identical radiative force profile, the critical point is located farther. With a self-consistent coupling, however, the radiative force will be different at each step of iterations (see below).

A numerical solution can be found using a shooting-splitting method described by Berruyer & Frisch (1983). We used this method in the case of optically thick media.

Two gradients (grains and gas velocities) are required. The gradient of the grains velocity in $R_{\rm c}$ can be deduced from Eq. (92):

$$\left(\frac{{\rm d}V_{\rm gr}}{{\rm d}R}\right)_{R_{\rm c}} ... ...1}{R^2_{\rm c}}\left(\beta_2 + \Phi_4 F_{\rm c}\right)\right].$$ (97)

Using this equation and Eq. (5), with the l'Hospital's rule we have the second degree equation for the gradient velocity:

$$\left(\frac{{\rm d}V}{{\rm d}R}\right)_{R_{\rm c}}^2 + B'(R_... ...rm d}R}\right)_{R_{\rm c}} + D'(R_{\rm c}, V_{\rm gr_c}) = 0,$$ (98)

where:

$$B'(R_{\rm c}, V_{\rm gr_c}) = - \frac{\beta_3}{2}\frac{H_{\rm c}}{R_{\rm c}^2 V_{\rm gr_c}},$$ (99)

$$D'(R_{\rm c}, V_{\rm gr_c}) = -\frac{1}{R_{\rm c}^2}\left( \frac{... ...}V_{\rm gr}}{{\rm d}R}\right)_{R_{\rm c}} -\frac{2F_{\rm c}}{R_{\rm c}}\right].$$ (100)

With a given $R_{\rm c}$, the gradient velocity of the gas at the critical point is:

$$\left(\frac{{\rm d}V}{{\rm d}R}\right)_{R_{\rm c}} = \frac{... ...V_{\rm gr_c}) + \sqrt{\Delta(R_{\rm c}, V_{\rm gr_c})}\right),$$ (101)

with:

$$\Delta(R_{\rm c}, V_{\rm gr_c}) =B^{'2}(R_{\rm c}, V_{\rm gr_c}) - 4 D'(R_{\rm c}, V_{\rm gr_c}) \ge 0.$$ (102)

At first, we look for that part of the critical line where $\Delta(R_{\rm c}, V_{\rm gr_c}) \ge 0$. For each $R_{\rm c}$ in this line, we derive $V_{\rm gr_c}$ from Eq. (81), $\left(\frac{{\rm d}V_{\rm gr}}{{\rm d}R}\right)_{R_{\rm c}}$ from Eq. (95) and $\left( \frac{{\rm d}V}{{\rm d}R}\right)_{R_{\rm c}}$ from Eq. (98). It is possible to perform integration inwards. The method of shooting-splitting enables us to find out the transsonic solution among these curves, for each possible choice of $R_{\rm c}$.

The shooting-splitting method is carried out until the relation (82) between the gas and the grain velocities is satisfied, with an accuracy roughly equal to that of the shooting-splitting method. Henceforth, we shall refer to the point where this limit is reached as $R_{\rm D}$.

When dynamics is coupled with the radiative transfer, the main difficulty encountered is the uncertainty of the existence of a critical point and a transsonic solution at each iteration. It is possible, for instance, to find a critical point which has a saddle configuration, but a negative gradient velocity for grains or gas in the subsonic part of the wind.

The Momentum Coupling Hypothesis can be used to determine an initial profile for the full problem and allow a preliminary research of models with solutions.

3.2 Solution of radiative transfer (iteration)

Details of the numerical procedure used to solve the radiative transfer equations developed in sect. 3.1.1 are described by Leung (1976). The special form of the linearized equations allows us to use an efficient algorithm devised by Rybicki (1971). Explanations on the radiative transport code can be found in Spagna & Leung (1983) and Egan et al. (1988).

Let us summarize the main steps of the solution. The differential equations described in Sect. 3.1.1 are solved numerically using the finite differences method. In particular, since the moment and energy balance equations are nonlinear, they are solved using an iterative scheme similar to the Newton-Raphson method, which includes the following steps:

1.

Using initial estimates of $T_{\rm gr}(r)$, $f^{\rm B}_{\nu}$, $f_{\nu}(r)$ and $\zeta_{\nu}(r)$, Eq. (30) is solved for $J_{\nu}(r)$ and $H_{\nu}(r)$ frequency by frequency;

2.

From $J_{\nu}(r)$, $H_{\nu}(r)$ and $T_{\rm gr}(r)$, $\eta_{\nu}$ is evaluated and then the ray Eq. (38) is solved for $u_{\nu}$, from which new estimates for $f^{\rm B}_{\nu}$, $f_{\nu}$ and $\zeta_{\nu}$ are obtained;

3.

With this set of values for $T_{\rm gr}(r)$, $J_{\nu}(r)$, $f^{\rm B}_{\nu}(r)$, $f_{\nu}(r)$ and  $\zeta_{\nu}(r)$ as starting solution, the system of moment and energy equations is solved to obtain the temperature correction  $\Delta T(r)$, with which  $T_{\rm gr}(r)$ is updated and the convergence and the energy conservation are tested.

This iterative process is repeated until convergence is achieved. For the initial estimates of $f^{\rm B}_{\nu}$, $f_{\nu}(r)$ and $\zeta_{\nu}(r)$ the diffusion approximation is used, i.e. $f^{\rm B}_{\nu}=\frac{1}{2}$, $f_{\nu}(r)=\frac{1}{3}$ and $\zeta_{\nu}(r)=1$. The initial temperature profile is approximated by a power law:

$$T_{\rm gr}(r)=\beta \left( \frac{r_{\ast}}{r}\right)^{\alpha} T_{\ast}$$ (103)

where $r_{\ast}$ and $T_{\ast }$ are the radius and effective temperature of the star. The values of the constants $\alpha$ and $\beta$ depend on the grains.

An improvement is obtained by starting with the profile of $T_{\rm gr}$ in an optically thin circumstellar envelope. For subsequent more opaque models, we simply use the temperature distribution calculated in the previous model as the starting solution. This provides an extremely smooth and efficient scheme.

3.3 Numerical scheme (iteration)

Coupling grain-gas dynamics with radiative transfer is difficult because of the different natures of the two physical problems. More specifically, dynamics is a local problem whereas radiative transfer is a global one. Thus, the solution of the radiative transfer uses an iterative scheme with an imposed radial grid, while the solution of dynamics requires the generation of a variable grid (with variable steps) together with the integration. A code was built in order to couple dynamics with radiative transfer through an iterative scheme. This not only allowed dynamical results to be used as an input in the transfer, but also, conversely, transfer results to be used as an input for dynamics. An initial radiative pressure profile is given to the dynamical code which provides a velocity profile, and then from this a number density profile of dust grains as an input of the radiative transfer code. Then we calculate opacities necessary to deduce a new radiative pressure profile, and so on, until convergence of the global iterative scheme is obtained.

This scheme is sketched in Fig. 1 which shows that, in order to couple the two codes, we calculate the density profile  $n_{\rm {gr}}(r)$, obtained in the radial grid of the dynamics by interpolation, into the radial grid of the transfer. We also want to obtain the profile of radiative force  $F_{\rm rad}$ at each new step in the grid of dynamics as results from the integration.

For these two calculations, we use the piecewise Cubic Splines numerical method (Mathews 1992). The principle is to construct a cubic function on each interval so that the resulting piecewise curve and its first and second derivatives are all continuous on the larger interval which includes the dimension of the circumstellar envelope. This method is used to "piece together'' the graphs of lower-degree polynomials. This allows for good accuracy in our calculations, with a no wiggling graph. Another advantage is to have one Spline function on each interval, and the values of the radiative force at each point of the envelope without having to define the radial grid of dynamics before the integration. The calculus of the coefficients S(i,j) (indicated in Fig. 1) of these cubic polynomes can be found in Mathews (1992).

Another numerical difficulty is encountered in the generation of the radiative transfer's fixed radial grid. Indeed, with most finite-difference methods, the succes of the transfer code's iterative scheme depends on the proper choice of spatial grid. As we are solving a boundary value problem, the boundary conditions must be satisfied accurately. In order that information be transmitted correctly at the interior points, it is necessary that near the boundaries, $\Delta \tau_{\nu}$, the optical depth in a layer of the dimension of the considered step be markedly less than one for the most opaque part of the spectrum. Away from the boundaries, no such constraint exists and larger  $\Delta \tau_{\nu}$ can be used. It is also essential that the radial increment $\Delta r$ be small enough to ensure uniform accuracy of the results. In regions where the source function varies rapidly, the number of radial steps must also be increased (Leung 1976).

Thus, at each iteration of the iterative system coupling dynamics with transfer, we must conduct numerical tests to check whether our previous choice of spatial grid was suitable for the new iteration.

	Figure 1: General numerical scheme concerning our three models (position coupling, momentum coupling and full problem) coupled with the transfer.
Open with DEXTER

4 Discussion

4.1 Comparison of Position Coupling with Momentum Coupling (optically thin case)

Table 1 indicates the input parameters for two models, the Position Coupling Hypothesis (PCH) and Momentum Coupling Hypothesis (MCH), for which only the mass loss rate ( $\dot{M_{\rm {gr}}}$) of the dust grains is different. The inner and the outer radii are expressed in units of the stellar radius. This radius is deduced from the stellar luminosity $L_{\ast }$ and the stellar effective temperature $T_{\ast }$ with the expression $L_{\ast}= 4 \pi r_{\ast}^2 \sigma T_{\ast}^4$ where $\sigma$ is the Stefan-Boltzmann constant. We use the stellar radius for the purpose of the discussion, but in order to improve numerical stability the numerical calculations for the dynamics use the inner radius R₀ as a unit. Table 2 lists some output physical quantities for both models (PCH and MCH).

With the same  $\dot{M_{\rm {gr}}}$ for PCH and MCH, it is impossible in both cases to find a solution where the limit $r_{\rm c}$ of the critical point  $r^n_{\rm c}$, if this limit exists, remains greater than r₀. This is a very good first example of the much marked differences between PCH and MCH.

   

Table 1: Input parameters for models under the Position Coupling Hypothesis (PCH) and the Momentum Coupling Hypothesis (MCH): $R_{\ast }$ is the stellar radius in units of the solar radius ( ${R_{\odot }}$), $L_{\ast }$ is the stellar luminosity in units of the solar luminosity ( ${L_{\odot }}$), $T_{\ast }$ the effective temperature of the star in Kelvins (K), $T_{\rm gas}$ the temperature of the isothermal gas (K), $C_{\rm s}$ the sound velocity deduced from $T_{\rm gas}$, R₀ the inner radius of the envelope, $R_{\rm f}$ the outer radius of the envelope, $\dot{M}$ the mass loss rate of the gas in units of a solar mass per year ( $M_{\odot }/$yr) and $\dot{M_{\rm {gr}}}$ the mass loss rate of the dust grains (same units).

INPUT

$R_{\ast}= 252.8~{R_{\odot}}$

$L_{\ast}= 4.6\times 10^{3}~{L_{\odot}}$

$T_{\ast}= 3000~{\rm K}$

$T_{\rm gaz}= 1300~{\rm K}$

$C_{\rm s}= 2.326~{\rm km\,s^{-1}}$

$R_0= 1.35\times 10^3~{R_{\odot}}=5.34~{R_{\ast}}$

$R_{\rm f}= 27\,10^3~{R_{\odot}}=106.80~{R_{\ast}}$

$\dot{M}= 10^{-7}~{M_{\odot}/{\rm yr}}$

PCH:

$\dot{M}_{\rm gr}= 10^{-11}~{M_{\odot}/{\rm yr}}$

$\dot{M} / \dot{M}_{\rm gr}= 10^4$

MCH:

$\dot{M}_{\rm gr}= 2.22\times 10^{-10}\,M_{\odot}/{\rm yr}$

$\dot{M} / \dot{M}_{\rm gr}\simeq 450$

	Figure 2: Velocity profiles for the grains ( $V_{\rm gr}$) and the gas (V) versus the radial distance (in units of the star radius: $R=r/r_{\ast }$).
Open with DEXTER

   

Table 2: A few output physical quantities for PCH and MCH. $\rho _0$ and $\rho _{\rm gr_0}$ (V₀ and $V_{\rm gr_0}$) are respectively the gas and dust densities (velocities) at the inner radius, $\rho _{\rm c}$ and $\rho _{\rm gr_c}$ the gas and dust densities at the critical point, $V_{\rm gr_c}$ is the grain velocity at the (gas) sonic point, $R_{\rm c}$ is the position of this sonic point and  $\tau _{\nu _{\rm ref}}$ the optical depth for suitable frequency reference ( $\nu _{\rm ref}=1.3\times 10^{15}\,{\rm Hz}$).

OUTPUT	PCH	MCH
$\rho_0 \,({\rm g/cm^3})$	$0.32\times 10^{-15}$	$0.33\times10^{-14}$
$\rho_{\rm gr_0}\,({\rm g/cm^{-3}})$	$0.3\times10^{-19}$	$0.96\times10^{-19}$
$\rho_{\rm c}\,({\rm g/cm^3})$	$0.14\times10^{-15}$	$0.14\times10^{-15}$
$\rho_{\rm gr_c}\,({\rm g/cm^{-3}})$	$0.14\times10^{-19}$	$0.14\times10^{-19}$
$V_0\,({C_{\rm s}})$	0.75	$7.40\times10^{-2}$
$V_{\rm f}\,({C_{\rm s}})$	3.21	5.93
$V_{\rm gr_0}\,({C_{\rm s}})$	0.75	5.60
$V_{\rm gr_f}\,({C_{\rm s}})$	3.21	57.02
$V_{\rm gr_c}\,({C_{\rm s}})$	1.00	22.00
$R_{\rm c}\,({R_{\ast}})$	7.00	7.00
$\tau _{\nu _{\rm ref}}$	0.090	0.107

Before discussing approximations for grain-gas dynamics, we can see that these two self-consistent models are optically thin (see Table 2, $\tau_{\nu_{\rm ref}}= 0.107$ for MCH and $\tau_{\nu_{\rm ref}}=0.090$ for PCH). Without the self-consistent models, it is impossible to make conclusions about the optical depth of the circumstellar envelope, precisely because the true optical depth is a result of the self-consistent models obtained when convergence occurs. We shall see hereafter that there are optically thick solutions which can be obtained by iterations in which, at the first step of iterations (starting from the R^-2 law for the radiative force) the solution of the transfer equations looks optically thin, and only the self-consistent model can predict the true optical depth.

First, let us note that in order to find the same limit $r_{\rm c}$ in both cases, the mass loss rate of dust grains ( $\dot{M_{\rm {gr}}}$) in PCH must be taken very low compared to the mass loss rate of the gas ($\dot{M}$); indeed, the mass loss rate ratio $\dot{M} / \dot{M}_{\rm gr}= 10^4$ is very high (whereas one would expect a few 10²). Furthermore, to estimate the minimum of this ratio in PCH (such that $R_{\rm c}=R_0$), we use the R^-2 law (64) for radiative force in an optically thin medium. For carbon grains, all the numerical computations show that the solution ${\cal F}_{\rm rad}(R) < F_{\rm rad_0}/R^2$, for each R. Since ${\cal F}_{\rm rad}(R)= F_{\rm rad_0}/R^2$, for $R_{\rm c}=R_0=1$, the critical radius Eq. (69) is:

$$R_{\rm C_{Park}} - \frac{\beta_1}{2}F_{\rm rad_0}=1.$$ (104)

Using the expression (101) and $\beta_1=\frac{\displaystyle \dot{M}_{\rm gr}}{\displaystyle \dot{M}} (\frac{\displaystyle r_0}{\displaystyle C_{\rm s}^2 m_{\rm gr}})$, we can infer the value of the maximum mass loss rates for the dust grains: $\dot{M}_{\rm gr}=1.015\times 10^{-11}~M_{\odot}/$yr. If $R_{\rm c}=R_0=R_{\ast}$, we find for the maximum $\dot{M}_{\rm gr}=1.083\times 10^{-11}~M_{\odot}/$yr. Thus it is impossible to obtain realistic ratios of mass loss rates with PCH, and the position $R_{\rm c}$ of the critical (sonic) point is very sensitive to small variations of $\dot{M_{\rm {gr}}}$.

The model using MCH enables us to explain this undersized value of $\dot{M_{\rm {gr}}}$ when we use PCH. To do this, we solve the self-consitent problem with the two assumptions MCH and PCH in such a way that the solutions give rise to the same sonic point (cf. Fig. 2). To obtain this result, we use the same set of input parameters (cf. Table 1), except  $\dot{M_{\rm {gr}}}$. $\dot{M_{\rm {gr}}}$ is found equal to $2.22\times 10^{-10}\,{\rm M_{\odot}/yr}$ with MCH, which correspond to a mass loss rate ratio of $\dot{M} / \dot{M}_{\rm gr}\simeq 450$ (cf. Table 1). This mass loss rate ratio is much more realistic, emphasizing that the drift velocity between the grains and the gas must be taken into account when dynamics is coupled with radiative transfer.

The difference between the two models (PCH and MCH) is obvious when comparing the sonic point Eq. (69) with Eqs. (84) and (86). Since, with our choice of parameters, both solutions are optically thin, the two radiative force profiles are very similar to the R^-2 law. In this case, the only difference between the two equations comes from the ratio of densities $\frac{\displaystyle n_{\rm gr}}{\displaystyle \rho}$, which corresponds exactly to  $\frac{\displaystyle 1}{\displaystyle m_{\rm gr}} \frac{\displaystyle \dot{M}_{\rm gr}}{\displaystyle \dot{M}}$ for PCH and to $\frac{\displaystyle 1}{\displaystyle m_{\rm gr}} \frac{\displaystyle \dot{M}_{... ...gr}}{\displaystyle \dot{M}} \frac{\displaystyle 1}{\displaystyle V_{\rm gr_c}}$ for MCH. Since we have the same critical point, the ratio of densities is also the same at $R=R_{\rm c}$, which implies different mass loss rate ratios. With the value of $V_{\rm gr_c}$ of table 2, $\dot{M}_{\rm gr}= 10^{-11} V_{\rm gr_c}$, which is the value of the mass loss rate of the grains in PCH.

Furthermore, if the two ratios of densities are identical at $R=R_{\rm c}$, they can not be equal at another point in the circumstellar envelope due to the increase of the drift velocity with R (see Fig. 2 where velocities of gas and grains are plotted). For the PCH model, this ratio is constant. This indicates that the density of grains relative to that of the gas varies with the radial distance, which suggests that drift velocity changes not only the results for velocity profiles but also all the physical conditions in the whole envelope. This has substantial consequences for thermodynamical conditions for instance, but also in the study of nucleation and growth of dust grains, or the propagation of shock waves.

At the outer limit of the envelope, the ratio of the terminal velocities of grains (noted  $V_{\rm gr_f}$ in Table 2) for MCH and PCH is comparable to the ratio of the mass loss rates of dust grains. This indicates that the number densities of grains in the two models become closer when the radial distance increases, contrarily to those of the gas, since the terminal velocity of the gas ($V_{\rm f}$) under MCH is lower than that of grains by a factor ten.

It is also interesting to make a further comparison. In Table 4, we replace the mass loss rate for dust under MCH ( $\dot{M}_{\rm gr}= 2.22\times 10^{-10}~M_{\odot}/$yr) with that previously chosen under PCH ( $\dot{M}_{\rm gr}= 10^{-11}~M_{\odot}/$yr). In this case, convergence is impossible because as early as at the first step in the MCH model, the optical depth is very high and unrealistic (about 5000). The critical point is close to that of the thermal wind model. In fact, if this model could converge, the position of the critical point would be farther than in PCH because of the decrease in radiative force for a thick medium. This confirms the large difference between the two approximations. With MCH, the medium will be more opaque.

It is, nevertheless, much easier with PCH than with MCH to explore coupling between dynamics and radiative transfer. With this assumption we can determine many fields of investigation and speculate about effects of dust grains on radiative transfer and mass loss. This interesting study will be the subject of our next paper. It will allow us to improve significantly the exploration of self-consistent models that take into account drift velocity.

   

Table 3: Input physical parameters for two models: One solved with MCH and one solved with FPB.

INPUT

$\dot{M}= 0.2625\times 10^{-7}~M_{\odot}/$yr

$\dot{M}_{\rm gr}= 0.675\times 10^{-10}~M_{\odot}/$yr

$\dot{M} / \dot{M}_{\rm gr}\simeq 450$

$R_{\ast}= 263.5~{R}_{\odot}$

$L_{\ast}= 5\,10^{3}~{L}_{\odot}$

$M_{\ast}=1.1~{M}_{\odot}$

$T_{\ast}= 3000~{\rm K}$

$T_{\rm gas}= 1300~{\rm K}$

$C_{\rm s}= 2.326~{\rm km\,s^{-1}}$

$R_0= 8.18~{R}_{\ast}$

$R_{\rm f}= 163.57~{R}_{\ast}$

   

Table 4: A few output physical quantities for MCH. With the R^-2 for the radiative pressure, the solution is optical thin. The corresponding parameters are in the column "thin''. With the self-consistent calculation, the solution is optically thick. The corresponding parameters are in the column "thick''.

OUTPUT	thin	thick
$\rho_0 \,({\rm g/cm^3})$	$2.42\times 10^{-13}$	$5.65\times 10^{-13}$
$\rho_{\rm gr_0}\,({\rm g/cm^{-3}})$	$7.06\times 10^{-19}$	$2.03\times 10^{-18}$
$\rho_{\rm c}\,({\rm g/cm^3})$	$1.36\times 10^{-18}$	$9.08\times 10^{-19}$
$\rho_{\rm gr_c} \,({\rm g/cm^3})$	$7.92\times 10^{-23}$	$6.57\times 10^{-23}$
$V_0\,({C_{\rm s}})$	$1.04\times 10^{-4}$	$4.46\times 10^{-5}$
$V_{\rm gr_0}\,({C_{\rm s}})$	$9.17\times 10^{-2}$	$3.2\times 10^{-2}$
$V_{\rm gr_c}\,({C_{\rm s}})$	44.30	35.54
$V_{\rm f}\,({C_{\rm s}})$	2.62	2.37
$V_{\rm gr_f}\,({C_{\rm s}})$	72.77	55.54
$R_{\rm c}\,({R_{\ast}})$	35	43
$\tau _{\nu _{\rm ref}}$	0.40	1.25

	Figure 3: Solutions of the MCH model. We note A the optically thin solutions, i.e. obtained with the R-2 for the radiative force. Solutions B are obtained with the self-consistent calculation. The radial distance is expressed in units of the stellar radius ( $R=r/r_{\ast }$).
Open with DEXTER

4.2 Momentum Coupling in optically thick media

Let us now continue looking at self-consistent models using MCH for which the solution is optically thick.

In Table 3 physical input parameters are listed for a new self-consistent model that will be compared to the model which ignores radiative transfer (where $F_{\rm rad}(R)=\frac{F_{\rm rad_0}}{R^2}$). This model corresponds to the first step of iterations. Hereafter we shall call it the "R^-2 model''.

Table 4 gives a few output physical quantities: the column "thin'' (respectively "thick'') displays results of the R^-2 model (respectively the self-consistent one). It is worth noting that the R^-2 model suggests that the medium may be optically thin since the optical depth $\tau_{\nu_{\rm ref}}=0.40 < 1$, but that after this first step the final solution of the self-consistent model is $\tau_{\nu_{\rm ref}}=1.25 > 1$. The circumstellar envelope is in fact optically thick. The R^-2 model, therefore, does not accurately approximate the self-consistent model. The differences between the two approaches are highly significant. Most of the results in Table 4 confirm this.

In Fig. 3 we compare the self-consistent optically thick solution (indexed by B) with the R^-2 model (indexed by A). Profiles of grain velocity  $V_{\rm gr}$, gas velocity V, dust number density  $n_{\rm gr}$ and radiative force  ${\cal F}_{\rm rad}$ are sketched. $n_{\rm gr}$ and  ${\cal F}_{\rm rad}$ are drawn solely in the subsonic part of the wind to highlight the zone where differences between models A and B for these two quantities are greater.

We see that the relative diminution of  ${\cal F}_{\rm rad}$ in model B compared to model A is amplified with increasing R in the innermost part of the envelope. In this region, $n_{\rm gr}$ of model B is higher than  $n_{\rm gr}$ of model A, especially at the inner boundary where, after coupling with radiative transfer, the density is multiplied by a factor of about 3. The density in the inner part of the wind is very sensitive to the variations in radiative force. The effect of dust grains on radiative pressure is an attenuation by absorption and scattering of the radiation. Thermal emission by the grains is not efficient because the maximum emission is in the infra-red spectral region, where the undimensioned cross section of radiative pressure  $Q_{\rm pr}(\nu)$ is low and also because the field is very isotropic. On the other hand, multiple scattering may be able to increase the radiative force, but our results show that absorption does inhibit such an effect.

Furthermore, at R=R₀, ${\cal F}_{\rm rad}$ is attenuated by the radiation coming from the inner layers of the dust envelope. In our model, we make a sharp separation between the thick medium, where dust grains exist, and the thin medium where we have just gas between the star and the circumstellar envelope. The radiation coming from the envelope at R=R₀, therefore, is no more isotropic and can play a significant role when cumulated with the effect of all grains in a sphere of radius larger than R₀ by about the mean free path. We may be inclined to expect that for an optically thicker medium this mean free path would be small and that the net Eddington flux at R=R₀ would decrease, but the number density of grains is actually greater and allows this increase of the radiative flux. In fact, the amount of inward flux close to R₀ depends on the number of grains rather than the size of the region emitting photons, all the more so as when R increases, temperature and thermal emission of the dust grains decrease.

Furthermore, beyond a few stellar radii in the envelope, and after a strong decrease, the number densities of both models are very close, but the influence of the inner part of the wind (in a zone lower than the subsonic region) remains effective over the whole envelope (cf. the velocities profiles). Indeed, grain and gas velocities are lower for solution B: the dust velocity profile is roughly about 1/3 below the solution A. Precisely the terminal velocities $V_{\rm gr_f}=72.77$ for solution A and $V_{\rm gr_f}=55.54$ for solution B.

Grains increase the wind velocity, driving the gas by friction, but they also attenuate the effect of radiative force through the effects of radiative transfer. This highlights two antagonistic mechanisms involving the dust grains.

Another effect on the wind structure appears in the self-consistent problem. The decrease in grain velocity between models A and B is greater than that of the gas, as is shown by the terminal velocities values for the two components (cf. Table 4). Let us remember that MCH consists of identifying radiative force and frictional force, so that frictional force decreases as soon as radiative force does so. Frictional force is proportional to the gas density and a function that increases with drift velocity. Because of the simultaneous decrease in gas velocity when the radiative force decreases, the density of the gas increases for a given position. Thus the drift velocity $V_{\rm d}$ must also decrease in a self-consistent model.

So when the medium becomes more opaque, not only do  $V_{\rm gr}$ and V decrease but also $V_{\rm d}$. This is why the effect of radiative transfer on dynamics is higher for the grain velocity than for that of the gas. This shows that values of the terminal velocities of gas in models A and B are very close.

   

Table 5: A few output physical quantities for the full problem (FPB) solved with the R^-2 law for the radiative force (column "thin'') and with the self-consitent model (column "thick'').

OUTPUT	thin	thick
$\rho_0 \,({\rm g/cm^3})$	$2.57\times 10^{-13}$	$6.91\times 10^{-13}$
$\rho_{\rm gr_0}\,({\rm g/cm^{-3}})$	$7.50\times 10^{-19}$	$2.65\times 10^{-18}$
$\rho_{\rm c}\,({\rm g/cm^3})$	$1.31\times 10^{-18}$	$8.10\times 10^{-19}$
$\rho_{\rm gr_c} \,({\rm g/cm^3})$	$7.76\times 10^{-23}$	$6.36\times 10^{-23}$
$V_0\,({C_{\rm s}})$	$9.78\times 10^{-5}$	$3.64\times 10^{-5}$
$V_{\rm gr_0}\,({C_{\rm s}})$	$8.63\times 10^{-2}$	$2.45\times 10^{-2}$
$V_{\rm gr_c}\,({C_{\rm s}})$	43.58	32.78
$V_{\rm f}\,({C_{\rm s}})$	2.59	2.29
$V_{\rm gr_f}\,({C_{\rm s}})$	69.28	49.11
$R_{\rm c}\,({R_{\ast}})$	35.79	45.55
$R_{\rm D} \,({R}_{\ast})$	17.58	20.45
$\tau _{\nu _{\rm ref}}$	0.42	1.69

	Figure 4: Velocity profiles of dust grain ( $V_{\rm gr}$) and gas (V) versus the radial distance (in units of the star radius: $R=r/r_{\ast }$). MCH and FPB are compared for the optically thin solutions (obtained with radiative force law in R-2) and for the optically thick solutions (obtained with the self-consistent models).
Open with DEXTER

4.3 Comparison of Momentum Coupling with the Full Problem

Table 5 is as Table 4, but for the full problem (FPB). We can note that $V_{\rm gr_0} \neq V_0$. In fact, the initial values of  $V_{\rm gr}$ and V are the same, but they rapidly differ through a boundary layer (see the paper by Berruyer & Frisch 1983). Since the radial dimension $\delta$ of this boundary layer is very small, r₀ and $r_0 + \delta$ are practically identical, when compared with other useful scales. Differences between the thin model and the thick model are amplified. The shift in the sonic points is larger in the full problem than in the case of MCH and the increase in the optical depth is greater still. The first conclusion we can draw, therefore, is that with a self-consistent treatment, dynamics requires the solution of the full problem.

In Fig. 4 this conclusion is still more obvious. When we compare thin models (MCH and FPB) with thick models, we see that the differences between velocities increase. When the medium is more opaque, the effects of the inertial force become marked compared with other terms of the momentum transfer for dust grains.

Furthermore, the radius $R_{\rm D}$ increases but the dimension $R_{\rm c}-R_{\rm D}$ of that region where the equation is not degenerated is greater and the inertial force is efficient over a greater part of the subsonic zone of the wind. The zone where MCH can be used also increases in size.

We can highlight the greater difference between the two models when the medium is more opaque by considering the double effects of the inertial force. At first, this force decreases the dust grains and gas velocity. Moreover, according to Eq. (4) (mass conservation for dust grains), when dust grains velocity decreases, their density increases. An increased density in radiative transfer leads to a lower radiative force and then to a still lower velocity. In the self-consistent model, both mechanisms due to the inertial force work to reduce the velocity.

Finally, the effect of radiative transfer on the drift velocity described for MCH is amplified under FPB conditions. A part of the momentum given to the grains by the radiative force is transferred to the gas while another part is used to overcome inertia and gravitationnal forces. For a given radiative force, we have a lower frictional force, and, since radiative force is also lower, velocity is lower, density is greater, and drift velocity is lower at each point of the envelope. Furthermore, since when MCH is relaxed both the gas and the grains have smaller velocities, and the sonic point is located further away from the star, FPB amplifies the faster diminution of grain velocity compared with gas velocity with MCH.

5 Conclusions

Whatever the choice for the treatment of dynamics, there are cases where the coupling with radiative transfer can produce fairly strong effects. The presented two-fluid, self-consitent models for dust driven winds show that the optical depth of the circumstellar envelope, even if very small, plays a highly significant role in dynamics. Depending on whether the medium is optically thick or not, solutions with different features are obtained. The roles of the drift velocity and the inertial force appear emphasized. In particular, the inclusion of drift velocities changes the mass ratio of dust grains to gas with the radial distance. Nonetheless, models using position coupling are very useful as tools for numerical experiments.

Wind structure, especially at the base of the envelope, is very sensitive to change in the radiative force profile, in particular through the number density of dust grains which determines the radiative force and the optical depth. When the medium becomes more opaque, the outwards flux is weakened and so the radiative force is weakened.

Grains accelerate the wind by friction with the gas but they can also weaken the radiative force, so that the presence of dust grains has two antagonistic effects on the wind structure.

The grain velocity is more sensitive than that of gas to the coupling with the radiative transfer. This effect is amplified under conditions of the full problem.

Including inertial force into the models decreases velocities and increases densities essentially in the inner part of the envelope. Those effects are amplified with the self-consistent model.

Acknowledgements

Part of this work has been performed using the computing facilities provided by the program "Simulations Interactives et Visualisation en Astronomie et Mécanique (SIVAM)''.

At each step of iterations, we determined the radiative force profile from the solution of the radiative transfer obtained using the code of Leung (1976) who can be thanked for leaving this code for public use.

The authors would like also to thank Dr. Rubens Freire for several helpful discussions.

References

 
• Auer, L. H. 1971, J. Quant. Spectrosc. Rad. Transf., 11, 573 In the text  
• Auer, L. H., & Mihalas, D. 1970, MNRAS, 149, 65+ In the text NASA ADS  
• Baines, M. J., Williams, I. P., & Asebiono, A. S. 1965, MNRAS, 130, 63 In the text  
• Batchelor, G. K. 1967, An Introduction to Fluid Dynamics (Cambridge University Press) In the text  
• Berruyer, N. 1991, A&A, 249, 181 In the text NASA ADS  
• Berruyer, N., & Frisch, H. 1983, A&A, 126, 269 NASA ADS  
• Biermann, P., & Harwitt, M. 1980, ApJ, 241, L105 In the text NASA ADS  
• Bulirsch, R., & Stoer, J. 1966, Numerische Mathematik, 8, 1 In the text  
• Cassinelli, J. P. 1979, ARA&A, 17, 275 In the text  
• Deguchi, S. 1980, ApJ, 236, 567 In the text NASA ADS  
• Draine, B. T. 1979, Ap&SS, 65, 313 In the text  
• Draine, B. T. 1985, ApJS, 57, 587 In the text NASA ADS  
• Egan, M. P., Leung, C. M., & Spagna, G. F. 1988, Computer Phys. Commun., 48, 271 In the text  
• Feautrier, P. 1964, C. R. Acad. Sci. Paris, 258, 3189 In the text  
• Gilman, R. C. 1972, ApJ, 178, 423 In the text NASA ADS  
• Hummer, D. G., & Rybicki, G. B. 1971, MNRAS, 152, 1+ In the text NASA ADS  
• Krüger, D., Gauger, A., & Sedlmayr, E. 1994, A&A, 290, 573 In the text NASA ADS  
• Krüger, D., & Sedlmayr, E. 1997, A&A, 321, 557 In the text NASA ADS  
• Kwok, S. 1975, ApJ, 198, 583 In the text NASA ADS  
• Lafon, J.-P. 1991, The Infrared Spectral Region of Stars, Proc. Int. Colloq. held in Montpellier (France), 16-19 october 1990, ed. C. Jaschek, & Y. Andrillat (Cambridge University Press), 193 In the text  
• Lafon, J. P. J., & Berruyer, N. 1991, A&ARv, 2, 249 NASA ADS  
• Leung, C. M. 1975, ApJ, 199, 340 In the text NASA ADS  
• Leung, C. M. 1976, J. Quant. Spectrosc. Rad. Transf., 16, 559 In the text NASA ADS  
• Lucy, L. B. 1971, ApJ, 163, 95+ In the text NASA ADS  
• Lucy, L. B. 1976, ApJ, 205, 482 In the text NASA ADS  
• Mathews, J. H. 1992, Numerical Methods for mathematics, science, and engineering (Prentice Hall International Editions) In the text  
• Menietti, J. D., & Fix, J. D. 1978, ApJ, 224, 961 In the text NASA ADS  
• Mihalas, D. 1978, Stellar Atmospheres (W. H. Freeman and Company, San Francisco) In the text  
• Netzer, N., & Elitzur, M. 1993, ApJ, 410, 701 In the text NASA ADS  
• Parker, E. N. 1958, ApJ, 128, 664+ In the text NASA ADS  
• Parker, E. N. 1960a, ApJ, 132, 175+ In the text NASA ADS  
• Parker, E. N. 1960b, ApJ, 132, 821+ In the text NASA ADS  
• Patzer, A. B. C., Gauger, A., & Sedlmayr, E. 1998, A&A, 337, 847 In the text NASA ADS  
• Phillips, J. P. 1979, A&A, 71, 115 In the text NASA ADS  
• Rybicki, G. B. 1971, J. Quant. Spectrosc. Rad. Transf., 11, 589 In the text  
• Salpeter, E. E. 1974a, ApJ, 193, 579 In the text NASA ADS  
• Salpeter, E. E. 1974b, ApJ, 193, 585 In the text NASA ADS  
• Scalo, J. M. 1977, A&A, 55, 253 In the text NASA ADS  
• Schaaf, S. A. 1963, Encyclopedia of Physics, Fluid Dynamics II, vol. VIII/2, ed. S. Flügge, & C. Truesdell (Springer-Verlag) In the text  
• Sedlmayr, E., & Dominik, C. 1995, Space Sci. Rev., 73, 211+ In the text NASA ADS  
• Simis, Y. J. W., Icke, V., & Dominik, C. 2001, A&A, 371, 205 In the text NASA ADS  
• Simpson, I. C., Simons, S., & Williams, I. P. 1980, Ap&SS, 71, 3 In the text  
• Spagna, G. F., & Leung, C. M. 1983, Computer Phys. Commun., 28, 337 In the text  
• Tabak, R., Hirth, J. P., Meyrick, G., & Roark, T. P. 1975, ApJ, 196, 457 In the text NASA ADS  
• Tielens, A. G. G. M. 1983, ApJ, 271, 702 In the text NASA ADS  
• Weymann, R. 1960, ApJ, 132, 380+ In the text NASA ADS  
• Winters, J. M., Le Bertre, T., Jeong, K. S., Helling, C., & Sedlmayr, E. 2000, A&A, 361, 641 In the text NASA ADS  
• Zuckerman, B. 1980, ARA&A, 18, 263 In the text NASA ADS