Open Access
Issue
A&A
Volume 667, November 2022
Article Number A75
Number of page(s) 16
Section Stellar structure and evolution
DOI https://doi.org/10.1051/0004-6361/202243540
Published online 09 November 2022

© L. Siess et al. 2022

Licence Creative CommonsOpen Access article, published by EDP Sciences, under the terms of the Creative Commons Attribution License (https://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

This article is published in open access under the Subscribe-to-Open model. Subscribe to A&A to support open access publication.

1. Introduction

During the thermally pulsing asymptotic giant branch (AGB) phase, low- and intermediate-mass stars develop strong mass loss and eject their envelope to become a short-lived (≲104 yr) post-AGB star before ending their life as a white dwarf.

Stellar winds are key ingredients that dictate the evolution of the star but also contribute to the chemical evolution of the Galaxy by releasing the by-products of stellar nucleosynthesis into the interstellar medium (ISM). These stellar winds manifest as nebulae surrounding their progenitor star. In binary systems the interaction of the wind of the evolved star with its companion leads to complex interactions, often resulting in a large-scale shaping of the nebula into specific morphologies such as spirals, discs, arcs or bipolar outflows (e.g., Mauron & Huggins 2006; Maercker et al. 2012; Kervella et al. 2016; Homan et al. 2020, 2021). These interactions can also affect the composition of the companion star as it potentially accretes a fraction of the wind material lost by the primary. It is commonly accepted that Carbon Enhanced Metal Poor (CEMP), CH, barium or yellow symbiotic stars to list only a few (for a review of chemically peculiar low-mass binaries, see Jorissen 2003) inherited their peculiar surface abundances from a the wind pollution from an AGB companion.

With the advent of high resolution imaging, we now have access to fine structural details of the circumstellar environments of AGB stars and their interpretation requires advanced modelling. The past years have seen a renewed interest in the study of the hydrodynamics of interacting winds in binary systems (e.g., Mohamed & Podsiadlowski 2012; van der Helm et al. 2019; de Val-Borro et al. 2017; Chen et al. 2018; El Mellah et al. 2020; MacLeod & Loeb 2020; Bermúdez-Bustamante et al. 2020; Maes et al. 2021; Malfait et al. 2021; Aydi & Mohamed 2022). However, in these studies the dust formation process was not taken into account. To our knowledge, dust nucleation has only been included in a single 3D hydrodynamics code (Mohamed & Podsiadlowski 2012), and grain growth is now treated in the CO5BOLD star-in-a-box models by Höfner & Freytag (2019).

The increase in opacity associated with dust formation is a key ingredient to provide the required acceleration that carries the wind beyond the escape velocity of the AGB star. Detailed 1D models of dust-driven AGB winds have been developed over the past decades with increasing refinement (e.g., Bowen 1988; Fleischer et al. 1992; Höfner et al. 1995, 2003; Andersen et al. 2003; Jeong et al. 2003; Höfner 2008; Mattsson et al. 2010; Bladh et al. 2015, 2019) and provided a clear picture of the phenomenon. It is now widely accepted that mass ejection proceeds in two steps. Stellar surface pulsations actuated by large-scale convective motions of the stellar interior (Freytag et al. 2017) provide the necessary kinetic energy to lift the stellar matter up to regions favourable for dust formation. The opaque dust material absorbs the stellar photons and is radially accelerated outward. The collision of the grains with the surrounding gas drags it along, until it reaches the stellar escape velocity, resulting in a large-scale outflow freed from the gravitational attraction of the AGB star. For this scenario to be successful, condensation into the solid state should be fast enough to ensure that the gas does not fall back onto the star, requiring efficient cooling mechanisms to pull the condensation radius closer to the stellar surface. This classical picture also assumes that gas and dust are strongly coupled which remains a valid assumption for high mass loss rates (e.g., Sandin & Höfner 2004).

In this paper, we present the implementation of mass ejection by dust-driven stellar winds and the treatment of dust formation including equilibrium chemistry, into the smoothed particle hydrodynamics code PHANTOM. The chemical reactions governing molecule formation as well as the type of dust that forms strongly depend on the C to O abundance ratio (C/O). This work focuses on carbon-rich mixtures (C/O > 1), whose chemistry is better known and has been extensively studied in the context of late-type stars (see, e.g., Gail et al. 1984; Gail & Sedlmayr 1985, 1987a, 1988; Dorfi & Höfner 1991).

This paper is the first in a series where a set of critical ingredients for the modelling of companion-perturbed dusty stellar winds will be implemented in the smoothed particle hydrodynamics code PHANTOM. The subsequent papers will revolve around the following topics: the treatment of gas-dust drift using the “one-fluid approximation” of hydrodynamics for small grains available in the code (Laibe & Price 2014b; Price & Laibe 2015), the radiative transfer and the acceleration of the dust-component of the wind using by the attenuated stellar radiation field, the incorporation of realistic cooling and heating rate of the gas and dust mixture, and a pulsating stellar atmosphere.

The current paper is organized as follows: in Sect. 2 we present the modelling of dust nucleation and dust growth and in Sect. 3 the physics involved in the SPH hydrodynamics of dusty winds. The numerical treatment of particle ejection in PHANTOM is described Sect. 4 and the validation of our approach is presented in Sect. 5 by comparing our 3D simulations with analytic profiles under various physical conditions and setups. We then discuss future prospects is Sect. 6 and conclude in Sect. 7.

2. Dust creation

The physical ingredients and methods used to simulate dust-driven winds closely follows the description presented in the reference textbook of Gail & Sedlmayr (2014). We consider a single-fluid approach where the grains are at rest with the gas particles and focus on mass loss from carbon-rich stars (C/O > 1). For a consistent calculation of the fluid dynamics, the abundance and opacity of the dust grains has to be determined at every point in space. The chemistry of the wind and the calculation of the dust properties are presented in Sects. 2.2 and 2.3. Knowledge of these quantities will then allow us to calculate the radiative acceleration as described in Sect. 3.

2.1. Grain growth and destruction

Classically, dust formation proceeds in two steps. First, the formation of critical clusters which serve as condensation seeds and then the growth of dust grains to macroscopic sizes by the addition of small entities called monomers. The nucleation process accounts for the formation of the seeds and involves species in chemical equilibrium. Such an equilibrium prevails as long as the clusters are small that is if their size N ≤ N0. Grain growth thus proceeds through the addition of small clusters of size N ≤ N0. In our calculations of C-based dust, the growth of graphite is assumed to proceed homogeneously by the addition of clusters of the same type (Ci-molecules where we consider i = 1, 2)

C N + C i C N + i , $$ \begin{aligned} \mathrm{C} _\mathrm{N} + \mathrm{C} _i \rightarrow \mathrm{C} _{\mathrm{N} +i}, \end{aligned} $$(1)

and by the reactions involving C-based molecules of a different type (heteromolecular growth)

C N + C 2 H 2 C N + 2 + H 2 , $$ \begin{aligned}&\mathrm{C} _\mathrm{N} + \mathrm{C} _2\mathrm{H} _2 \rightarrow \mathrm{C} _{\mathrm{N} +2} + \mathrm{H} _2, \end{aligned} $$(2)

C N + C 2 H C N + 2 + H . $$ \begin{aligned}&\mathrm{C} _\mathrm{N} + \mathrm{C} _2\mathrm{H} \ \rightarrow \mathrm{C} _{\mathrm{N} +2} + \mathrm{H} . \end{aligned} $$(3)

This set of reactions corresponds to N0 = 4. In the following presentation of the theory, for the calculation of the density of dust grains, it is necessary to define the minimum size Nl of a dust grains. This quantity is not easy to determine and following Gail & Sedlmayr (1988) we will use Nl = 1000.

For a growth process to be energetically favourable, the cluster size must exceed a critical value N*. The equilibrium abundance of clusters smaller than N* decreases with increasing N, they are thus unstable and do not participate to the dust production. Conversely, once the seed has reached the critical size, growth can proceed. It results that the cluster with size N* is the least abundant and because it has the slowest reaction rate, it will set the nucleation rate. The existence of this critical size is the reason why grain formation is a two step process. The formation of critical clusters is considered as a stationary process while the growth process has to be treated as a time-dependent problem.

The formation of C-based dust grains depends on the amount of atomic carbon and carbon-containing molecules available in the gas phase. We present the method for calculating these atomic and molecular abundances in the next section. The grain formation and growth is treated with the method of moments and described in Sect. 2.3.

2.2. Molecular abundances

Apart from atomic carbon, only the molecules C2, C2H, and C2H2 were considered to contribute to the formation of carbonaceous dust (Gail et al. 1984). However, the formation of these molecules competes with the formation of other molecules containing carbon, e. g., HCN and CH4 (Gail & Sedlmayr 2014, ch. 10.4), whose formation in turn competes with the formation of all molecules containing H, C, and N. Thus, for a consistent determination of the abundances of molecules contributing to the formation of C-based dust, a network of chemical reactions must be considered. In our code, we take into account 26 molecules H2, OH, H2O, CO, CO2, CH4, C2H, C2H2, N2, NH3, CN, HCN, Si2, Si3, SiO, Si2C, SiH4, S2, HS, H2S, SiS, SiH, TiO, TiS, TiO2 and TiS and 26 chemical reactions accounting for their formation (details about the chemical network is provided in Appendix A). We assume that all these reactions are in chemical equilibrium. For example, to calculate the abundance of H2, the law of mass action of the reaction 2 H → H2 in chemical equilibrium is used, which reads

P H 2 = P H 2 K H 2 , $$ \begin{aligned} P_{\mathrm{H} _{2}} = P_\mathrm{H} ^2 K_{\mathrm{H} _{2}}, \end{aligned} $$(4)

(cf. Eq. (10.16) of Gail & Sedlmayr 2014), where PH2 and PH are the partial pressures of H2 and H, respectively and KH2 is the dissociation constant of the reaction. The latter is given by

K H 2 = e Δ G H 2 / ( R T g ) , $$ \begin{aligned} K_{\mathrm{H} _2} = \mathrm{e} ^{-\Delta G_{\mathrm{H} _2} / (\mathcal{R} T_{\mathrm{g} })}, \end{aligned} $$(5)

where ℛ is the universal gas constant, Tg the gas temperature, and ΔGH2 is the change in Gibbs energy of the reaction,

Δ G H 2 = Δ G f ( H 2 ) 2 Δ G f ( H ) . $$ \begin{aligned} \Delta G_{\mathrm{H} _2} = \Delta G_\mathrm{f} (\mathrm{H} _2) - 2 \Delta G_\mathrm{f} (\mathrm{H} ). \end{aligned} $$(6)

We take the values of the Gibbs energies of formation ΔGf for all atoms and molecules from the JANAF tables (Chase 1986).

By use of relations of the form of Eq. (4), the partial pressures of the different molecules can be expressed as a function of the partial pressures of the atomic species, so that, e. g., the total number of hydrogen atoms per unit volume n⟨H⟩, atomic and bound in molecules, can be expressed as

n H k T g = P H + 2 P H 2 K H 2 , $$ \begin{aligned} n_{\langle \mathrm{H} \rangle } k T_{\mathrm{g} }= P_\mathrm{H} + 2 P_\mathrm{H} ^2 K_{\mathrm{H} _2}, \end{aligned} $$(7)

assuming that H2 is the dominant hydrogenous molecule. Similarly, the total number density of nitrogen atoms n⟨N⟩ = ϵNn⟨H⟩, where ϵN is the nitrogen abundance by number relative to hydrogen, reads

n N k T g = P N + P N P H 3 K NH 3 + P N P C K CN + P N P C P H K HCN + 2 P N 2 K N 2 , $$ \begin{aligned}&n_{\langle \mathrm{N} \rangle } k T_{\mathrm{g} }= P_\mathrm{N} + P_\mathrm{N} P_\mathrm{H} ^3 K_{\mathrm{NH} _3} + P_\mathrm{N} P_\mathrm{C} K_\mathrm{CN} \nonumber \\&\qquad \qquad + P_\mathrm{N} P_\mathrm{C} P_\mathrm{H} K_\mathrm{HCN} + 2 P_\mathrm{N} ^2 K_{\mathrm{N} _2}, \end{aligned} $$(8)

where KNH3, KHCN, and KN2 are the dissociation constants of the reactions of formation of NH3, HCN, and N2, respectively. The constants are calculated by use of Eq. (5), substituting ΔGH2 by the change in Gibbs energy of the respective reaction, which is calculated according to relations of the form of Eq. (6).

Equations (7) and (8), along with the equations for the total number densities of the other considered elements, n⟨C⟩, n⟨O⟩, etc. (see Table 1 for the list of considered species), form a system of equations whose solution gives the atomic partial pressures Pel of the elements “el”. These atomic partial pressures are then used to calculate the molecular partial pressures Pmol of the molecules “mol” or, equivalently, their number densities nmol = Pmol/(kTg), k being the Boltzmann constant, by means of Eq. (4) and similar equations for the other molecules. These equations are equivalent to the mass conservation equations for the species involved in the nucleation phase.

Table 1.

Elemental abundances by number relative to hydrogen ϵ used in our models (Ferrarotti & Gail 2002).

2.3. Moment equations

For the simulation of grain growth and destruction, we use the moment equations described by Gauger et al. (1990; see also Gail & Sedlmayr 1988; Fleischer et al. 1992; Gail & Sedlmayr 2014). This method does not require the calculation of the grain size distribution f(N, t) itself, but uses its moments. For spherical grains, the moment of order i of the grain size distribution reads

K i = N = N l N i / 3 f ( N , t ) , $$ \begin{aligned} \mathcal{K} _i = \sum _{N = N_l}^\infty N^{i / 3} f(N, t), \end{aligned} $$(9)

where N is the grain size expressing the number of monomers in a grain, f(N, t) the number density of grains of size N and Nl ∼ 1000 the limit below which the cluster is not considered as a dust grain. The surface area of a grain composed of N monomers can be estimated as

A N = 4 π a 0 2 N 2 / 3 , $$ \begin{aligned} \mathcal{A} _N = 4 \pi a_0^2\,N^{2 / 3} , \end{aligned} $$(10)

where a0 is the radius of a monomer inside a grain, which is derived from the density ρcarb of the grain material, i.e.

a 0 = [ 3 A carb m u / ( 4 π ρ carb ) ] 1 / 3 , $$ \begin{aligned} a_0 = [3 A_\mathrm{carb} m_\mathrm{u} / (4 \pi \rho _\mathrm{carb} )]^{1 / 3}, \end{aligned} $$(11)

where Acarb = 12.011 denotes the atomic weight of carbon in our case. For a density ρcarb = 2.25 g cm−3 (Smithsonian Physical Tables, Forsythe 2003) we find a0 = 1.28 × 10−8 cm. Global dust properties can be estimated from the knowledge of the moments. In particular,

  • a)

    the mean grain density nd = 𝒦0, that is the density of all dust grains with size N ≥ Nl.

  • b)

    the average grain radius ⟨a⟩=a0𝒦1/𝒦0.

  • c)

    the average grain surface A N = 4 π a 0 2 K 2 / K 0 $ \langle \mathcal{A}_N\rangle= 4 \pi a_0^2 \mathcal{K}_2/\mathcal{K}_0 $.

  • d)

    the average particle size ⟨N⟩=𝒦3/𝒦0.

  • e)

    the number density of monomers of size N ≥ Nl condensed into grains ncond = 𝒦3.

