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/00046361/202243540  
Published online  09 November 2022 
3D simulations of AGB stellar winds
I. Steady winds and dust formation
^{1}
Institut d’Astronomie et d’Astrophysique, Université libre de Bruxelles (ULB), CP 226, 1050 Brussels, Belgium
email: lionel.siess@ulb.be
^{2}
School of Physics and Astronomy, Monash University, Clayton, VIC 3800, Australia
Received:
14
March
2022
Accepted:
12
August
2022
Aims. We present the implementation of the treatment of particle ejection and dust nucleation in the smoothed particle hydrodynamics (SPH) code PHANTOM. These developments represent the first step toward a more complete modeling of dustdriven winds emanating from asymptotic giant branch (AGB) stars that can be used for comparison with high resolution imaging of these stars.
Methods. The AGB outflow is modeled by injecting the SPH particles from a spherical inner boundary. This boundary is a series of concentric shells, with the AGB star at its center, and the particles are positioned on these shells on the vertices of an isocahedron geodesic surface. The outermost shell is ejected with a predefined radial velocity, and subsequent lower shells replenish the ejected ones, all rotated randomly to improve the isotropy of the outflow. The physical properties of the particles on these shells are set by solving the 1D analytic steady wind equations. The formation of dust is calculated starting from a compact chemical network for carbonrich material, which creates the building blocks of the solidstate particles. Subsequently, the theory of the moments is used to obtain dust growth rates, without requiring knowledge on the grain size distribution.
Results. We tested our implementation against a series of 1D reference solutions. We demonstrate that our method is able to reproduce Parkertype wind solutions. For the transsonic solution, small oscillations are present in the vicinity of the sonic point, but these do not impact the transsonic passage or terminal wind velocity. Supersonic solutions always compare nicely with 1D analytic profiles. We also tested our implementation of dust using two formalisms: an analytic prescription for the opacity devised by Bowen and the full treatment of carbondust formation. Both simulations reproduce the 1D analytic solution displaying the expected additional acceleration when the gas temperature falls below the condensation temperature.
Key words: stars: winds, outflows / methods: numerical / hydrodynamics / stars: AGB and postAGB / astrochemistry
© L. Siess et al. 2022
Open 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 SubscribetoOpen model. Subscribe to A&A to support open access publication.
1. Introduction
During the thermally pulsing asymptotic giant branch (AGB) phase, low and intermediatemass stars develop strong mass loss and eject their envelope to become a shortlived (≲10^{4} yr) postAGB 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 byproducts 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 largescale 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 lowmass 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 ValBorro et al. 2017; Chen et al. 2018; El Mellah et al. 2020; MacLeod & Loeb 2020; BermúdezBustamante 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 starinabox 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 dustdriven 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 largescale 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 largescale 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 dustdriven 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 carbonrich mixtures (C/O > 1), whose chemistry is better known and has been extensively studied in the context of latetype 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 companionperturbed 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 gasdust drift using the “onefluid approximation” of hydrodynamics for small grains available in the code (Laibe & Price 2014b; Price & Laibe 2015), the radiative transfer and the acceleration of the dustcomponent 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 dustdriven winds closely follows the description presented in the reference textbook of Gail & Sedlmayr (2014). We consider a singlefluid approach where the grains are at rest with the gas particles and focus on mass loss from carbonrich 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 ≤ N_{0}. Grain growth thus proceeds through the addition of small clusters of size N ≤ N_{0}. In our calculations of Cbased dust, the growth of graphite is assumed to proceed homogeneously by the addition of clusters of the same type (C_{i}molecules where we consider i = 1, 2)
and by the reactions involving Cbased molecules of a different type (heteromolecular growth)
This set of reactions corresponds to N_{0} = 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 N_{l} of a dust grains. This quantity is not easy to determine and following Gail & Sedlmayr (1988) we will use N_{l} = 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 timedependent problem.
The formation of Cbased dust grains depends on the amount of atomic carbon and carboncontaining 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 C_{2}, C_{2}H, and C_{2}H_{2} 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 CH_{4} (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 Cbased dust, a network of chemical reactions must be considered. In our code, we take into account 26 molecules H_{2}, OH, H_{2}O, CO, CO_{2}, CH_{4}, C_{2}H, C_{2}H_{2}, N_{2}, NH_{3}, CN, HCN, Si_{2}, Si_{3}, SiO, Si_{2}C, SiH_{4}, S_{2}, HS, H_{2}S, SiS, SiH, TiO, TiS, TiO_{2} 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 H_{2}, the law of mass action of the reaction 2 H → H_{2} in chemical equilibrium is used, which reads
(cf. Eq. (10.16) of Gail & Sedlmayr 2014), where P_{H2} and P_{H} are the partial pressures of H_{2} and H, respectively and K_{H2} is the dissociation constant of the reaction. The latter is given by
where ℛ is the universal gas constant, T_{g} the gas temperature, and ΔG_{H2} is the change in Gibbs energy of the reaction,
We take the values of the Gibbs energies of formation ΔG_{f} 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
assuming that H_{2} is the dominant hydrogenous molecule. Similarly, the total number density of nitrogen atoms n_{⟨N⟩} = ϵ_{N}n_{⟨H⟩}, where ϵ_{N} is the nitrogen abundance by number relative to hydrogen, reads
where K_{NH3}, K_{HCN}, and K_{N2} are the dissociation constants of the reactions of formation of NH_{3}, HCN, and N_{2}, respectively. The constants are calculated by use of Eq. (5), substituting ΔG_{H2} 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 P_{el} of the elements “el”. These atomic partial pressures are then used to calculate the molecular partial pressures P_{mol} of the molecules “mol” or, equivalently, their number densities n_{mol} = P_{mol}/(kT_{g}), 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.
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
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 N_{l} ∼ 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
where a_{0} is the radius of a monomer inside a grain, which is derived from the density ρ_{carb} of the grain material, i.e.
where A_{carb} = 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 a_{0} = 1.28 × 10^{−8} cm. Global dust properties can be estimated from the knowledge of the moments. In particular,

a)
the mean grain density n_{d} = 𝒦_{0}, that is the density of all dust grains with size N ≥ N_{l}.