As demonstrated in Gauger et al. (1990), an elegant way of writing the moment equations is to normalize the moments by n⟨H⟩, the number of H atoms per unit volume

K ^ i = K i n H = K i m ¯ H ρ , $$ \begin{aligned} \widehat{\mathcal{K} }_i = \frac{\mathcal{K} _i}{n_{\langle \mathrm{H} \rangle }}= \frac{\mathcal{K} _i \, \bar{m}_\mathrm{H} }{\rho }, \end{aligned} $$(12)

where ρ is the density and m ¯ H = m H i A i ϵ i $ \bar{m}_{\mathrm{H}}= m_{\mathrm{H}} \sum_i A_i \epsilon_i $ the mean mass per H-atom calculated using the abundance by number relative to hydrogen ϵi, atomic weight Ai of species i (see Table 1 for a list of the selected species) and mH the mass of an hydrogen atom. An important consequence of the previous relation is that at each time, the abundance of free carbon atoms (relative to H) is known and given by

ϵ C = ϵ C 0 K ^ 3 , $$ \begin{aligned} \epsilon _\mathrm{C} = \epsilon _\mathrm{C} ^0-\widehat{\mathcal{K} }_3, \end{aligned} $$(13)

where ϵ C 0 $ \epsilon_\mathrm{C}^0 $ is the initial abundance at the beginning of the nucleation. With these normalized variables, the evolution of the first four moments reads

d J ^ d t = J ^ s J ^ τ , $$ \begin{aligned}&\frac{\mathrm{d} \widehat{J}_{*}}{\mathrm{d} t} = \frac{\widehat{J}_{*}^s - \widehat{J}_{*}}{\tau _{*}}, \end{aligned} $$(14)

d K ^ 0 d t = J ^ , $$ \begin{aligned}&\frac{\mathrm{d} \widehat{\mathcal{K} }_0}{\mathrm{d} t} = \widehat{J}_{*}, \end{aligned} $$(15)

d K ^ i d t = i K ^ i 1 3 τ + N l i / 3 J ^ , $$ \begin{aligned}&\frac{\mathrm{d} \widehat{\mathcal{K} }_i}{\mathrm{d} t} = \frac{i\, \widehat{\mathcal{K} }_{i-1}}{3 \tau } +N_l^{i/3}\widehat{J}_{*}, \end{aligned} $$(16)

where d/dt is the material derivative. τ−1 is the net rate of growth or evaporation of the grains (see next section), J ^ $ \widehat{J}_{*} $ is the normalized rate of formation of critical clusters per unit volume, J ^ s $ \widehat{J}_{*}^{\mathrm{s}} $ the normalized quasi-stationary rate, and τ* the relaxation time towards equilibrium (see Sect. 2.3.2). The last term in Eq. (16) describes the transition from the (microscopic) newly formed clusters to the (macroscopic) grains with size N >  Nl. As shown in Höfner & Dorfi (1992), the inclusion of this term has a small impact. The solutions of Eqs. (14)–(16) are computed analytically assuming that during the integration time step Δt, the quasi-stationary rate J ^ s $ \widehat{J}_{*}^{\mathrm{s}} $ is constant. We stress that Δt is constrained so this assumption holds.

2.3.1. Net rate of grain growth and destruction

For our C-rich chemistry, as described by Gauger et al. (1990), the net rate of grain growth and destruction can be expressed as

1 τ = A 1 π { α 1 P C m C v rel , C ( 1 1 S α ( C ) b C ) + 2 α 2 P C 2 m C 2 v rel , C 2 ( 1 1 S 2 α ( C 2 ) b C 2 ) + 2 α 2 , 1 c P C 2 H m C 2 H v rel , C 2 H ( 1 1 S 2 α c ( C 2 H ) b C 2 H c ) + 2 α 2 , 2 c P C 2 H 2 m C 2 H 2 v rel , C 2 H 2 ( 1 1 S 2 α c ( C 2 H 2 ) b C 2 H 2 c ) } , $$ \begin{aligned}&\frac{1}{\tau } = \frac{\mathcal{A} _1}{\sqrt{\pi }} \Biggl \{\alpha _1 \frac{P_\mathrm{C} }{m_\mathrm{C} \, {v}_{\mathrm{rel} ,\mathrm{C} }} \left(1 - \frac{1}{S}\frac{\alpha _{*}(\mathrm{C} )}{b_\mathrm{C} } \right) \nonumber \\&\qquad + 2 \alpha _2\frac{P_{\mathrm{C} _2}}{m_{\mathrm{C} _2}\, {v}_{\mathrm{rel} ,\mathrm{C} _2}} \left(1-\frac{1}{S^2} \frac{\alpha _{*}(\mathrm{C} _2)}{b_{\mathrm{C} _2}}\right) \nonumber \\&\qquad + 2 \alpha ^c_{2,1} \frac{P_{\mathrm{C} _2\mathrm{H} }}{m_{\mathrm{C} _2\mathrm{H} } \, {v}_{\mathrm{rel} ,\mathrm{C} _2\mathrm{H} }}\left(1-\frac{1}{S^2}\frac{\alpha ^c_{*}(\mathrm{C} _2\mathrm{H} )}{b^c_{\mathrm{C} _2\mathrm{H} }} \right) \nonumber \\&\qquad +2 \alpha ^c_{2,2} \frac{P_{\mathrm{C} _2\mathrm{H} _2}}{ m_{\mathrm{C} _2\mathrm{H} _2}\, {v}_{\mathrm{rel} ,\mathrm{C} _2\mathrm{H} _2}} \left(1-\frac{1}{S^2} \frac{\alpha ^c_{*}(\mathrm{C} _2\mathrm{H} _2)}{b^c_{\mathrm{C} _2\mathrm{H} _2}} \right) \Biggr \} , \end{aligned} $$(17)

where we only consider the chemical reactions given by Eqs. (1)–(3) involving C, C2, C2H and C2H2. The average relative velocity vrel, i between the reacting species i (mass mi) and the dust particle, in the absence of drift, is given by

v rel , i = 2 k T g m i . $$ \begin{aligned} {v}_{\mathrm{rel} ,i} = \sqrt{\frac{2 kT_{\mathrm{g} }}{m_i}} . \end{aligned} $$(18)

The quantities α1 = 0.37 and α2 = 0.34 are the sticking coefficients of C and C2 (Thorn & Winslow 1957) and α 2 , i = 1 , 2 c $ \alpha^c_{2,i=1,2} $ the reaction efficiency of the considered molecules (C2H or C2H2) with carbon. Following Gail & Sedlmayr (1988), we set α 2 , 1 c = α 2 , 2 c = α 2 $ \alpha^c_{2,1}=\alpha^c_{2,2}=\alpha_2 $. The coefficients α*(i) quantify the deviation from a thermal equilibrium population of states and are given by α*(i)=(Td/Tg)1/2, α c = 1 $ \alpha^{c}_{*} = 1 $ where Td and Tg refer to the dust and gas temperatures, respectively. The quantity bi accounts for departures of the particle densities from chemical and thermodynamical equilibrium between the solid (dust) and gas phase. Assuming chemical equilibrium, the coefficients bi are given by

b C = T d T g , $$ \begin{aligned}&b_\mathrm{C} = \frac{T_{\mathrm{d} }}{T_{\mathrm{g} }}, \end{aligned} $$(19)

b C 2 = K C 2 ( T g ) K C 2 ( T d ) T d T g , $$ \begin{aligned}&b_{\mathrm{C} _2} = \frac{K_{\mathrm{C} _2}(T_{\mathrm{g} })}{K_{\mathrm{C} _2}(T_{\mathrm{d} })}\frac{T_{\mathrm{d} }}{T_{\mathrm{g} }}, \end{aligned} $$(20)

b C 2 H c = K C 2 H ( T g ) K C 2 H ( T d ) , $$ \begin{aligned}&b_{\mathrm{C} _2\mathrm{H} }^{c} = \frac{K_{\mathrm{C} _2\mathrm{H} }(T_{\mathrm{g} })}{K_{\mathrm{C} _2\mathrm{H} }(T_{\mathrm{d} })}, \end{aligned} $$(21)

b C 2 H 2 c = K C 2 H 2 ( T g ) K H 2 ( T d ) K C 2 H 2 ( T d ) K H 2 ( T g ) . $$ \begin{aligned}&b_{\mathrm{C} _2\mathrm{H} _2}^{c} = \frac{K_{\mathrm{C} _2\mathrm{H} _2}(T_{\mathrm{g} }) K_{\mathrm{H} _2}(T_{\mathrm{d} })}{K_{\mathrm{C} _2\mathrm{H} _2}(T_{\mathrm{d} }) K_{\mathrm{H} _2}(T_{\mathrm{g} })}. \end{aligned} $$(22)

The supersaturation ratio S is defined as the ratio of the partial pressure of carbon in the gas phase divided by the vapour pressure of solid carbon

S = P C ( T g ) / P C , solid v ( T d ) , $$ {S}= P_\mathrm{C} (T_{\mathrm{g} }) / P^\mathrm{v} _{\mathrm{C,solid}} (T_{\mathrm{d} }), $$(23)

where the saturation vapour pressure of carbon in the solid phase

P C , solid v = K C , solid 1 = e Δ G C , solid / ( R T ) , $$ P^\mathrm{v} _\mathrm{C,solid} = K_\mathrm{C,solid}^{-1} = \mathrm{e} ^{\Delta G_\mathrm{C,solid} / (\mathcal{R} T)}, $$(24)

is derived using the dust temperature (T = Td). The free enthalpies of formation are taken from the JANAF tables and according to Eq. (6) reads ΔGC, solid = ΔGf(C, solid)−ΔGf(C). See Appendix B for a more detailed explanation of Eq. (17). When gas and dust are thermally coupled (Tg = Td), which is the assumption adopted in the current paper, the expressions for the growth rate simplifies since in this case α*(i)=bi = 1. The full equations are however given so when radiative transfer is implemented, the expressions remain valid. Deviations from chemical and/or thermal equilibrium may be important in 3D simulations because of the development of structures (spirals, shocks, ...) that affect both the chemistry and dust temperature and also because of the dependence of the grain growth rate on the gas-dust drift velocity (see e.g. Krüger & Sedlmayr 1997; Sandin & Höfner 2004).

2.3.2. Formation rate of seed particles

The quasi stationary rate of formation of critical clusters, i.e., clusters of size N*, per unit volume reads (chap. 13.7.4 of Gail & Sedlmayr 2014; Gail et al. 1984)

J s = β A N Z n ˚ d ( N ) . $$ \begin{aligned} J_{*}^{\mathrm{s} } = \beta \mathcal{A} _{N_{*}} Z\, \mathring{n}_{\mathrm{d} }(N_{*}). \end{aligned} $$(25)

The quantity n ˚ d ( N ) $ \mathring{n}_{\mathrm{d}}(N) $ is the density of grains of size N in thermal equilibrium,

n ˚ d ( N ) = n C exp { Δ F k T g } = n C exp { ( N 1 ) ln S ~ θ N ( N 1 ) 2 / 3 T g } , $$ \begin{aligned} \mathring{n}_{\mathrm{d} }(N) = n_{\mathrm{C} } \exp \left\{ -\frac{\Delta F}{kT_{\mathrm{g} }}\right\} = n_{\mathrm{C} } \exp {\left\{ (N -1) \ln {\tilde{{S}}} - \frac{\theta _{N} (N - 1)^{2 / 3}}{T_{\mathrm{g} }}\right\} }, \end{aligned} $$(26)

where ΔF is the free energy of formation of the grain from the vapour (Draine & Salpeter 1977). The saturation ratio of the gas phase S ~ $ \tilde{{{S}}} $ is computed using T = Tg in Eq. (24), nC is the density of monomers, i.e., of C atoms, and the quantity θN is related to the surface tension σ of the carbon grain material,

θ N = σ A 1 k { 1 + [ N h / ( N 1 ) ] 1 / 3 } = θ 1 + ( N h / ( N 1 ) ) 1 / 3 , $$ \begin{aligned} \theta _{N} = \frac{\sigma \, \mathcal{A} _{1} }{k\left\{ 1 + \left[N_{h} / (N - 1)\right]^{1 / 3}\right\} } = \frac{\theta _{\infty }}{1+(N_{h}/(N-1))^{1/3}}, \end{aligned} $$(27)

where we use a constant surface tension of σ = 1400 erg cm−2 for carbon grains (Tabak et al. 1975). The surface of a monomer (N = 1) is given by 𝒜1 and Nh is the particle size for which the surface tension is reduced to one half of its value. Following Gail et al. (1984) (see also ch. 13.7 of Gail & Sedlmayr (2014) for a longer discussion), we set Nh = 5.

The size of the critical cluster N* is determined by the root of the equation n ˚ d / N = 0 $ \partial \mathring{n}_{\mathrm{d}} / \partial N = 0 $. Its value is given by

N = 1 + N , 8 { 1 + [ 1 + 2 ( N h N , ) 1 / 3 ] 1 / 2 2 ( N h N , ) 1 / 3 } 3 , $$ \begin{aligned} N_{*} = 1+\frac{N_{*,\infty }}{8}\left\{ 1+\left[1+2\left(\frac{N_{h}}{N_{*,\infty }}\right)^{1/3}\right]^{1/2}-2\left(\frac{N_{h}}{N_{*,\infty }}\right)^{1/3}\right\} ^{3}, \end{aligned} $$(28)

where

N , = ( 2 θ / 3 T g ln S ~ ) 3 . $$ \begin{aligned} N_{*,\infty } = (2\theta _\infty /3T_{\mathrm{g} }\ln \tilde{{S}})^{3} . \end{aligned} $$(29)

The Zeldovich factor Z entering Eq. (25) measures the deviation of the grain distribution from equilibrium and is defined by

Z = 1 2 π 2 ln n ˚ d N 2 | N $$ \begin{aligned} Z = \sqrt{\frac{1}{2 \pi } \left.\frac{\partial ^{2} \ln {\mathring{n}_{\mathrm{d} }}}{\partial N^{2}}\right|{_{N_{*}}}} \end{aligned} $$(30)

and β, which gives the flux of monomers colliding with the grains, by

β = 1 2 π k T g [ α 1 P C m C + 4 α 2 ( P C 2 m C 2 + P C 2 H m C 2 H + P C 2 H 2 m C 2 H 2 ) ] , $$ \begin{aligned} \tiny \beta = \frac{1}{\sqrt{2\pi k T_{\mathrm{g} }}} \left[\alpha _\mathrm{1} \frac{P_\mathrm{C} }{\sqrt{m_\mathrm{C} }} + 4 \alpha _2 \left(\frac{P_{\mathrm{C} _2}}{\sqrt{m_{\mathrm{C} _2}}} + \frac{P_{\mathrm{C} _2\mathrm{H} }}{\sqrt{m_{\mathrm{C} _2\mathrm{H} }}} + \frac{P_{\mathrm{C} _2\mathrm{H} _2}}{\sqrt{m_{\mathrm{C} _2\mathrm{H} _2}}}\right) \right] , \end{aligned} $$(31)

(Gail et al. 1984). With these quantities, the relaxation time towards equilibrium writes

τ = ( 2 π Z 2 β A N ) 1 . $$ \begin{aligned} \tau _{*} = \left(2 \pi Z^{2} \beta \mathcal{A} _{N_{*}}\right)^{-1}. \end{aligned} $$(32)

As an illustration, we show in Fig. 1 the nucleation rate calculated with Eq. (25) as a function of temperature and pressure for a C/O ratio of 10. We use the chemical mixture of Ferrarotti & Gail (2002) which is very similar to that used by Gail et al. (1984). The computed nucleation rates are nearly identical to the results presented in Fig. 5 of Gail et al. (1984). Again, we find that nucleation mainly occurs around the carbon grain condensation temperature near 1500 K where the most abundant molecules are H2, C2H2, and CO. With increasing temperature, the dust starts to be dissociated and the nucleation rate drops. In the regime of low temperature on the other hand, dust formation is almost prohibited because the growth timescale becomes exceedingly long due to the very low density.

thumbnail Fig. 1.

Nucleation rate per hydrogen atom ( J s / n H $ J_{*}^{\mathrm{s}}/n_{\langle\mathrm{H}\rangle} $) as a function of temperature and pressure for a C/O ratio of 10. The computed values are nearly identical to the results presented in Fig. 5 of Gail et al. (1984).

3. The physics of dust-driven winds

3.1. 3D hydrodynamics

For our 3D hydrodynamic simulations we use the parallel smoothed particle hydrodynamics (SPH) code PHANTOM (Price & Federrath 2010; Lodato & Price 2010; Price et al. 2018). In the framework of SPH, the hydrodynamic equations for a discretised fluid surrounding a central gravitational potential read (see also the review by Price 2012)

ρ i = j m j W ( | r i r j | , h i ) , $$ \begin{aligned}&\rho _i = \sum _j m_j W(|\boldsymbol{r}_i - \boldsymbol{r}_j|, h_i), \end{aligned} $$(33)

D v i D t = j m j [ P i + q i Ω i ρ i 2 i W ij ( | r i r j | , h i ) + P j + q j Ω j ρ j 2 i W ij ( | r i r j | , h j ) ] G M r i 2 ( 1 Γ ) , $$ \begin{aligned}&\frac{\mathrm{D} \boldsymbol{v}_i}{\mathrm{D} t} = - \sum _j m_j \left[\frac{P_i + q_i}{\Omega _i \rho _i^2} \nabla _i W_{ij}(|\boldsymbol{r}_i - \boldsymbol{r}_j|, h_i)\right. \nonumber \\&\qquad \left. + \frac{P_j + q_j}{\Omega _j \rho _j^2} \nabla _i W_{ij}(|\boldsymbol{r}_i - \boldsymbol{r}_j|, h_j)\right] - \frac{G M_{*}}{r_i^2}(1 - \Gamma ), \end{aligned} $$(34)

D e i D t = P i Ω i ρ i 2 j m j ( v i v j ) · i W ij ( | r i r j | , h i ) + Λ , $$ \begin{aligned}&\frac{\mathrm{D} e_i}{\mathrm{D} t} = \frac{P_i}{\Omega _i \rho _i^2} \sum _j m_j (\boldsymbol{v}_i - \boldsymbol{v}_j) \cdot \nabla _i W_{ij}(|\boldsymbol{r}_i - \boldsymbol{r}_j|, h_i) + \Lambda , \end{aligned} $$(35)

where the indices i and j denote an SPH particle and its neighbours, respectively. Each SPH particle carries the following properties: the local density ρi, its (constant) mass mi, its position ri, smoothing length hi, velocity vi, pressure Pi, artificial viscosity qi, and specific internal energy ei. W(|ri − rj|, h) is the smoothing kernel, ∇iWij its gradient with respect to the location ri of particle i, and Ω accounts for the gradient of the smoothing length.

The last term in Eq. (35) represents the net cooling/heating rate of the SPH fluid. The last term in Eq. (34) corresponds to the net acceleration due to the contributions from the gravitational attraction by the stellar mass M* with G the gravitational constant and from an outward force, where Γ is the ratio of the outward and gravitational accelerations. For a dust-driven stellar wind, the outward force exerted on the gas is generated by a two-stage mechanism. Firstly, the dust in the wind absorbs the incident stellar radiation, which deposits its momentum and propels the dust outwards. Secondly, the dust collides with the surrounding gas, and drags it along. Assuming a full degree of mechanical coupling (also known as position coupling1), it can be shown that

Γ = ( κ d + κ g ) L 4 π c G M , $$ \begin{aligned} \Gamma = \frac{(\kappa _\mathrm{d} + \kappa _\mathrm{g} ) L_{*}}{4 \pi c G M_{*}}, \end{aligned} $$(36)

where κd and κg are the Planck mean opacity coefficients of the dust (see following section) and gas, respectively, L* is the luminosity of the star, and c the speed of light. Summing the gas and dust opacity contributions is consistent with our assumption of position coupling where the gas moves with the same speed as the dust. This remains a good approximation as long as the mass-loss rate is higher than ∼10−8M yr−1 (Jorissen & Knapp 1998).

The system of Eqs. (33)–(35) is supplemented with an equation of state (EOS) of the form P = f(ρ, e). We assume an ideal gas equation, i.e.,

P = ρ μ m u k T g , $$ \begin{aligned}&P = \frac{\rho }{\mu m_\mathrm{u} } k T_{\mathrm{g} }, \end{aligned} $$(37)

= ( γ 1 ) ρ e , $$ \begin{aligned} &\ \ = (\gamma - 1) \rho e, \end{aligned} $$(38)

where P is the total gas pressure, and mu the atomic mass unit. For the calculation of the mean molecular weight μ and the polytropic index γ of the gas it is sufficient to take into account only the most abundant species, i.e. H, He, and H2, giving

μ = ( 1 + 4 ϵ He ) n H k T g P H + P H 2 + ϵ He n H k T g , $$ \begin{aligned} \mu = \frac{(1 + 4 \epsilon _\mathrm{He} )\, n_{\langle \mathrm{H} \rangle } k T_{\mathrm{g} }}{P_\mathrm{H} + P_{\mathrm{H} _2} + \epsilon _\mathrm{He} n_{\langle \mathrm{H} \rangle } k T_{\mathrm{g} }}, \end{aligned} $$(39)

and following Grassi et al. (2014)

γ = 5 P H + 7 P H 2 + 5 ϵ He n H k T g 3 P H + 5 P H 2 + 3 ϵ He n H k T g . $$ \begin{aligned} \gamma = \frac{5 P_\mathrm{H} +7 P_{\mathrm{H} _2}+5\epsilon _\mathrm{He} n_{\langle \mathrm{H} \rangle } k T_{\mathrm{g} }}{3 P_\mathrm{H} +5 P_{\mathrm{H} _2}+3\epsilon _\mathrm{He} n_{\langle \mathrm{H} \rangle } k T_{\mathrm{g} }}. \end{aligned} $$(40)

In these expressions, n⟨H⟩ is the total number of hydrogen atoms per unit volume in the medium and ϵHe is the total number of helium atoms per hydrogen atom. In this formula it is assumed that H2 is the dominant molecule carrying hydrogen, and that all helium remains monoatomic. The contribution from the other species is also neglected but this remains a good approximation given that they do not exceed 1% of the total mass fraction. This also implies that the variations in μ and γ will only result from the formation of molecular hydrogen. Below 500 K, we use a simplified chemical network where all hydrogen is in the form of H2 and all carbon atoms locked in molecules (see Fig. A.1). This implies that dust formation is extremely inefficient at low temperatures yielding a very small value for the super-saturation ratio S.

We neglect the dust component in the EOS because the dust-to-gas ratio is at most 1 % in our simulations. Equation (37) is exact when the interaction energy of the particles is small compared to their kinetic energy, which is the case in circumstellar dust shells (Gail & Sedlmayr 2014, ch. 3.4).

In PHANTOM, the information carried by each SPH particles i are by default ri, vi, ei and hi (ρi is a simple function of hi). To these variables, we add J*, i, 𝒦0, i, 𝒦1, i, 𝒦2, i, 𝒦3, i, the mean molecular weight μi, and the polytropic index γi. Because the chemistry couples the mean molecular weight of the gas to its temperature, and the temperature is obtained by assuming a certain value for μi (i.e. μi = μ(T) and T = T(μi)), our models initially suffered from numerical oscillations around the recombination temperature of H2. These instabilities were eliminated after we introduced a sub-loop over μi, γi and T to iteratively obtain continuous self-consistent values. Every other quantity (molecular abundances, dust fraction, opacity, etc.) can be recovered from these variables.

3.2. Opacities

3.2.1. Gas opacity

For the mean gas opacity κg, we use the constant value proposed by Bowen (1988), κg = 2 × 10−4 cm2 g−1. To test the relevance of this approximation, we carried out detailed NLTE microphysics simulations over a wide range of densities and temperatures and confirm that the gas opacity is indeed of the order of 10−4 cm2 g−1 in our models. For these simulations, we used version 13.03 of the Cloudy code (Ferland et al. 2013), with a C abundance defined by C/O = 1.4, and solar abundances for the other elements. Figure 2 shows a contour plot of the Planck mean gas opacity as a function of H number density and temperature under a mean intensity of Jν = Bν(3000 K) representative of the radiation field of a central AGB star. The solid red line shows the trajectory of a wind particle wind in the nH − T plane and demonstrates that the constant value of κg proposed by Bowen is a good approximation for our simulations. Winters et al. (2000) analyzed in detail the impact of molecular opacities and showed that using a lower opacity leads to higher mass loss rates and terminal velocities.

thumbnail Fig. 2.

Planck mean values of the gas opacity κg (logarithmic scale), derived from a grid of models calculated with the Cloudy code using solar abundances with C/O = 1.4, and Jν = Bν(3000 K). The thick red line traces a wind particle and illustrates the range of densities and temperatures encountered in our models.

3.2.2. Dust opacity

In the framework of the Mie theory, the total opacity of a collection of spherical grains depends on the wavelength λ and grain size a and can be formulated as

κ d = 1 ρ 0 π a 2 Q ext ( a , λ ) n ( a ) d a , $$ \begin{aligned} \kappa _\mathrm{d} = \frac{1}{\rho } \int _0^\infty \pi a^2 Q_\mathrm{ext} (a,\lambda ) \, n(a)\, \mathrm{d} a , \end{aligned} $$(41)

where n(a) is the grain size distribution function (cm−4) and Qext the extinction efficiency. The determination of the opacity thus requires the knowledge of the distribution function. However, if the grain size is small compared to the wavelength of the stellar radiation field (2πa/λ ≪ 1), in the so-called small particle limit (SPL), the extinction efficiency is dominated by absorption and in this Rayleigh regime it can be expressed as Q e x t ( a , λ ) = Q e x t ( λ ) a $ Q_\mathrm{ext}(a,\lambda) = {Q}^{\prime}_\mathrm{ext}(\lambda)\,a $ allowing to separate in Eq. (41) the wavelength and grain size dependence leading to