b)
the average grain radius ⟨a⟩=a_{0}𝒦_{1}/𝒦_{0}.

c)
the average grain surface .

d)
the average particle size ⟨N⟩=𝒦_{3}/𝒦_{0}.

e)
the number density of monomers of size N ≥ N_{l} condensed into grains n_{cond} = 𝒦_{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
where ρ is the density and the mean mass per Hatom calculated using the abundance by number relative to hydrogen ϵ_{i}, atomic weight A_{i} of species i (see Table 1 for a list of the selected species) and m_{H} 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
where is the initial abundance at the beginning of the nucleation. With these normalized variables, the evolution of the first four moments reads
where d/dt is the material derivative. τ^{−1} is the net rate of growth or evaporation of the grains (see next section), is the normalized rate of formation of critical clusters per unit volume, the normalized quasistationary 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 > N_{l}. 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 quasistationary rate is constant. We stress that Δt is constrained so this assumption holds.
2.3.1. Net rate of grain growth and destruction
For our Crich chemistry, as described by Gauger et al. (1990), the net rate of grain growth and destruction can be expressed as
where we only consider the chemical reactions given by Eqs. (1)–(3) involving C, C_{2}, C_{2}H and C_{2}H_{2}. The average relative velocity v_{rel, i} between the reacting species i (mass m_{i}) and the dust particle, in the absence of drift, is given by
The quantities α_{1} = 0.37 and α_{2} = 0.34 are the sticking coefficients of C and C_{2} (Thorn & Winslow 1957) and the reaction efficiency of the considered molecules (C_{2}H or C_{2}H_{2}) with carbon. Following Gail & Sedlmayr (1988), we set . The coefficients α_{*}(i) quantify the deviation from a thermal equilibrium population of states and are given by α_{*}(i)=(T_{d}/T_{g})^{1/2}, where T_{d} and T_{g} refer to the dust and gas temperatures, respectively. The quantity b_{i} accounts for departures of the particle densities from chemical and thermodynamical equilibrium between the solid (dust) and gas phase. Assuming chemical equilibrium, the coefficients b_{i} are given by
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
where the saturation vapour pressure of carbon in the solid phase
is derived using the dust temperature (T = T_{d}). The free enthalpies of formation are taken from the JANAF tables and according to Eq. (6) reads ΔG_{C, solid} = ΔG_{f}(C, solid)−ΔG_{f}(C). See Appendix B for a more detailed explanation of Eq. (17). When gas and dust are thermally coupled (T_{g} = T_{d}), which is the assumption adopted in the current paper, the expressions for the growth rate simplifies since in this case α_{*}(i)=b_{i} = 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 gasdust 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)
The quantity is the density of grains of size N in thermal equilibrium,
where ΔF is the free energy of formation of the grain from the vapour (Draine & Salpeter 1977). The saturation ratio of the gas phase is computed using T = T_{g} in Eq. (24), n_{C} 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,
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 N_{h} 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 N_{h} = 5.
The size of the critical cluster N_{*} is determined by the root of the equation . Its value is given by
where
The Zeldovich factor Z entering Eq. (25) measures the deviation of the grain distribution from equilibrium and is defined by
and β, which gives the flux of monomers colliding with the grains, by
(Gail et al. 1984). With these quantities, the relaxation time towards equilibrium writes
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 H_{2}, C_{2}H_{2}, 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.
Fig. 1. Nucleation rate per hydrogen atom () 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 dustdriven 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)
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 m_{i}, its position r_{i}, smoothing length h_{i}, velocity v_{i}, pressure P_{i}, artificial viscosity q_{i}, and specific internal energy e_{i}. W(r_{i} − r_{j}, h) is the smoothing kernel, ∇_{i}W_{ij} its gradient with respect to the location r_{i} 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 dustdriven stellar wind, the outward force exerted on the gas is generated by a twostage 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 coupling^{1}), it can be shown that
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 massloss rate is higher than ∼10^{−8} M_{⊙} 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.,
where P is the total gas pressure, and m_{u} 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 H_{2}, giving
and following Grassi et al. (2014)
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 H_{2} 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 H_{2} 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 supersaturation ratio S.
We neglect the dust component in the EOS because the dusttogas 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 r_{i}, v_{i}, e_{i} and h_{i} (ρ_{i} is a simple function of h_{i}). 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 H_{2}. These instabilities were eliminated after we introduced a subloop over μ_{i}, γ_{i} and T to iteratively obtain continuous selfconsistent 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} cm^{2} 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} cm^{2} 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 n_{H} − 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.
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
where n(a) is the grain size distribution function (cm^{−4}) and Q_{ext} 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 socalled small particle limit (SPL), the extinction efficiency is dominated by absorption and in this Rayleigh regime it can be expressed as allowing to separate in Eq. (41) the wavelength and grain size dependence leading to
where a_{0} 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 nongrey 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 radiativetogravitational 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 (Orich) Mtype 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)
Figure 3 shows the Planck mean absorption efficiencies (Q_{ext}) of amorphous carbon grains from various sources. The Q_{ext}’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 Draine^{2}). 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 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.
Fig. 3. Planck mean values of the dust extinction efficiency Q_{ext} 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 R_{inj}. 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
where c_{s} is the adiabatic sound speed defined by . Assuming a stationary wind, Eq. (44) reduces to the equation of mass conservation
where Ṁ is the constant massloss rate. In this case, DP/Dt = v dP/dr and Eq. (46) can be used to replace the pressure gradient in Eq. (45), giving
The ideal gas law (Eq. (37)) can be used to convert Eq. (46) into an equation for the temperature
Finally the time evolution is given by
We integrate Eqs. (48)–(50) using the fourthorder RungeKutta method, to obtain the physical properties of the ejected SPH particles as a function of both distance and time. For transsonic 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 transsonic 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 breezetype 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 < c_{s}, that is when the wind is subsonic, the breeze is obtained as long as the outward acceleration factor Γ is higher than the threshold value
where v_{esc} = (2GM_{*}(1 − Γ)/R_{*})^{1/2} is the escape velocity.
Our SPH simulations are unable to reproduce breezetype solutions. Even if particles are launched from a hot corona with a velocity below the transsonic value, they ultimately converge on the transsonic 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 nonzero outer pressure could be implemented using a method similar to the “handledshells” 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 R_{inj}. Taking a constant number p of particles per shell, the time interval δt between two subsequent shells is given by the relation
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 . 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 t^{n + 1} is the time to be reached and Δt^{n} the value of the current hydro (SPH) time step, then the shells with index i_{min} to i_{max} have to be ejected, where
(“⌊ ⌋” is the floor function). When Δt < δt, i_{min} may become greater than i_{max}, in which case no shell is released.
However, in general the hydrodynamical time t^{n + 1} does not coincide with the precise time t̃_{i} at which injection should occur. Consequently, if t^{n+1} > t̃_{i} the properties of the released particles have to be synchronized and estimated at the “current” time, that is at a time after their ejection. The velocity is given by the solution v(t_{e}) of Eq. (48), and the density by Eq. (47), converted to ρ(t_{e}) 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(t_{e}) 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 R_{inj}. 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.
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 shell^{3}. 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.
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 (N_{p} = 10(2q − 1)^{2} + 2).
where c_{i} = cosϑ_{i} and s_{i} = sinϑ_{i}. 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 .
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
where is the golden ratio, and q can be any natural number, which represents the tangential resolution setup, 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
where w_{ss} 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 lengthscale of the system must be resolved. We address how to approach both these requirements in wind systems.
The smoothing length h is given by h_{fact} n^{−1/3} (Price et al. 2018), where n is the particle number density, and h_{fact} 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 rewritten as
where d_{⊥} is the tangential particle distance given by Eq. (56), and d_{r} ≈ δr = w_{ss}d_{⊥} is the radial distance between particle on adjacent shells. Hence, the first resolution requirement is satisfied when both the h ≥ d_{r} and h ≥ d_{⊥} conditions are simultaneously fulfilled. Because d_{r} depends on d_{⊥}, this translates to a resolution condition on w_{ss} given by
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 lengthscale of the system. For a Parker wind (Lamers & Cassinelli 1999), this lengthscale can be identified to the sonic radius. For a windRochelobeoverflow binary (Mohamed & Podsiadlowski 2012), the representative lengthscale 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 massloss 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
where v_{inj} is the wind velocity at the injection radius R_{inj} (also available in the PHANTOM parameter card as wind_velocity), and N_{p} is the number of particles^{4} of mass m_{p}. Given that the intershell distance is also given by Eq. (57), this implies that the mass of the particles is determined by the setup of the boundary conditions and
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 abovementioned 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 dustfree 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 w_{ss} = 1.0. Finally, to position the sonic radius beyond the radius of the handled shells an initial sound speed c_{s, iso} = 13.12 km s^{−1} was assumed (corresponding to T_{iso} = 1.25 × 10^{4} K). In all cases presented here and in the following sections, conditions resulting in a breezetype solution (see Sect. 4.1) were deliberately avoided because the SPH solver was unable to recover it properly.
1D model parameters. T_{iso} 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 sonicpoint 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 quasitransonic 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.
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 T_{iso} = 1.25 × 10^{4} 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 supersonic 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 sonicpoint oscillations can be seen for the subsonic portion of the transsonic simulation.
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 lengthscale 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.
Fig. 7. Effect of the resolution setup on the predicted dynamics of an transonic adiabatic wind, for a fixed value of w_{ss}. 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
where κ_{max} is the maximum value the dust opacity can reach, δ is a measure for the temperature range over which condensation occurs, and T_{cond} the dust condensation temperature which depends on the star’s composition (C/O ratio). T_{eq} is the dust equilibrium temperature which value depends on how the gas and dust are thermally coupled. A correct determination of T_{eq} would require radiative transfer calculations which is outside the scope of this paper (see Sect. 6) so for the moment we assume that T_{eq} is set to the gas temperature.
Though the Bowen parameterization is a severe approximation of the opacityforming 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 × 10^{4} L_{⊙} 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 followup 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 adhoc 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
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. 
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 v_{inj} = 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 cm^{2} g^{−1} is chosen, yielding Γ(T_{eq} = T_{cond})≃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 k_{max} 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 RungeKutta 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 dustnucleation described in Sect. 2. The numerical setup is nearly identical to the one used with the Bowen dustopacity 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 rhoinvariant 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 massloss rate of 10^{−5} M_{⊙} yr^{−1}, and an initial velocity of v_{inj} = 6 km s^{−1} was assumed. This initial velocity differs slightly from the Bowentype dustopacitysetup in the previous section to avoid the breezetype solution.
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 dustfree profile (black dashed) in Fig. 9. While the SPH solution reaches a terminal velocity of nearly 30 km s^{−1}, its dustfree counterpart (accelerated only via the pressuregradient at the stellar boundary and α = 1) barely reaches a speed of 12 km s^{−1}. The dusty wind also exhibits a doublestaged acceleration. The innermost velocity profile of the dustfree 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 dustfree curve, and is therefore solely associated with dust formation.
The formation of dust can be quantified either by the rate of condensation from the gasphase, 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 ∼10^{6} 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.
Fig. 10. Radial profiles of chemical and dustrelated 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 gasphase carbon exceeds that of solidstate 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 windinjection boundary, as well as the calculations involving dust opacity can be used to model cool and dusty stellar winds. However, to levelup the modelling of dustdriven winds to match the stateoftheart 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 largescale 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 windcompanion 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., KesselDeynet & 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 protoplanetary 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 dustgas drift. This process is important because conceptually the absence of drift is inconsistent with the dustdriven 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 massloss rates and wind velocities. Fortunately, PHANTOM allows the modelling of dustgas mixtures in either the standard twofluid approach where each species is regarded as a separate fluid (Laibe & Price 2012) or in the socalled onefluid 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 massfraction of the dustgrain mixture. This “diffusion approximation for dust” holds as long as the stopping timestep 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 twofluid 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 stepwise, 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 setup, we have calculated a number of different benchmarking simulations that test the ability of the code to recover the stationary Parkerwind 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 carbonrich chemical network assuming localthermodynamicalequilibrium. 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 (T_{g} = T_{d} 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 carboncompounds 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 setup, and the nucleation and dust growth prescriptions in the PHANTOM code can be used to successfully model dusty stellar outflows.
A comprehensive description of this tiling is provided (in french) at https://fr.wikipedia.org/wiki/Géode_(géométrie)
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 FRSF.N.R.S. research associate. W.H. acknowledges support from a FRSF.N.R.S grant. D.J.P. is grateful for Australian Research Council funding via DP180104235 and FT130100034.
References
 Altay, G., & Theuns, T. 2013, MNRAS, 434, 748 [NASA ADS] [CrossRef] [Google Scholar]
 Andersen, A. C., Höfner, S., & GautschyLoidl, R. 2003, A&A, 400, 981 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Aubert, D., & Teyssier, R. 2008, MNRAS, 387, 295 [NASA ADS] [CrossRef] [Google Scholar]
 Aydi, E., & Mohamed, S. 2022, MNRAS, 513, 4405 [NASA ADS] [CrossRef] [Google Scholar]
 Bate, M. R., Bonnell, I. A., & Price, N. M. 1995, MNRAS, 277, 362 [Google Scholar]
 BermúdezBustamante, L. C., GarcíaSegura, G., Steffen, W., & Sabin, L. 2020, MNRAS, 493, 2606 [CrossRef] [Google Scholar]
 Bladh, S., Höfner, S., Aringer, B., & Eriksson, K. 2015, A&A, 575, A105 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Bladh, S., Liljegren, S., Höfner, S., Aringer, B., & Marigo, P. 2019, A&A, 626, A100 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Bohren, C. F., & Huffman, D. R. 1983, Absorption and Scattering of Light by Small Particles (New York: Wiley) [Google Scholar]
 Boulangier, J., Clementel, N., van Marle, A. J., Decin, L., & de Koter, A. 2018, MNRAS, 482, 5052 [Google Scholar]
 Bowen, G. H. 1988, ApJ, 329, 299 [NASA ADS] [CrossRef] [Google Scholar]
 Chan, T. K., Theuns, T., Bower, R., & Frenk, C. 2021, MNRAS, 505, 5784 [NASA ADS] [CrossRef] [Google Scholar]
 Chase, M. W. 1986, in NISTJANAF 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]
 Chen, Z., Blackman, E. G., Nordhaus, J., Frank, A., & CarrollNellenback, J. 2018, MNRAS, 473, 747 [NASA ADS] [CrossRef] [Google Scholar]
 Chen, Z., Ivanova, N., & CarrollNellenback, J. 2020, ApJ, 892, 110 [NASA ADS] [CrossRef] [Google Scholar]
 de ValBorro, M., Karovska, M., Sasselov, D. D., & Stone, J. M. 2017, MNRAS, 468, 3408 [NASA ADS] [CrossRef] [Google Scholar]
 Dorfi, E. A., & Höfner, S. 1991, A&A, 248, 105 [NASA ADS] [Google Scholar]
 Draine, B. T. 1981, ApJ, 245, 880 [NASA ADS] [CrossRef] [Google Scholar]
 Draine, B. T., & Salpeter, E. E. 1977, J. Chem. Phys., 67, 2230 [NASA ADS] [CrossRef] [Google Scholar]
 El Mellah, I., Bolte, J., Decin, L., Homan, W., & Keppens, R. 2020, A&A, 637, A91 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Ferland, G. J., Porter, R. L., van Hoof, P. A. M., et al. 2013, Rev. Mex. Astron. Astrofis., 49, 137 [Google Scholar]
 Ferrarotti, A. S., & Gail, H.P. 2002, A&A, 382, 256 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Fleischer, A. J., Gauger, A., & Sedlmayr, E. 1992, A&A, 266, 321 [NASA ADS] [Google Scholar]
 Forsythe, W. E. 2003, Smithsonian Physical Tables, 9th edn. (Norwich, New York: Knovel) [Google Scholar]
 Freytag, B., Liljegren, S., & Höfner, S. 2017, A&A, 600, A137 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Gail, H.P., & Sedlmayr, E. 1985, A&A, 148, 183 [Google Scholar]
 Gail, H. P., & Sedlmayr, E. 1987a, A&A, 171, 197 [NASA ADS] [Google Scholar]
 Gail, H.P., & Sedlmayr, E. 1987b, A&A, 177, 186 [NASA ADS] [Google Scholar]
 Gail, H.P., & Sedlmayr, E. 1988, A&A, 206, 153 [NASA ADS] [Google Scholar]
 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]
 Gail, H.P., Keller, R., & Sedlmayr, E. 1984, A&A, 133, 320 [NASA ADS] [Google Scholar]
 Gauger, A., Sedlmayr, E., & Gail, H.P. 1990, A&A, 235, 345 [Google Scholar]
 Gnedin, N. Y., & Abel, T. 2001, New A, 6, 437 [NASA ADS] [CrossRef] [Google Scholar]
 Grassi, T., Bovino, S., Schleicher, D. R. G., et al. 2014, MNRAS, 439, 2386 [Google Scholar]
 Hasegawa, K., & Umemura, M. 2010, MNRAS, 407, 2632 [NASA ADS] [CrossRef] [Google Scholar]
 Haworth, T. J., Ilee, J. D., Forgan, D. H., et al. 2016, PASA, 33, e053 [NASA ADS] [CrossRef] [Google Scholar]
 Höfner, S. 2008, A&A, 491, L1 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Höfner, S., & Dorfi, E. A. 1992, A&A, 265, 207 [Google Scholar]
 Höfner, S., & Freytag, B. 2019, A&A, 623, A158 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Höfner, S., Feuchtinger, M. U., & Dorfi, E. A. 1995, A&A, 297, 815 [Google Scholar]
 Höfner, S., GautschyLoidl, R., Aringer, B., & Jørgensen, U. G. 2003, A&A, 399, 589 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Homan, W., Montargès, M., Pimpanuwat, B., et al. 2020, A&A, 644, A61 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Homan, W., Pimpanuwat, B., Herpin, F., et al. 2021, A&A, 651, A82 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Hutchison, M., Price, D. J., & Laibe, G. 2018, MNRAS, 476, 2186 [Google Scholar]
 Iliev, I. T., Ciardi, B., Alvarez, M. A., et al. 2006, MNRAS, 371, 1057 [NASA ADS] [CrossRef] [Google Scholar]
 Iliev, I. T., Whalen, D., Mellema, G., et al. 2009, MNRAS, 400, 1283 [NASA ADS] [CrossRef] [Google Scholar]
 Jahanara, B., Mitsumoto, M., Oka, K., et al. 2005, A&A, 441, 589 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Jeong, K. S., Winters, J. M., Le Bertre, T., & Sedlmayr, E. 2003, A&A, 407, 191 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Jorissen, A. 2003, in Asymptotic giant branch stars, eds. H. J. Habing, & H. Olofsson (Berlin: Springer), Astron. Astrophys. Lib., 461 [Google Scholar]
 Jorissen, A., & Knapp, G. R. 1998, A&AS, 129, 363 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Kervella, P., Homan, W., Richards, A. M. S., et al. 2016, A&A, 596, A92 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 KesselDeynet, O., & Burkert, A. 2000, MNRAS, 315, 713 [NASA ADS] [CrossRef] [Google Scholar]
 Kim, H., & Taam, R. E. 2012, ApJ, 759, L22 [NASA ADS] [CrossRef] [Google Scholar]
 Krüger, D., & Sedlmayr, E. 1997, A&A, 321, 557 [NASA ADS] [Google Scholar]
 Krüger, D., Gauger, A., & Sedlmayr, E. 1994, A&A, 290, 573 [NASA ADS] [Google Scholar]
 Laibe, G., & Price, D. J. 2012, MNRAS, 420, 2345 [NASA ADS] [CrossRef] [Google Scholar]
 Laibe, G., & Price, D. J. 2014a, MNRAS, 440, 2147 [NASA ADS] [CrossRef] [Google Scholar]
 Laibe, G., & Price, D. J. 2014b, MNRAS, 444, 1940 [Google Scholar]
 Lamers, H. J. G. L. M., & Cassinelli, J. P. 1999, Introduction to Stellar Winds (Cambridge University Press), 452 [Google Scholar]
 Liu, Z.W., Stancliffe, R. J., Abate, C., & Matrozis, E. 2017, ApJ, 846, 117 [Google Scholar]
 Lodato, G., & Price, D. J. 2010, MNRAS, 405, 1212 [NASA ADS] [Google Scholar]
 MacLeod, M., & Loeb, A. 2020, ApJ, 902, 85 [NASA ADS] [CrossRef] [Google Scholar]
 Maercker, M., Mohamed, S., Vlemmings, W. H. T., et al. 2012, Nature, 490, 232 [NASA ADS] [CrossRef] [Google Scholar]
 Maes, S., Homan, W., Malfait, J., et al. 2021, A&A, 653, A25 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Malfait, J., Homan, W., Maes, S., et al. 2021, A&A, 652, A51 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Mastrodemos, N., & Morris, M. 1998, ApJ, 497, 303 [Google Scholar]
 Mattsson, L. 2016, Rev. Mex. Astron. Astrofis., 87, 249 [Google Scholar]
 Mattsson, L., & Höfner, S. 2011, A&A, 533, A42 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Mattsson, L., & Sandin, C. 2021, Universe, 7, 113 [NASA ADS] [CrossRef] [Google Scholar]
 Mattsson, L., Wahlin, R., & Höfner, S. 2010, A&A, 509, A14 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Mauron, N., & Huggins, P. J. 2006, A&A, 452, 257 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Mohamed, S., & Podsiadlowski, P. 2012, Balt. Astron., 21, 88 [Google Scholar]
 Neugebauer, M., & Snyder, C. W. 1962, Science, 138, 1095 [Google Scholar]
 Pawlik, A. H., & Schaye, J. 2008, MNRAS, 389, 651 [NASA ADS] [CrossRef] [Google Scholar]
 Petkova, M., & Springel, V. 2009, MNRAS, 396, 1383 [NASA ADS] [CrossRef] [Google Scholar]
 Petkova, M., & Springel, V. 2011, MNRAS, 415, 3731 [NASA ADS] [CrossRef] [Google Scholar]
 Petkova, M. A., Vandenbroucke, B., Bonnell, I. A., & Kruijssen, J. M. D. 2021, MNRAS, 507, 858 [NASA ADS] [CrossRef] [Google Scholar]
 Pinte, C., van der Plas, G., Ménard, F., et al. 2019, Nat. Astron., 3, 1109 [NASA ADS] [CrossRef] [Google Scholar]
 Preibisch, T., Ossenkopf, V., Yorke, H. W., & Henning, T. 1993, A&A, 279, 577 [NASA ADS] [Google Scholar]
 Price, D. J. 2012, J. Comput. Phys., 231, 759 [NASA ADS] [CrossRef] [Google Scholar]
 Price, D. J., & Federrath, C. 2010, MNRAS, 406, 1659 [NASA ADS] [Google Scholar]
 Price, D. J., & Laibe, G. 2015, MNRAS, 451, 813 [NASA ADS] [CrossRef] [Google Scholar]
 Price, D. J., Wurster, J., Tricco, T. S., et al. 2018, PASA, 35, e031 [Google Scholar]
 Rijkhorst, E. J., Plewa, T., Dubey, A., & Mellema, G. 2006, A&A, 452, 907 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Sandin, C. 2008, MNRAS, 385, 215 [NASA ADS] [CrossRef] [Google Scholar]
 Sandin, C., & Höfner, S. 2004, A&A, 413, 789 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Sandin, C., & Mattsson, L. 2020, MNRAS, 499, 1531 [NASA ADS] [CrossRef] [Google Scholar]
 Sharp, C. M., & Huebner, W. F. 1990, ApJS, 72, 417 [NASA ADS] [CrossRef] [Google Scholar]
 Skinner, M. A., & Ostriker, E. C. 2013, ApJS, 206, 21 [NASA ADS] [CrossRef] [Google Scholar]
 Struck, C., Smith, D. C., Willson, L. A., Turner, G., & Bowen, G. H. 2004, MNRAS, 353, 559 [NASA ADS] [CrossRef] [Google Scholar]
 Susa, H. 2006, PASJ, 58, 445 [NASA ADS] [Google Scholar]
 Tabak, R. G., Hirth, J. P., Meyrick, G., & Roark, T. P. 1975, ApJ, 196, 457 [NASA ADS] [CrossRef] [Google Scholar]
 Tegmark, M. 1996, ApJ, 470, L81 [NASA ADS] [CrossRef] [Google Scholar]
 Theuns, T., & Jorissen, A. 1993, MNRAS, 265, 946 [Google Scholar]
 Thorn, R. J., & Winslow, G. H. 1957, J. Chem. Phys., 26, 186 [NASA ADS] [CrossRef] [Google Scholar]
 Tsuji, T. 1973, A&A, 23, 411 [NASA ADS] [Google Scholar]
 van der Helm, E., Saladino, M. I., Portegies Zwart, S., & Pols, O. 2019, A&A, 625, A85 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Velli, M. 2001, Ap&SS, 277, 157 [NASA ADS] [CrossRef] [Google Scholar]
 Wang, N., & Lee, J. 2011, SIAM J. Sci. Comput., 33, 2536 [NASA ADS] [CrossRef] [Google Scholar]
 Whitehouse, S. C., Bate, M. R., & Monaghan, J. J. 2005, MNRAS, 364, 1367 [NASA ADS] [CrossRef] [Google Scholar]
 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 :
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
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 H_{2} and the second most abundant species is CO. As the temperature increases, H atoms become available and C_{2}H_{2} 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, C_{2}, C_{2}H, and C_{2}H_{2} that we consider for the evolution of carbon grains, the net growth rate (Eq. 17) writes
The first two terms of Eq. (17), which include P_{C} and P_{C2}, account for homogeneous grain growth and evaporation. In chemical equilibrium, the general form of the second term reads (Gauger et al. 1990),
where α(2) is an averaged sticking efficiency for a dimer, in this case for C_{2}, and the factor α_{*}(2) accounts for nonthermal equilibrium (nonTE) effects,
The quantity 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 and α(2). Similarly, we use α_{1} in the expression of . In chemical equilibrium, and when the gas has the same temperature as the dust, the correction factors b read b_{C} = b_{C2} = 1.
The third and fourth terms of Eq. (17) that include the contributions from the molecules C_{2}H and C_{2}H_{2}, 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
where the factor expresses nonTE effects,
α denotes the reaction efficiency of the growth reaction, and β the efficiency of the reverse reaction, being the TE value of the latter. We assume that both are equal (α^{c} = β^{c}) and set . Moreover, we set .
All Tables
Elemental abundances by number relative to hydrogen ϵ used in our models (Ferrarotti & Gail 2002).
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 (N_{p} = 10(2q − 1)^{2} + 2).
1D model parameters. T_{iso} corresponds to the temperature of the isothermal wind.
All Figures
Fig. 1. Nucleation rate per hydrogen atom () 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 
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 
Fig. 3. Planck mean values of the dust extinction efficiency Q_{ext} 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 
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 
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 T_{iso} = 1.25 × 10^{4} K. 

In the text 
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 
Fig. 7. Effect of the resolution setup on the predicted dynamics of an transonic adiabatic wind, for a fixed value of w_{ss}. 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 
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 
Fig. 9. Acceleration profile of a free wind for which the dust opacity is calculated using the nucleation formalism. 

In the text 
Fig. 10. Radial profiles of chemical and dustrelated 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 
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 (fulltext 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 4896 hours after online publication and is updated daily on week days.
Initial download of the metrics may take a while.