κ d = π ρ Q ext ( λ ) 0 a 3 n ( a ) d a = π a 0 3 ρ Q ext ( λ ) K 3 , $$ \begin{aligned} \kappa _\mathrm{d} = \frac{\pi }{\rho } Q^\prime _\mathrm{ext} (\lambda )\, \int _0^\infty a^3 n(a)\, \mathrm{d} a = \frac{\pi a_0^3 }{\rho }\, Q^\prime _\mathrm{ext} (\lambda )\, \mathcal{K} _3, \end{aligned} $$(42)

where a0 is the radius of the monomer (Eq. 11), and 𝒦3 is the third moment as described in Sect. 2.3. For an AGB star with an effective temperature of ∼2500 − 3000 K, most of the flux is emitted around 1 μm. This implies that the SPL approximation will hold if a ≲ 0.1 μm. Beyond this limit the scattering contribution must be included and the full Mie theory used. From their non-grey radiative transfer models, Höfner & Dorfi (1992) concluded that the SPL approximation remained valid. Later Mattsson & Höfner (2011) mitigated this conclusion, mentioning that for models where the radiative-to-gravitational acceleration ratio is close to one, the SPL approximation breaks down but they also proposed a new approach that alleviates the limitation of the SPL approximation by using the Mie theory for one representative mean grain size. In future developments, we will implement this less restrictive and physically more consistent approach but if the grain size distribution shows a large variance, using one characteristic grain size may not give such a satisfactory result. Finally, it should be mentioned that recent models including non grey effects showed that scattering plays an essential role in launching winds in (O-rich) M-type AGB stars (Höfner 2008; Bladh et al. 2015), further emphasizing the need to go beyond the SPL approximation.

Various prescriptions are available for the Planck mean of these extinctions (Gail & Sedlmayr 1987b) and in our simulations we use the fit given by Draine (1981)

Q ext = 6.7 T d . $$ \begin{aligned} Q^\prime _\mathrm{ext} = 6.7\, T_{\mathrm{d} }. \end{aligned} $$(43)

Figure 3 shows the Planck mean absorption efficiencies (Qext) of amorphous carbon grains from various sources. The Qext’s were calculated for three representative grains sizes (a = 0.01, 0.1 and 1 μm) using the Mie theory for spherical grains (program BHMIE from Bohren & Huffman 1983 adapted by Draine2). The optical data were taken from Preibisch et al. (1993). These values are also compared to those obtained when approximating the wavelength dependence of the extinction efficiency in the SPL regime by a power law of the form Q e x t = Q 0 ( λ 0 / λ ) n $ Q^\prime_\mathrm{ext} =Q^\prime_0 (\lambda_0/\lambda)^n $ where the exponent varies between n = 1.19, 1.38 or 1.46 depending on the chosen optical data (Winters et al. 2000). As can be seen, the quality of the fits deteriorates with increasing grain size, further illustrating the limitation of the SPL approximation. However, provided the grain size does not exceed ∼0.1 μm, the agreement between the different prescriptions is good.

thumbnail Fig. 3.

Planck mean values of the dust extinction efficiency Qext for amorphous carbon grains with size a = 0.01, 0.1 and 1 μm (black solid lines, from Preibisch et al. 1993). For comparison, the approximations by Draine (1981; blue) and Winters et al. (2000) with n = 1.19 (magenta), 1.38 (green) and 1.46 (red) are shown.

4. Numerical implementation in PHANTOM

4.1. Inner wind boundary condition

To simulate mass loss from a central source, particles are released in the numerical domain from a specified distance to the star’s center called the injection radius Rinj. The properties of the injected particles (density, velocity, internal energy and chemical composition) are set by solving the equations below for a radially expanding wind.

In this framework, the Lagrangian form of the hydrodynamic wind equations reads

D ρ D t = ρ r 2 d ( r 2 v ) d r , $$ \begin{aligned}&\frac{\mathrm{D} \rho }{\mathrm{D} t} = - \frac{\rho }{r^2} \frac{\mathrm{d} (r^2 {v})}{\mathrm{d} r}, \end{aligned} $$(44)

D v D t = 1 ρ d P d r G M r 2 ( 1 Γ ) , $$ \begin{aligned}&\frac{\mathrm{D} {v}}{\mathrm{D} t} = - \frac{1}{\rho } \frac{\mathrm{d} P}{\mathrm{d} r} - \frac{G M_{*}}{r^2}(1 - \Gamma ), \end{aligned} $$(45)

D P D t = ρ c s 2 r 2 d ( r 2 v ) d r + ( γ 1 ) ρ Λ , $$ \begin{aligned}&\frac{\mathrm{D} P}{\mathrm{D} t} = - \frac{\rho c_\mathrm{s} ^2}{r^2} \frac{\mathrm{d} (r^2 {v})}{\mathrm{d} r} + (\gamma - 1) \rho \Lambda , \end{aligned} $$(46)

where cs is the adiabatic sound speed defined by c s 2 = P / ρ = γ P / ρ $ c_\mathrm{s}^2 = \partial P / \partial \rho = \gamma P / \rho $. Assuming a stationary wind, Eq. (44) reduces to the equation of mass conservation

ρ ( r ) = M ˙ 4 π r 2 v , $$ \begin{aligned} \rho (r) = \frac{\dot{M}}{4 \pi r^2 {v}}, \end{aligned} $$(47)

where is the constant mass-loss rate. In this case, DP/Dt = v dP/dr and Eq. (46) can be used to replace the pressure gradient in Eq. (45), giving

d v d r = 2 c s 2 / r G M ( 1 Γ ) / r 2 ( γ 1 ) Λ / v v ( 1 c s 2 / v 2 ) . $$ \begin{aligned} \frac{\mathrm{d}{v}}{\mathrm{d}r} = \frac{2 c_\mathrm{s} ^2 / r - G M_{*} (1 - \Gamma ) / r^2 - (\gamma - 1) \Lambda / {v}}{{v} (1 - c_\mathrm{s} ^2 / {v}^2)}. \end{aligned} $$(48)

The ideal gas law (Eq. (37)) can be used to convert Eq. (46) into an equation for the temperature

d T g d r = ( 1 γ ) T g ( 2 r + 1 v d v d r ) + ( γ 1 ) μ m u k Λ v . $$ \begin{aligned} \frac{\mathrm{d}T_{\mathrm{g} }}{\mathrm{d}r} =(1 - \gamma ) T_{\mathrm{g} }\left(\frac{2}{r} + \frac{1}{{v}}\frac{\mathrm{d} {v}}{\mathrm{d} r}\right) + \frac{(\gamma - 1)\, \mu m_\mathrm{u} }{k}\frac{\Lambda }{{v}} . \end{aligned} $$(49)

Finally the time evolution is given by

d r d t = v , $$ \begin{aligned} \frac{\mathrm{d}r}{\mathrm{d}t} = {v}, \end{aligned} $$(50)

We integrate Eqs. (48)–(50) using the fourth-order Runge-Kutta method, to obtain the physical properties of the ejected SPH particles as a function of both distance and time. For trans-sonic wind solutions, the numerator and denominator of Eq. (45) vanish at the sonic radius leading to an undefined fraction that causes numerical problems. To find the trans-sonic trajectory, the wind velocity at the inner boundary is varied until the numerator and the denominator are smaller than a threshold value of 10−2 in the vicinity of the sonic point. Then, L’Hôpital’s rule is applied to obtain the correct velocity gradient.

A free wind is obtained when the velocity gradient given by Eq. (48) is greater than or equal to zero beyond the sonic radius. However, there exist a set of solutions to Eqs. (48)–(50) for which the velocity gradient is smaller than zero, which are known as breeze-type solutions (Lamers & Cassinelli 1999, see region 3 in their Fig. 3.1). These appear when the numerator and the denominator of Eq. (48) have opposite signs. When v <  cs, that is when the wind is subsonic, the breeze is obtained as long as the outward acceleration factor Γ is higher than the threshold value

Γ max = 1 4 c s 2 v esc 2 + ( γ 1 ) Λ r 2 G M v , $$ \begin{aligned} \Gamma _\mathrm{max} = 1-\frac{4c_\mathrm{s} ^2}{{v}_\mathrm{esc} ^2}+(\gamma -1)\Lambda \frac{r^2}{GM_{*}{v}}, \end{aligned} $$(51)

where vesc = (2GM*(1 − Γ)/R*)1/2 is the escape velocity.

Our SPH simulations are unable to reproduce breeze-type solutions. Even if particles are launched from a hot corona with a velocity below the trans-sonic value, they ultimately converge on the trans-sonic solution. We speculate that the reason for this behavior is due to the fact that a breeze solution has a non vanishing pressure at infinity in contrast to the trans- and supersonic cases. Since our SPH simulations assume a free boundary condition and therefore zero pressure at infinity, the particles tend to converge toward that solution. A fixed outer boundary providing a non-zero outer pressure could be implemented using a method similar to the “handled-shells” method for particle injection described in the following section. However, a subsonic breeze was not detected in the solar system by the Mariner 2 space craft (e.g., Neugebauer & Snyder 1962) and this solution is also formally unstable (Velli 2001). Given the limited interest of such types of wind, they are excluded from our study.

4.2. Injection of SPH particles

At the inner boundary, the injected SPH particles are homogeneously distributed on a shell of radius Rinj. Taking a constant number p of particles per shell, the time interval δt between two subsequent shells is given by the relation

δ t = p m / M ˙ , $$ \begin{aligned} \delta t = p m / \dot{M}, \end{aligned} $$(52)

where m is the constant mass of a SPH particle and the wind mass loss rate. In this framework, the k’th shell is injected at time t ~ k = k δ t $ \tilde{t}_k = k \delta t $. Hence, the number of shells injected at every time step is set by the fraction of δt over the SPH timestep Δt.

Using the age of the simulation as a way to index the shells, we can track how consecutive shells should be launched. If tn + 1 is the time to be reached and Δtn the value of the current hydro (SPH) time step, then the shells with index imin to imax have to be ejected, where

i min = t n + 1 Δ t n δ t + 1 , $$ \begin{aligned}&i_\mathrm{min} = \left\lfloor \frac{t^{n+1} - \Delta t^n}{\delta t}\right\rfloor + 1, \end{aligned} $$(53)

i max = t n + 1 δ t , $$ \begin{aligned}&i_\mathrm{max} = \left\lfloor \frac{t^{n+1}}{\delta t}\right\rfloor , \end{aligned} $$(54)

(“⌊ ⌋” is the floor function). When Δt <  δt, imin may become greater than imax, in which case no shell is released.

However, in general the hydrodynamical time tn + 1 does not coincide with the precise time i at which injection should occur. Consequently, if tn+1 > i the properties of the released particles have to be synchronized and estimated at the “current” time, that is at a time t e = t n + 1 t ~ i $ t_e=t^{n+1}-\tilde{t}_i $ after their ejection. The velocity is given by the solution v(te) of Eq. (48), and the density by Eq. (47), converted to ρ(te) using the solution of Eq. (50). The internal energy (see Eqs. (37) and (38)) is calculated using the temperature from Eq. (49). The radius r(te) at which the new shell is freed is given by Eq. (50).

In our SPH simulations, the central star is represented by a sink particle (Bate et al. 1995; Price et al. 2018), which can not exert pressure on its surroundings, and thus produces a central vacuum. This implies that any particle in the vicinity of the sink will feel a positive pressure gradient pulling material towards the sink. Consequently, simulating a wind by adding only one shell of SPH particles at a time will not result in an effective particle injection. To compensate for this effect, a number of radially fixed, consecutive shells are placed around the sink. At each time step the density, velocity, and internal energy of the particles on these boundary shells are set to the analytical radial wind solution derived in Sect. 4.1. This procedure ensures a correct behaviour of the physical variables at the inner boundary and in particular of the pressure gradient. The farthest handled shell dictates the radius beyond which the properties of the SPH particles are no longer imposed by the analytical solution and can evolve freely, as dictated by the SPH solver. This is the effective injection radius Rinj. The number of fixed boundary shells can be specified in the code with an integer parameter (iboundary_shells). Our tests show show that 5 shells are able to accurately provide the expected pressure gradient.

In Fig. 4, we show a schematic diagram of what happens when the hydrodynamical time step Δt is greater than the theoretical time between the release of consecutive shells δt. In this example three handled shells define the inner boundary, and shell with index number 12 is released.

thumbnail Fig. 4.

Schematic representation of the method for injecting SPH particles. The horizontal axes represent the time since ejection, or interchangeably the distance from the boundary. The numbered dots represent particle shells, numbered by their index. Older shells have a lower index. Top: before the hydro timestep the shell with index 12 lies below the injection boundary, and is therefore part of the handled shells, and has imposed physical properties. Bottom: after the hydro timestep, the shell with index 12 has moved beyond the injection boundary, and no longer has imposed physical properties.

4.3. Distributing particles on geodesic spheres

To distribute a finite number of particles as uniformly as possible on the surface of a sphere, we employ the icosahedron interpolation method (Wang & Lee 2011). Using this method the triangular faces of an icosahedron are tiled with smaller triangles. The vertices of all these triangles are then projected onto the platonic shell of the original icosahedron, and the SPH particles are placed on the projected vertex coordinates.

Each triangular face of the icosahedron is subdivided in such a way that every edge contains 2q vertices, including the corner points. Thus, there are 20(2q − 1)2 triangular faces and p = 10(2q − 1)2 + 2 vertices in each shell3. To ensure that the projected vertices are as equidistant as possible we use the mapping procedure described by Tegmark (1996), which ensures that every projected facet has the same area. This is necessary to reduce initial inhomogeneities on the sphere and numerical artifacts in the simulations.

As the geodesic spheres are not perfectly isotropic, it is necessary to rotate each sphere with respect to the previous one to prevent anisotropies from building up. The optimal rotation angles ϑi (i = 1, 2, 3) that will maximize the minimal distance between vertices of two subsequent geodesic spheres is given for various values of q in Table 2. We determined those angles empirically by testing thousands of random combinations.

o ̂ = ( c 2 c 3 c 2 s 2 s 2 s 1 s 2 c 3 + c 1 s 3 s 1 s 2 s 3 + c 1 c 3 s 1 c 2 c 1 s 2 c 3 + s 1 s 3 c 1 s 2 s 3 + s 1 c 3 c 1 c 2 ) o ̂ , $$ \begin{aligned} \hat{\boldsymbol{o}}^{\prime } = \begin{pmatrix} c_2 c_3&-c_2 s_2&-s_2 \\ -s_1 s_2 c_3 + c_1 s_3&s_1 s_2 s_3 + c_1 c_3&-s_1 c_2 \\ c_1 s_2 c_3 + s_1 s_3&-c_1 s_2 s_3 + s_1 c_3&c_1 c_2 \end{pmatrix} \hat{\boldsymbol{o}}, \end{aligned} $$(55)

Table 2.

Optimized rotation angles (in rad) for the distribution of SPH particles at the inner boundary as a function of the number of vertices or particles (Np = 10(2q − 1)2 + 2).

where ci = cosϑi and si = sinϑi. o ̂ $ \hat{\boldsymbol{o}} $ is defined as a unit vector pointing towards an SPH particle on a sphere, and the subsequent sphere will be identical to the previous but rotated in such a way that the vector changes into o ̂ $ \hat{\boldsymbol{o}}\prime $.

4.4. Wind resolution

There is only a finite number of ways in which particles can be arranged onto a shell using the icosahedron interpolation method described in the previous section. This inherent quantization implies that the tangential resolution of the inner boundary condition can only take on discrete configurations. It can be shown that the distance d between neighbouring particles on such a shell of unit radius is equal to

d = 2 ( 2 q 1 ) ( φ 5 ) 1 / 2 , $$ \begin{aligned} d_\mathrm \perp = \frac{2}{(2 q - 1)} \left(\varphi \sqrt{5}\right)^{-1/2}, \end{aligned} $$(56)

where φ = ( 1 + 5 ) / 2 $ \varphi = (1 + \sqrt{5}) / 2 $ is the golden ratio, and q can be any natural number, which represents the tangential resolution set-up, and appears in the PHANTOM parameter card as the iwind_resolution variable.

Conversely, there is no geometrical constraint on the resolution long the radial axis. Given the inherent length scale associated with Eq. (56), it is natural to express the resolution along the radial axis as a function of this distance. In PHANTOM, the radial distance between two consecutive shells is defined as

δ r = w ss d , $$ \begin{aligned} \delta r = { w}_\mathrm{ss} d_\mathrm \perp , \end{aligned} $$(57)

where wss can take on any value larger than zero, and can also be set in the PHANTOM parameter card though the variable wind_shell_spacing.

Ideally, the resolution parameters should be chosen such that two primary requirements are satisfied. Firstly, the smoothing length of the injected SPH particles should be larger than the distance between them, both radially and tangentially. This is necessary to ensure that the particle kernels overlap and effectively feel each other along these directions. Secondly, the representative length-scale of the system must be resolved. We address how to approach both these requirements in wind systems.

The smoothing length h is given by hfactn−1/3 (Price et al. 2018), where n is the particle number density, and hfact is a proportionality factor that can vary with the adopted kernel and which is set by default to 1.2. As a consequence, h can be re-written as

h = h fact d r 1 3 d 2 3 , $$ \begin{aligned} h = h_\mathrm{fact} \, d_r^{\frac{1}{3}} d_{\perp }^{\frac{2}{3}}, \end{aligned} $$(58)

where d is the tangential particle distance given by Eq. (56), and dr ≈ δr = wssd is the radial distance between particle on adjacent shells. Hence, the first resolution requirement is satisfied when both the h ≥ dr and h ≥ d conditions are simultaneously fulfilled. Because dr depends on d, this translates to a resolution condition on wss given by

h fact 3 w ss h fact 3 2 . $$ \begin{aligned} h_\mathrm{fact} ^{-3} \le { w}_{\rm ss} \le h_\mathrm{fact} ^\frac{3}{2}. \end{aligned} $$(59)

The second resolution requirement can be satisfied by choosing the parameter q in Eq. (56) such that the radial distance between consecutive shells is smaller than the representative length-scale of the system. For a Parker wind (Lamers & Cassinelli 1999), this length-scale can be identified to the sonic radius. For a wind-Roche-lobe-overflow binary (Mohamed & Podsiadlowski 2012), the representative length-scale would be set by the position of the first Lagrange point and for detached binaries by the orbital separation. Since a typical number of handled shells (see Sect. 4.2) of five is required to properly generate the inner boundary condition on the pressure, one should aim towards setting the parameter q such that ∼ ten shells can fit between the inner boundary and the relevant critical radius.

Given the above constraints on the spatial resolutions, the mass of the SPH particles (which is typically considered as the standard measure for the resolution in SPH), is fully constrained by the choice of the mass-loss rate and initial wind velocity which are input parameters of PHANTOM. This is because the distance δr between two adjacent shells can also be approximately given by

δ r v inj N p m / M ˙ , $$ \begin{aligned} \delta r \simeq {v}_\mathrm{inj} \,N_p m / \dot{M}, \end{aligned} $$(60)

where vinj is the wind velocity at the injection radius Rinj (also available in the PHANTOM parameter card as wind_velocity), and Np is the number of particles4 of mass mp. Given that the inter-shell distance is also given by Eq. (57), this implies that the mass of the particles is determined by the set-up of the boundary conditions and

m p = w ss d v inj M ˙ N p . $$ \begin{aligned} m_p = { w}_\mathrm{ss} \frac{d_\mathrm \perp }{{v}_\mathrm{inj} }\frac{\dot{M}}{N_p}. \end{aligned} $$(61)

The user can alternatively impose the particle mass and the code will find the closest value of the parameter q that fulfills the above resolution criteria.

5. Benchmarking

To validate the implementation of the above-mentioned methods, we calculate a sample of benchmarking simulations and compare the analytic solution to the 3D SPH wind profiles. In all simulations shown here, the net cooling rate Λ in Eqs. (35), (46) and (48) is set to zero. Furthermore, the quintic kernel is used; its larger compact support radius contributes to reduce the dispersion of the SPH solution.

5.1. Stationary dust-free winds

In this section, we discuss the performance of the implemented wind physics in different idealised thermodynamical contexts. In particular, we investigate the ability of the code to predict transonic and supersonic wind dynamics. Consequently, we focus the discussion on the velocity profile of the wind. For all studied cases, the physical parameters are described in Table 3. Furthermore, the resolution parameter q is set to 15 (see Eq. (56)), and the spacing of the shells satisfies the condition in Eq. (59), unless stated otherwise. In practice, we use wss = 1.0. Finally, to position the sonic radius beyond the radius of the handled shells an initial sound speed cs, iso = 13.12 km s−1 was assumed (corresponding to Tiso = 1.25 × 104 K). In all cases presented here and in the following sections, conditions resulting in a breeze-type solution (see Sect. 4.1) were deliberately avoided because the SPH solver was unable to recover it properly.

Table 3.

1D model parameters. Tiso corresponds to the temperature of the isothermal wind.

Figure 5 shows the velocity profiles of an isothermal SPH wind solution, in the transonic and the supersonic regimes. Shown in black is the analytic solution as derived in Section 4.1. The SPH solution shows good agreement with the analytic Parker wind prediction, exhibiting Gaussian 1σ dispersion values of 0.22 km s−1 and 0.35 km s−1 for the supersonic and transsonic solutions, respectively. However, in the subsonic regime of the transonic solution small oscillations can be seen in the SPH solution when the wind is released by the last boundary shell. We henceforth refer to these oscillations as sonic-point oscillations, which are the consequence of the high sensitivity of the transonic solution to the initial conditions. When the SPH particles are released from the last handled shell, the SPH particles rearrange and disperse, which puts the particles on a distribution of interacting quasi-transonic trajectories. However, these oscillations do not affect the overall behaviour of the SPH gas, since beyond the sonic point they die down and the SPH solution nicely follows the analytic solution.

thumbnail Fig. 5.

Isothermal velocity profiles of a transonic (bottom) and supersonic (top) SPH wind solution, in red. The analytic solution is given by the black curves for a wind temperature Tiso = 1.25 × 104 K.

Similar benchmarks were performed for an adiabatic wind, with a polytropic index γ = 1.4 and are illustrated in Fig. 6. As for the isothermal case, both the trans- and super-sonic simulations tightly follow the analytic Parker wind solution, exhibiting Gaussian 1σ dispersion of 0.16 km s−1 and 0.29 km s−1 for the supersonic and transsonic cases, respectively. Very weak sonic-point oscillations can be seen for the sub-sonic portion of the trans-sonic simulation.

thumbnail Fig. 6.

Adiabatic velocity profiles of a transonic (bottom) and supersonic (top) SPH wind solution, in red. The analytic solution is given by the black curves.

In Fig. 7, we show the effect of not properly adjusting the resolution to resolve the critical length-scale of the simulation. We demonstrate the effect by performing a resolution study on the adiabatic transsonic wind case. From the bottom to the top panel, the resolution parameter q is increased from 3 to 10 to 20, which corresponds to a factual decrease of the SPH particle mass with a factor ∼10 at each step. The number of handled shells is 5, which for the simulation with q = 3 (bottom panel) nearly completely fill up the subsonic region. Consequently, the transonic passage is not properly sampled, and the SPH solution sways around the analytic curve. For the simulations at higher resolution, the transonic passage is nicely resolved, and the SPH solution tightly follows the analytic curve.

thumbnail Fig. 7.

Effect of the resolution setup on the predicted dynamics of an transonic adiabatic wind, for a fixed value of wss. The value of q is chosen such that SPH particle mass is decreased by a factor ten from the top to the middle panel, and again from the middle to the bottom panel.

5.2. Effect of Bowen dust opacity prescription

In this section, we discuss the effect on the SPH wind dynamics of the activation of Γ in Eq. (34). As in the previous section, we focus solely on the velocity profile of the wind. For all studied cases, the resolution parameter q is set to 10, the number of handled shells is five, and the spacing of the shells satisfies the condition in Eq. (59), unless otherwise stated.

In his seminal paper, Bowen (1988) proposed a very simple parametrisation of the dust opacity as a function of temperature. In this formulation, the dust opacity writes

κ d ( T eq ) = κ max 1 + exp [ ( T eq T cond ) / δ ] , $$ \begin{aligned} \kappa _\mathrm{d} (T_\mathrm{eq} ) = \frac{\kappa _\mathrm{max} }{1+\exp [(T_\mathrm{eq} -T_\mathrm{cond} )/\delta ]}, \end{aligned} $$(62)

where κmax is the maximum value the dust opacity can reach, δ is a measure for the temperature range over which condensation occurs, and Tcond the dust condensation temperature which depends on the star’s composition (C/O ratio). Teq is the dust equilibrium temperature which value depends on how the gas and dust are thermally coupled. A correct determination of Teq would require radiative transfer calculations which is outside the scope of this paper (see Sect. 6) so for the moment we assume that Teq is set to the gas temperature.

Though the Bowen parameterization is a severe approximation of the opacity-forming processes, we use this prescription to showcase in Fig. 8 the activation of dust formation on the SPH dynamics. In this setup, a polytropic wind is launched with an initial speed of 5 km s−1 from a 1.5 M cool star with a surface temperature of 3000 K and a luminosity of 2 × 104L and with μ = 2.38, corresponding to hydrogen in molecular form. The injection radius is set to 1.1 au. Under these initial conditions, a free wind without pulsations and cooling, the particles would ballistically fall back onto the stellar surface. These critical ingredients will be the topics of follow-up papers, as discussed in Sect. 6. To circumvent this problem we introduced an additional acceleration parameter α in Eqs. (34) and (45) to the Eddington factor Γ. This parameter α, takes on a value between zero and one, and serves solely to counteract the gravitational potential. Although this parameter is used in an ad-hoc manner in our context, it could in principle be related to the energy injected by the pulsations and the propagation of a pressure wave (e.g., Struck et al. 2004; Mattsson 2016). This technique is commonly used in hydrodynamical wind simulations that do not calculate dust acceleration (e.g., Mastrodemos & Morris 1998; Kim & Taam 2012; Liu et al. 2017). Consequently, the gravity term in Eqs. (34) and (45) now reads

thumbnail Fig. 8.

Acceleration profile of a free, steady, polytropic wind with the Bowen dust opacity prescription activated around ∼ 2 au. The temperature is shown in blue.

G M r 2 ( 1 α Γ ) , $$ \begin{aligned} - \frac{G M_{*}}{r^2} (1 - \alpha - \Gamma ), \end{aligned} $$(63)

where the meaning of all other parameters has remained unchanged.

The profile in Fig. 8 was generated by setting α equal to unity and an initial velocity at the injection radius of vinj = 5 km s−1. We note that setting α = 1 will effectively overestimates the true acceleration of the wind. This is of no consequence for this exercise, since the aim is to illustrate the mechanism of dust acceleration in action. We set the dust condensation temperature to 1500 K, and the width of the condensation temperature range δ to 60 K and the polytropic index of the gas to γ = 1.4. Finally a value of κmax = 1.2 cm2 g−1 is chosen, yielding Γ(Teq = Tcond)≃0.6 (using Eqs. (36) and (62)).

The SPH gas (red in Fig. 8) undergoes a double acceleration. The first stage, which occurs between 1 and ∼2 au, is a consequence of the inner wind pressure gradient caused by the stellar temperature boundary condition. Thanks to the parameter α set to unity in Eq. (63), the gas can freely escape and accelerate outwards. At a radius of ∼2 au, the temperature of the gas drops below the condensation temperature which activates the dust opacity and raises Γ in Eq. (36). This results in a second acceleration that brings the velocity of the gas around 37 km s−1. The dynamics of the SPH particles nicely follows the analytic curve. Especially in the first acceleration stage, the analytic curve (black) is centered on the SPH solution. When dust starts forming and triggers the secondary acceleration, the velocity of the SPH is slightly underpredicted with respect to the analytical solution. It appears that when the radiative acceleration is reduced, for example when the value of kmax in Eq. (62) is lowered, the agreement with the theoretical line is even better. This effect is likely related to the accuracy of the time integration scheme which is 2nd order accurate in SPH compared to the 4th order Runge-Kutta scheme used for the 1D solution. But the overall agreement is within a few percent. As a final note, we remark that the low surface temperature of the star implies that the SPH particles are launched fully supersonically. Consequently, the instabilities associated with passing the sonic point do not arise in this simulation (see Figs. 5 and 6).

5.3. Steady wind with dust nucleation

In Fig. 9, we show the velocity profile of a free (α = 1), steady wind that includes the full computation of the dust-nucleation described in Sect. 2. The numerical setup is nearly identical to the one used with the Bowen dust-opacity prescription described above, except that now the composition of the gas, and consequently its mean molecular weight μ and polytropic index γ, are set by the solution of the chemical network described in Sect. 2.2. The chemistry is set to a carbon to oxygen ratio of two. And contrarily to the previous case studies, the velocity and temperature profiles now depend on the value of the mass loss rate because the chemical and nucleation processes are explicit functions of the density. Indeed, the independence of the stationary wind velocity profile on comes from the fact that in the absence of cooling and chemistry, the density ρ acts only as a scaling factor in Eqs. (48) and (49). Integration of these equations yield rho-invariant functions v(r) and T(r), with the density profile then deduced from Eq. (47) that directly scales with the mass loss rate. For the simulation showcased in Fig. 9, a mass-loss rate of 10−5M yr−1, and an initial velocity of vinj = 6 km s−1 was assumed. This initial velocity differs slightly from the Bowen-type dust-opacity-setup in the previous section to avoid the breeze-type solution.

thumbnail Fig. 9.

Acceleration profile of a free wind for which the dust opacity is calculated using the nucleation formalism.

The effect of including dust nucleation on the wind dynamics is significant, and can be readily assessed by comparing the SPH solution (red) and the dust-free profile (black dashed) in Fig. 9. While the SPH solution reaches a terminal velocity of nearly 30 km s−1, its dust-free counterpart (accelerated only via the pressure-gradient at the stellar boundary and α = 1) barely reaches a speed of 12 km s−1. The dusty wind also exhibits a double-staged acceleration. The inner-most velocity profile of the dust-free curve matches the SPH solution, indicating that the first acceleration stage is due to the inner boundary conditions similar to what was found in the previous section (see Fig. 8). However, around 3 R* a much steeper acceleration of the SPH gas can be seen, which is not reproduced by the dust-free curve, and is therefore solely associated with dust formation.

The formation of dust can be quantified either by the rate of condensation from the gas-phase, of by their effect on the dust opacity. Both these quantities are shown in the left panel of Fig. 10. The sudden increase in the nucleation rate 𝒦3 (blue) from about zero to ∼106 cm−3 around ∼3R* coincides with a dramatic increase in dust opacity (green), which explains the sudden acceleration of the SPH particles. At ∼4R*, the nucleation rate reaches a plateau, and the corresponding dust opacity decreases as a consequence of decreasing density.

thumbnail Fig. 10.

Radial profiles of chemical and dust-related quantities. The temperature of the SPH gas, with which the dust is assumed to be fully thermalised, is shown in red. Left: the logarithm of the third moment 𝒦3 (blue), which represents the number density of condensing monomers from the gas phase into the dust, shown together with the dust opacity profile (green). Right: the mean molecular weight μ (blue), polytropic index γ (green), and the base 10 logarithm of the supersaturation ratio S (magenta) of the gas in the wind. A supersaturation ratio above unity implies favourable conditions for dust formation.

The supersaturation ratio S, shown in purple in the right panel of Fig. 10, can also be used as a proxy for dust formation. For values above unity, the partial pressure of gas-phase carbon exceeds that of solid-state carbon, indicating favourable condensation conditions. Indeed, the sudden increase in relative condensation rate 𝒦3 seen in the left panel of Fig. 10 occurs as S becomes much larger than unity. At distances slightly larger than 4 R* the temperature reaches 500 K, below which the simplified chemistry is activated (all H and C atoms locked in molecules) producing the artificial discontinuity in S (Sect. 3.1).

The radial decrease in gas temperature causes the hydrogen in the wind to transition from being fully monatomic to fully fully diatomic (see also Fig. A.1). This is reflected by the rapid change in mean molecular weight μ and polytropic index γ of the gas, seen in the right panel of Fig. 10. These quantities directly affect the temperature of the gas. Given that these simulations do not include cooling, except for adiabatic compression and expansion, this results in a small bump in the radial temperature profile.

6. Discussion and future prospects

The simulations performed in Sect. 5 show that our implementation of the wind-injection boundary, as well as the calculations involving dust opacity can be used to model cool and dusty stellar winds. However, to level-up the modelling of dust-driven winds to match the state-of-the-art knowledge on the physical mechanisms governing these outflows, many additional physical and chemical ingredients have yet to be accounted for (a comprehensive summary of theoretical challenges can be found in Haworth et al. 2016).

The stellar surface pulsations, driven by convective motions are an essential ingredient for launching a dusty wind and have been circumvented in our approach by enforcing a “free” wind by setting α = 1. This is a common assumption in the modelling of 3D winds (e.g., Jahanara et al. 2005; Kim & Taam 2012; Liu et al. 2017) but when dust formation is involved a more correct treatment of the wind launching mechanism is needed (Freytag et al. 2017). We will adopt a dynamical model as originally done by Bowen (1988), where the pulsations are simulated by a radially oscillating inner boundary, which can be considered to act as a piston located at the base of the photosphere. A recent implementation of this approach in the SPH context is described in Aydi & Mohamed (2022). Considering the modelling of large-scale convection as the source of pulsation as done in Freytag et al. (2017) would require the implementation of an energy transport equation in PHANTOM and is beyond the scope of our ambitions for this project.

Another crucial issue is related to the chemistry and cooling processes (Boulangier et al. 2018). The use of a polytropic index does not provide an accurate estimate of the thermodynamics of the gas because several processes such as radiative cooling by line transitions, chemical, and continuum cooling, will significantly modify the internal energy. In turn, this will impact the dust formation process and consequently the wind acceleration. In the context of wind-companion interactions, the absence of cooling has been shown to prevent the formation of an accretion disk around an AGB companion (Theuns & Jorissen 1993), with a dramatic impact on the derived accretion and angular momentum transfer rates and potentially on the secular evolution of the binary system. Furthermore, it will impact the flow morphology, affecting the interpretation of key features in the observations.

Related to this issue, is the correct treatment of the radiative acceleration which requires the resolution of the radiative transfer equations. If shocks form in the outflow, the opacity can locally increase dramatically and, as we showed, provide an extra acceleration which will impact the wind morphology. The treatment of radiative transfer is a very challenging and computationally intensive endeavour but is not beyond reach. Different approaches has been explored to solve the radiative transfer equation and are compared in the context of cosmological simulations, in Iliev et al. (2006, 2009). They include ray tracing algorithms where each source casts a number of rays and the equation of radiative transfer (RT) is integrated along the ray (e.g., Kessel-Deynet & Burkert 2000; Susa 2006; Altay & Theuns 2013; Chen et al. 2020). If scattering is neglected, this method reduces to calculating the optical depth along the ray. Another option is to solve the angular moments of the radiative transfer equation using specific closure relations and/or flux limiters (e.g., Gnedin & Abel 2001; Whitehouse et al. 2005; Aubert & Teyssier 2008; Petkova & Springel 2009; Skinner & Ostriker 2013; Chan et al. 2021). Other methods involving transport on unstructured meshes such as Delaunay tessellations (e.g., Rijkhorst et al. 2006; Pawlik & Schaye 2008; Petkova & Springel 2011; Petkova et al. 2021) have also been proposed. A final option is to delegate the resolution of the RT problem to an external code. This has been done in PHANTOM in dusty proto-planetary disk simulations where the output of the SPH simulation were regularly processed by the RT code MCFOST (e.g., Pinte et al. 2019) providing updated values for the gas and dust temperatures of each particles. Considering that we are mostly interested in “simple” systems consisting of one or two sources, the ray tracing approach seems to be a promising avenue. We will explore this possibility in a the future following an approach similar to Hasegawa & Umemura (2010).

A great simplification of the presented model is the absence of dust-gas drift. This process is important because conceptually the absence of drift is inconsistent with the dust-driven wind mechanism, which says that the gas is dragged along by the faster, radiatively accelerated dust particles (see e.g., discussion in Sect. 4.1 of Mattsson & Sandin 2021). In addition, as shown by Krüger et al. (1994), Sandin & Höfner (2004), Sandin (2008), Sandin & Mattsson (2020), the effect of drift can potentially lead to an increase in the dust production and, in some cases, to a decrease of the mass-loss rates and wind velocities. Fortunately, PHANTOM allows the modelling of dust-gas mixtures in either the standard two-fluid approach where each species is regarded as a separate fluid (Laibe & Price 2012) or in the so-called one-fluid approximation (Laibe & Price 2014a,b; Price & Laibe 2015; Hutchison et al. 2018) where the evolution of the gas and dust velocities require the tracking of only one additional quantity, the mass-fraction of the dust-grain mixture. This “diffusion approximation for dust” holds as long as the stopping time-step is small compared to the hydrodynamical timestep which should be the case for AGB winds in a steady regime outside shocked regions (see Sect. 6.8.3 of Gail & Sedlmayr 2014). Using this particular prescription overcomes some of the limitation associated with the two-fluid approach, and allows the use of much larger timesteps, which speeds up the calculations. We are currently working on including this formalism in our dusty wind simulations.

The road to reach a more consistent picture is still long and will proceed step-wise, but even with some missing physics, our simulations (Maes et al. 2021; Malfait et al. 2021) can already help understand wind morphologies around AGB stars.

7. Summary

In this paper, we present the implementation of a dusty wind into the smoothed particle hydrodynamics code PHANTOM. The SPH particles are initially distributed on a geodesic sphere and their properties determined by solving the 1D stationary wind equations. To counteract the buildup of a positive pressure gradient at the inner boundary due to the presence of a sink particle at the center, we force a number of fixed boundary shells, whose properties are imposed by the 1D wind equations.

Using this set-up, we have calculated a number of different benchmarking simulations that test the ability of the code to recover the stationary Parker-wind solutions. We show that the SPH solution tightly follows the expected supersonic and transonic wind trajectories. For the transonic case, the SPH solution shows weak oscillations in the subsonic portion of the trajectory. This is a consequence of the fact that the transonic trajectory is very sensitive to the precise initial conditions of the launched particles. The mutual interaction between the SPH particles makes them deviate slightly from this very narrow path, leading to oscillations in their velocity. Nevertheless, when the sonic point is eventually crossed oscillations die down leading to a smooth SPH solution.

The dust component of the wind is implemented according to two different schemes. The first is via an analytical expression that describes how the dust opacity increases once a dust condensation temperature threshold is crossed, as described by Bowen (1988). In the second scheme dust nucleation is calculated from the gas composition of the outflow. This composition is obtained by solving a compact carbon-rich chemical network assuming local-thermodynamical-equilibrium. This leads to the formation of carbonaceous monomers, that nucleate into dust particles. The dust formation formalism follows the moment equations prescription described by Gail & Sedlmayr (2014). This provides a quick way to estimate the rate at which the monomers condense into the dust grains, which can be directly related to the dust opacity.

Assuming a free (α = 1) wind, in which the gas and the dust are thermally and mechanically coupled (Tg = Td and with no drift velocity), both schemes lead to the establishment of an significant additional opacity term that dramatically impacts the wind dynamics by providing an additional outward acceleration. For the nucleation models, we show how the partial pressures of the carbon-compounds affect the supersaturation ratio, how this directly impacts the rate at which the carbonaceous monomers condense into the dust, and how this further coincides with the vast opacity increase that accelerates the wind.

As such, we illustrate that our implementation of the SPH particle injection set-up, and the nucleation and dust growth prescriptions in the PHANTOM code can be used to successfully model dusty stellar outflows.


1

In the position coupling regime, the dust and gas particles have the same position at all time and the gas-dust mixture is then one fluid. This also implies that gas and dust move at the same speed though the flow.

3

A comprehensive description of this tiling is provided (in french) at https://fr.wikipedia.org/wiki/Géode_(géométrie)

4

The number of particles on the isocahedron sphere is given by Np = 10(2q − 1)2 + 2 where the parameter q corresponds to the variable iwind_resolution in the PHANTOM input file.

Acknowledgments

The authors are in debt to H.-P. Gail and C. Dreyer who provided us a long time ago the source of their code that were extremely useful to check and validate our dust routines. We also thank the anonymous referee for very constructive and encouraging comments. L.S. is senior FRS-F.N.R.S. research associate. W.H. acknowledges support from a FRS-F.N.R.S grant. D.J.P. is grateful for Australian Research Council funding via DP180104235 and FT130100034.

References

  1. Altay, G., & Theuns, T. 2013, MNRAS, 434, 748 [NASA ADS] [CrossRef] [Google Scholar]
  2. Andersen, A. C., Höfner, S., & Gautschy-Loidl, R. 2003, A&A, 400, 981 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  3. Aubert, D., & Teyssier, R. 2008, MNRAS, 387, 295 [NASA ADS] [CrossRef] [Google Scholar]
  4. Aydi, E., & Mohamed, S. 2022, MNRAS, 513, 4405 [NASA ADS] [CrossRef] [Google Scholar]
  5. Bate, M. R., Bonnell, I. A., & Price, N. M. 1995, MNRAS, 277, 362 [Google Scholar]
  6. Bermúdez-Bustamante, L. C., García-Segura, G., Steffen, W., & Sabin, L. 2020, MNRAS, 493, 2606 [CrossRef] [Google Scholar]
  7. Bladh, S., Höfner, S., Aringer, B., & Eriksson, K. 2015, A&A, 575, A105 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  8. Bladh, S., Liljegren, S., Höfner, S., Aringer, B., & Marigo, P. 2019, A&A, 626, A100 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  9. Bohren, C. F., & Huffman, D. R. 1983, Absorption and Scattering of Light by Small Particles (New York: Wiley) [Google Scholar]
  10. Boulangier, J., Clementel, N., van Marle, A. J., Decin, L., & de Koter, A. 2018, MNRAS, 482, 5052 [Google Scholar]
  11. Bowen, G. H. 1988, ApJ, 329, 299 [NASA ADS] [CrossRef] [Google Scholar]
  12. Chan, T. K., Theuns, T., Bower, R., & Frenk, C. 2021, MNRAS, 505, 5784 [NASA ADS] [CrossRef] [Google Scholar]
  13. Chase, M. W. 1986, in NIST-JANAF Thermochemical Tables, 4th edn. (New York: American Chemical Society and American Institute of Physics for the National Bureau of Standards and Technology), Monograph, 9 [Google Scholar]
  14. Chen, Z., Blackman, E. G., Nordhaus, J., Frank, A., & Carroll-Nellenback, J. 2018, MNRAS, 473, 747 [NASA ADS] [CrossRef] [Google Scholar]
  15. Chen, Z., Ivanova, N., & Carroll-Nellenback, J. 2020, ApJ, 892, 110 [NASA ADS] [CrossRef] [Google Scholar]
  16. de Val-Borro, M., Karovska, M., Sasselov, D. D., & Stone, J. M. 2017, MNRAS, 468, 3408 [NASA ADS] [CrossRef] [Google Scholar]
  17. Dorfi, E. A., & Höfner, S. 1991, A&A, 248, 105 [NASA ADS] [Google Scholar]
  18. Draine, B. T. 1981, ApJ, 245, 880 [NASA ADS] [CrossRef] [Google Scholar]
  19. Draine, B. T., & Salpeter, E. E. 1977, J. Chem. Phys., 67, 2230 [NASA ADS] [CrossRef] [Google Scholar]
  20. El Mellah, I., Bolte, J., Decin, L., Homan, W., & Keppens, R. 2020, A&A, 637, A91 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  21. Ferland, G. J., Porter, R. L., van Hoof, P. A. M., et al. 2013, Rev. Mex. Astron. Astrofis., 49, 137 [Google Scholar]
  22. Ferrarotti, A. S., & Gail, H.-P. 2002, A&A, 382, 256 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  23. Fleischer, A. J., Gauger, A., & Sedlmayr, E. 1992, A&A, 266, 321 [NASA ADS] [Google Scholar]
  24. Forsythe, W. E. 2003, Smithsonian Physical Tables, 9th edn. (Norwich, New York: Knovel) [Google Scholar]
  25. Freytag, B., Liljegren, S., & Höfner, S. 2017, A&A, 600, A137 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  26. Gail, H.-P., & Sedlmayr, E. 1985, A&A, 148, 183 [Google Scholar]
  27. Gail, H. P., & Sedlmayr, E. 1987a, A&A, 171, 197 [NASA ADS] [Google Scholar]
  28. Gail, H.-P., & Sedlmayr, E. 1987b, A&A, 177, 186 [NASA ADS] [Google Scholar]
  29. Gail, H.-P., & Sedlmayr, E. 1988, A&A, 206, 153 [NASA ADS] [Google Scholar]
  30. Gail, H. P., & Sedlmayr, E. 2014, in Physics and Chemistry of Circumstellar Dust Shells, (New York: Cambridge University Press), Cambridge Astrophysics Series, 52 [Google Scholar]
  31. Gail, H.-P., Keller, R., & Sedlmayr, E. 1984, A&A, 133, 320 [NASA ADS] [Google Scholar]
  32. Gauger, A., Sedlmayr, E., & Gail, H.-P. 1990, A&A, 235, 345 [Google Scholar]
  33. Gnedin, N. Y., & Abel, T. 2001, New A, 6, 437 [NASA ADS] [CrossRef] [Google Scholar]
  34. Grassi, T., Bovino, S., Schleicher, D. R. G., et al. 2014, MNRAS, 439, 2386 [Google Scholar]
  35. Hasegawa, K., & Umemura, M. 2010, MNRAS, 407, 2632 [NASA ADS] [CrossRef] [Google Scholar]
  36. Haworth, T. J., Ilee, J. D., Forgan, D. H., et al. 2016, PASA, 33, e053 [NASA ADS] [CrossRef] [Google Scholar]
  37. Höfner, S. 2008, A&A, 491, L1 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  38. Höfner, S., & Dorfi, E. A. 1992, A&A, 265, 207 [Google Scholar]
  39. Höfner, S., & Freytag, B. 2019, A&A, 623, A158 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  40. Höfner, S., Feuchtinger, M. U., & Dorfi, E. A. 1995, A&A, 297, 815 [Google Scholar]
  41. Höfner, S., Gautschy-Loidl, R., Aringer, B., & Jørgensen, U. G. 2003, A&A, 399, 589 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  42. Homan, W., Montargès, M., Pimpanuwat, B., et al. 2020, A&A, 644, A61 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  43. Homan, W., Pimpanuwat, B., Herpin, F., et al. 2021, A&A, 651, A82 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  44. Hutchison, M., Price, D. J., & Laibe, G. 2018, MNRAS, 476, 2186 [Google Scholar]
  45. Iliev, I. T., Ciardi, B., Alvarez, M. A., et al. 2006, MNRAS, 371, 1057 [NASA ADS] [CrossRef] [Google Scholar]
  46. Iliev, I. T., Whalen, D., Mellema, G., et al. 2009, MNRAS, 400, 1283 [NASA ADS] [CrossRef] [Google Scholar]
  47. Jahanara, B., Mitsumoto, M., Oka, K., et al. 2005, A&A, 441, 589 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  48. Jeong, K. S., Winters, J. M., Le Bertre, T., & Sedlmayr, E. 2003, A&A, 407, 191 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  49. Jorissen, A. 2003, in Asymptotic giant branch stars, eds. H. J. Habing, & H. Olofsson (Berlin: Springer), Astron. Astrophys. Lib., 461 [Google Scholar]
  50. Jorissen, A., & Knapp, G. R. 1998, A&AS, 129, 363 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  51. Kervella, P., Homan, W., Richards, A. M. S., et al. 2016, A&A, 596, A92 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  52. Kessel-Deynet, O., & Burkert, A. 2000, MNRAS, 315, 713 [NASA ADS] [CrossRef] [Google Scholar]
  53. Kim, H., & Taam, R. E. 2012, ApJ, 759, L22 [NASA ADS] [CrossRef] [Google Scholar]
  54. Krüger, D., & Sedlmayr, E. 1997, A&A, 321, 557 [NASA ADS] [Google Scholar]
  55. Krüger, D., Gauger, A., & Sedlmayr, E. 1994, A&A, 290, 573 [NASA ADS] [Google Scholar]
  56. Laibe, G., & Price, D. J. 2012, MNRAS, 420, 2345 [NASA ADS] [CrossRef] [Google Scholar]
  57. Laibe, G., & Price, D. J. 2014a, MNRAS, 440, 2147 [NASA ADS] [CrossRef] [Google Scholar]
  58. Laibe, G., & Price, D. J. 2014b, MNRAS, 444, 1940 [Google Scholar]
  59. Lamers, H. J. G. L. M., & Cassinelli, J. P. 1999, Introduction to Stellar Winds (Cambridge University Press), 452 [Google Scholar]
  60. Liu, Z.-W., Stancliffe, R. J., Abate, C., & Matrozis, E. 2017, ApJ, 846, 117 [Google Scholar]
  61. Lodato, G., & Price, D. J. 2010, MNRAS, 405, 1212 [NASA ADS] [Google Scholar]
  62. MacLeod, M., & Loeb, A. 2020, ApJ, 902, 85 [NASA ADS] [CrossRef] [Google Scholar]
  63. Maercker, M., Mohamed, S., Vlemmings, W. H. T., et al. 2012, Nature, 490, 232 [NASA ADS] [CrossRef] [Google Scholar]
  64. Maes, S., Homan, W., Malfait, J., et al. 2021, A&A, 653, A25 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  65. Malfait, J., Homan, W., Maes, S., et al. 2021, A&A, 652, A51 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  66. Mastrodemos, N., & Morris, M. 1998, ApJ, 497, 303 [Google Scholar]
  67. Mattsson, L. 2016, Rev. Mex. Astron. Astrofis., 87, 249 [Google Scholar]
  68. Mattsson, L., & Höfner, S. 2011, A&A, 533, A42 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  69. Mattsson, L., & Sandin, C. 2021, Universe, 7, 113 [NASA ADS] [CrossRef] [Google Scholar]
  70. Mattsson, L., Wahlin, R., & Höfner, S. 2010, A&A, 509, A14 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  71. Mauron, N., & Huggins, P. J. 2006, A&A, 452, 257 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  72. Mohamed, S., & Podsiadlowski, P. 2012, Balt. Astron., 21, 88 [Google Scholar]
  73. Neugebauer, M., & Snyder, C. W. 1962, Science, 138, 1095 [Google Scholar]
  74. Pawlik, A. H., & Schaye, J. 2008, MNRAS, 389, 651 [NASA ADS] [CrossRef] [Google Scholar]
  75. Petkova, M., & Springel, V. 2009, MNRAS, 396, 1383 [NASA ADS] [CrossRef] [Google Scholar]
  76. Petkova, M., & Springel, V. 2011, MNRAS, 415, 3731 [NASA ADS] [CrossRef] [Google Scholar]
  77. Petkova, M. A., Vandenbroucke, B., Bonnell, I. A., & Kruijssen, J. M. D. 2021, MNRAS, 507, 858 [NASA ADS] [CrossRef] [Google Scholar]
  78. Pinte, C., van der Plas, G., Ménard, F., et al. 2019, Nat. Astron., 3, 1109 [NASA ADS] [CrossRef] [Google Scholar]
  79. Preibisch, T., Ossenkopf, V., Yorke, H. W., & Henning, T. 1993, A&A, 279, 577 [NASA ADS] [Google Scholar]
  80. Price, D. J. 2012, J. Comput. Phys., 231, 759 [NASA ADS] [CrossRef] [Google Scholar]
  81. Price, D. J., & Federrath, C. 2010, MNRAS, 406, 1659 [NASA ADS] [Google Scholar]
  82. Price, D. J., & Laibe, G. 2015, MNRAS, 451, 813 [NASA ADS] [CrossRef] [Google Scholar]
  83. Price, D. J., Wurster, J., Tricco, T. S., et al. 2018, PASA, 35, e031 [Google Scholar]
  84. Rijkhorst, E. J., Plewa, T., Dubey, A., & Mellema, G. 2006, A&A, 452, 907 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  85. Sandin, C. 2008, MNRAS, 385, 215 [NASA ADS] [CrossRef] [Google Scholar]
  86. Sandin, C., & Höfner, S. 2004, A&A, 413, 789 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  87. Sandin, C., & Mattsson, L. 2020, MNRAS, 499, 1531 [NASA ADS] [CrossRef] [Google Scholar]
  88. Sharp, C. M., & Huebner, W. F. 1990, ApJS, 72, 417 [NASA ADS] [CrossRef] [Google Scholar]
  89. Skinner, M. A., & Ostriker, E. C. 2013, ApJS, 206, 21 [NASA ADS] [CrossRef] [Google Scholar]
  90. Struck, C., Smith, D. C., Willson, L. A., Turner, G., & Bowen, G. H. 2004, MNRAS, 353, 559 [NASA ADS] [CrossRef] [Google Scholar]
  91. Susa, H. 2006, PASJ, 58, 445 [NASA ADS] [Google Scholar]
  92. Tabak, R. G., Hirth, J. P., Meyrick, G., & Roark, T. P. 1975, ApJ, 196, 457 [NASA ADS] [CrossRef] [Google Scholar]
  93. Tegmark, M. 1996, ApJ, 470, L81 [NASA ADS] [CrossRef] [Google Scholar]
  94. Theuns, T., & Jorissen, A. 1993, MNRAS, 265, 946 [Google Scholar]
  95. Thorn, R. J., & Winslow, G. H. 1957, J. Chem. Phys., 26, 186 [NASA ADS] [CrossRef] [Google Scholar]
  96. Tsuji, T. 1973, A&A, 23, 411 [NASA ADS] [Google Scholar]
  97. van der Helm, E., Saladino, M. I., Portegies Zwart, S., & Pols, O. 2019, A&A, 625, A85 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  98. Velli, M. 2001, Ap&SS, 277, 157 [NASA ADS] [CrossRef] [Google Scholar]
  99. Wang, N., & Lee, J. 2011, SIAM J. Sci. Comput., 33, 2536 [NASA ADS] [CrossRef] [Google Scholar]
  100. Whitehouse, S. C., Bate, M. R., & Monaghan, J. J. 2005, MNRAS, 364, 1367 [NASA ADS] [CrossRef] [Google Scholar]
  101. Winters, J. M., Le Bertre, T., Jeong, K. S., Helling, C., & Sedlmayr, E. 2000, A&A, 361, 641 [NASA ADS] [Google Scholar]

Appendix A: Chemical network

The chemical network used to compute the carbon abundance needed to estimate the opacity includes the following 26 reactions :

H + H H 2 O + H OH O + 2 H H 2 O C + C C 2 C + O CO C + 2 O CO 2 C + 4 H CH 4 C + H C 2 H H + H C 2 H 2 N + N N 2 N + 3 H NH 3 C + N CN CN + H HCN Si + Si Si 2 Si + 2 Si Si 3 Si + O SiO 2 Si + C Si 2 C Si + 4 H SiH 4 S + S S 2 S + H HS S + 2 H H 2 S Si + S SiS Si + H SiH Ti + O TiO Ti + 2 O TiO 2 Ti + S TiS $$ \begin{aligned} \mathrm{H} + \mathrm{H}&\rightarrow \mathrm{H} _2 \\ \mathrm{O} + \mathrm{H}&\rightarrow \mathrm{OH} \\ \mathrm{O} + 2\mathrm{H}&\rightarrow \mathrm{H} _2\mathrm{O} \\ \mathrm{C} + \mathrm{C}&\rightarrow \mathrm{C} _2 \\ \mathrm{C} + \mathrm{O}&\rightarrow \mathrm{CO} \\ \mathrm{C} + 2\mathrm{O}&\rightarrow \mathrm{CO} _2 \\ \mathrm{C} + 4\mathrm{H}&\rightarrow \mathrm{CH} _4 \\ \mathrm{C} + \mathrm{H}&\rightarrow \mathrm{C} _2\mathrm{H} \\ \mathrm{H} + \mathrm{H}&\rightarrow \mathrm{C} _2\mathrm{H} _2 \\ \mathrm{N} + \mathrm{N}&\rightarrow \mathrm{N} _2 \\ \mathrm{N} + 3\mathrm{H}&\rightarrow \mathrm{NH} _3 \\ \mathrm{C} + \mathrm{N}&\rightarrow \mathrm{CN} \\ \mathrm{CN} + \mathrm{H}&\rightarrow \mathrm{HCN} \\ \mathrm{Si} + \mathrm{Si}&\rightarrow \mathrm{Si} _2 \\ \mathrm{Si} + 2\mathrm{Si}&\rightarrow \mathrm{Si} _3 \\ \mathrm{Si} + \mathrm{O}&\rightarrow \mathrm{SiO} \\ 2\mathrm{Si} + \mathrm{C}&\rightarrow \mathrm{Si} _2\mathrm{C} \\ \mathrm{Si} + 4\mathrm{H}&\rightarrow \mathrm{SiH} _4 \\ \mathrm{S} + \mathrm{S}&\rightarrow \mathrm{S} _2 \\ \mathrm{S} + \mathrm{H}&\rightarrow \mathrm{HS} \\ \mathrm{S} + 2\mathrm{H}&\rightarrow \mathrm{H} _2\mathrm{S} \\ \mathrm{Si} + \mathrm{S}&\rightarrow \mathrm{SiS} \\ \mathrm{Si} + \mathrm{H}&\rightarrow \mathrm{SiH} \\ \mathrm{Ti} + \mathrm{O}&\rightarrow \mathrm{TiO} \\ \mathrm{Ti} + 2\mathrm{O}&\rightarrow \mathrm{TiO} _2 \\ \mathrm{Ti} + \mathrm{S}&\rightarrow \mathrm{TiS} \end{aligned} $$

thumbnail Fig. A.1.

Evolution of the number density of the most abundant species as a function of temperature for the envelope composition given in Table 1 and for a density ρ = 2 × 10−14 g cm−3 and a carbon to oxygen ratio C/O=2.

The dissociation constants for the different reactions are calculating using Eqs. (5) and (6) where the Gibbs energies are taken from the fits of the JANAF tables (Chase 1986). They are expressed in terms of polynomials of the temperature T of the form

Δ G f = a T + b + c T + d T 2 + e T 3 , $$ \begin{aligned} \Delta G_\mathrm{f} = \frac{a}{T}+b+cT+dT^2+eT^3, \end{aligned} $$(A.1)

where the coefficients a,  b,  c,  d and e can be found in Sharp & Huebner (1990). For TiS, the fit is taken from Tsuji (1973).

Figure A.1 illustrates how the abundance of the main species vary as a function of temperature. At low temperatures (T ≲ 1400 K) all hydrogen is in the form of H2 and the second most abundant species is CO. As the temperature increases, H atoms become available and C2H2 dissociates. At higher temperature (T >  2500 K) all the material returns to an atomic state. We note that CO is the most resistant molecule and at this selected density it can survive up to T ≈ 3000 K.

Appendix B: Computing the net growth rate

With the four chemical species C, C2, C2H, and C2H2 that we consider for the evolution of carbon grains, the net growth rate (Eq. 17) writes

1 τ = 1 τ C + 1 τ C 2 + 1 τ C 2 H + 1 τ C 2 H 2 . $$ \begin{aligned} \frac{1}{\tau } = \frac{1}{\tau _\mathrm{C} } + \frac{1}{\tau _{\mathrm{C} _2}} + \frac{1}{\tau _{\mathrm{C} _2\mathrm{H} }} + \frac{1}{\tau _{\mathrm{C} _2\mathrm{H} _2}}. \end{aligned} $$(B.1)

The first two terms of Eq. (17), which include PC and PC2, account for homogeneous grain growth and evaporation. In chemical equilibrium, the general form of the second term reads (Gauger et al. 1990),

1 τ C 2 = 2 A 1 α ( 2 ) P C 2 2 π k T g m C 2 [ 1 α ( 2 ) S 2 b C 2 ] , $$ \begin{aligned} \frac{1}{\tau _{\mathrm{C} _2}} = \frac{2 \mathcal{A} _1 \alpha (2) P_{\mathrm{C} _2}}{\sqrt{2 \pi k T_{\mathrm{g} }\, m_{\mathrm{C} _2}}} \left[1 - \frac{\alpha _{*}(2)}{{S}^2 b_{\mathrm{C} _2}}\right], \end{aligned} $$(B.2)

where α(2) is an averaged sticking efficiency for a dimer, in this case for C2, and the factor α*(2) accounts for non-thermal equilibrium (non-TE) effects,

α ( 2 ) = α ˚ ( 2 ) α ( 2 ) . $$ \begin{aligned} \alpha _{*}(2) = \frac{\mathring{\alpha }(2)}{\alpha (2)}. \end{aligned} $$(B.3)

The quantity α ˚ $ \mathring{\alpha} $ denotes the averaged thermal equilibrium (TE) reaction efficiency when the gas temperature equals the grain temperature. We assume TE for the reaction efficiency and use the sticking coefficient α2 (see sect. 2.3.1) both for α ˚ ( 2 ) $ \mathring{\alpha}(2) $ and α(2). Similarly, we use α1 in the expression of τ C 1 $ \tau_\mathrm{C}^{-1} $. In chemical equilibrium, and when the gas has the same temperature as the dust, the correction factors b read bC = bC2 = 1.

The third and fourth terms of Eq. (17) that include the contributions from the molecules C2H and C2H2, account for heteromolecular growth and chemical sputtering. The general form of these terms has been presented by Gauger et al. (1990), and for the third term it reads

1 τ C 2 H = 2 A 1 α C 2 H c ( 2 ) P C 2 H 2 π k T g m C 2 H [ 1 α c ( 2 , C 2 H ) S 2 b C 2 H c ] , $$ \begin{aligned} \frac{1}{\tau _{\mathrm{C} _2\mathrm{H} }} = \frac{2 \mathcal{A} _1 \alpha _{\mathrm{C} _2\mathrm{H} }^\mathrm{c} (2) P_{\mathrm{C} _2\mathrm{H} }}{\sqrt{2 \pi k T_{\mathrm{g} }\, m_{\mathrm{C} _2\mathrm{H} }}} \left[1 - \frac{\alpha _{*}^\mathrm{c} (2, \mathrm{C} _2\mathrm{H} )}{{S}^2 b_{\mathrm{C} _2\mathrm{H} }^\mathrm{c} }\right], \end{aligned} $$(B.4)

where the factor α c ( 2 , C 2 H ) $ \alpha_{*}^\mathrm{c}(2, \mathrm{C}_2\mathrm{H}) $ expresses non-TE effects,

α c ( 2 , C 2 H ) = α ˚ C 2 H c ( 2 ) β C 2 H c ( 2 ) α C 2 H c ( 2 ) β ˚ C 2 H c ( 2 ) . $$ \begin{aligned} \alpha _{*}^\mathrm{c} (2, \mathrm{C} _2\mathrm{H} ) = \frac{\mathring{\alpha }_{\mathrm{C} _2\mathrm{H} }^\mathrm{c} (2) \beta _{\mathrm{C} _2\mathrm{H} }^\mathrm{c} (2)}{\alpha _{\mathrm{C} _2\mathrm{H} }^\mathrm{c} (2) \mathring{\beta }_{\mathrm{C} _2\mathrm{H} }^\mathrm{c} (2)}. \end{aligned} $$(B.5)

α denotes the reaction efficiency of the growth reaction, and β the efficiency of the reverse reaction, β ˚ $ \mathring{\beta} $ being the TE value of the latter. We assume that both are equal (αc = βc) and set α c ( 2 , C 2 H ) = α c ( 2 , C 2 H 2 ) = 1 $ \alpha_{*}^\mathrm{c}(2, \mathrm{C}_2\mathrm{H}) = \alpha_{*}^\mathrm{c}(2, \mathrm{C}_2\mathrm{H}_2) = 1 $. Moreover, we set α C 2 H c ( 2 ) = α C 2 H 2 c ( 2 ) = α 2 $ \alpha_{\mathrm{C}_2\mathrm{H}}^\mathrm{c}(2) = \alpha_{\mathrm{C}_2\mathrm{H}_2}^\mathrm{c}(2) = \alpha_2 $.

All Tables

Table 1.

Elemental abundances by number relative to hydrogen ϵ used in our models (Ferrarotti & Gail 2002).

Table 2.

Optimized rotation angles (in rad) for the distribution of SPH particles at the inner boundary as a function of the number of vertices or particles (Np = 10(2q − 1)2 + 2).

Table 3.

1D model parameters. Tiso corresponds to the temperature of the isothermal wind.

All Figures

thumbnail Fig. 1.

Nucleation rate per hydrogen atom ( J s / n H $ J_{*}^{\mathrm{s}}/n_{\langle\mathrm{H}\rangle} $) as a function of temperature and pressure for a C/O ratio of 10. The computed values are nearly identical to the results presented in Fig. 5 of Gail et al. (1984).

In the text
thumbnail Fig. 2.

Planck mean values of the gas opacity κg (logarithmic scale), derived from a grid of models calculated with the Cloudy code using solar abundances with C/O = 1.4, and Jν = Bν(3000 K). The thick red line traces a wind particle and illustrates the range of densities and temperatures encountered in our models.

In the text
thumbnail Fig. 3.

Planck mean values of the dust extinction efficiency Qext for amorphous carbon grains with size a = 0.01, 0.1 and 1 μm (black solid lines, from Preibisch et al. 1993). For comparison, the approximations by Draine (1981; blue) and Winters et al. (2000) with n = 1.19 (magenta), 1.38 (green) and 1.46 (red) are shown.

In the text
thumbnail Fig. 4.

Schematic representation of the method for injecting SPH particles. The horizontal axes represent the time since ejection, or interchangeably the distance from the boundary. The numbered dots represent particle shells, numbered by their index. Older shells have a lower index. Top: before the hydro timestep the shell with index 12 lies below the injection boundary, and is therefore part of the handled shells, and has imposed physical properties. Bottom: after the hydro timestep, the shell with index 12 has moved beyond the injection boundary, and no longer has imposed physical properties.

In the text
thumbnail Fig. 5.

Isothermal velocity profiles of a transonic (bottom) and supersonic (top) SPH wind solution, in red. The analytic solution is given by the black curves for a wind temperature Tiso = 1.25 × 104 K.

In the text
thumbnail Fig. 6.

Adiabatic velocity profiles of a transonic (bottom) and supersonic (top) SPH wind solution, in red. The analytic solution is given by the black curves.

In the text
thumbnail Fig. 7.

Effect of the resolution setup on the predicted dynamics of an transonic adiabatic wind, for a fixed value of wss. The value of q is chosen such that SPH particle mass is decreased by a factor ten from the top to the middle panel, and again from the middle to the bottom panel.

In the text
thumbnail Fig. 8.

Acceleration profile of a free, steady, polytropic wind with the Bowen dust opacity prescription activated around ∼ 2 au. The temperature is shown in blue.

In the text
thumbnail Fig. 9.

Acceleration profile of a free wind for which the dust opacity is calculated using the nucleation formalism.

In the text
thumbnail Fig. 10.

Radial profiles of chemical and dust-related quantities. The temperature of the SPH gas, with which the dust is assumed to be fully thermalised, is shown in red. Left: the logarithm of the third moment 𝒦3 (blue), which represents the number density of condensing monomers from the gas phase into the dust, shown together with the dust opacity profile (green). Right: the mean molecular weight μ (blue), polytropic index γ (green), and the base 10 logarithm of the supersaturation ratio S (magenta) of the gas in the wind. A supersaturation ratio above unity implies favourable conditions for dust formation.

In the text
thumbnail Fig. A.1.

Evolution of the number density of the most abundant species as a function of temperature for the envelope composition given in Table 1 and for a density ρ = 2 × 10−14 g cm−3 and a carbon to oxygen ratio C/O=2.

In the text

Current usage metrics show cumulative count of Article Views (full-text article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.

Data correspond to usage on the plateform after 2015. The current usage metrics is available 48-96 hours after online publication and is updated daily on week days.

Initial download of the metrics may take a while.