Three-component modelling of C-rich AGB star winds V. Effects of frequency-dependent radiative transfer including drift

Stellar winds of cool carbon stars enrich the interstellar medium with significant amounts of carbon and dust. We present a study of the influence of two-fluid flow on winds where we add descriptions of frequency-dependent radiative transfer. Our radiation hydrodynamic models in addition include stellar pulsations, grain growth and ablation, gas-to-dust drift using one mean grain size, dust extinction based on both the small particle limit and Mie scattering, and an accurate numerical scheme. We calculate models at high spatial resolution using 1024 gridpoints and solar metallicities at 319 frequencies, and we discern effects of drift by comparing drift models to non-drift models. Our results show differences of up to 1000 per cent in comparison to extant results. Mass-loss rates and wind velocities of drift models are typically, but not always, lower than in non-drift models. Differences are larger when Mie scattering is used instead of the small particle limit. Amongst other properties, the mass-loss rates of the gas and dust, dust-to-gas density ratio, and wind velocity show an exponential dependence on the dust-to-gas speed ratio. Yields of dust in the least massive winds increase by a factor four when drift is used. We find drift velocities in the range 10-67 km/s, which is drastically higher than in our earlier works that use grey radiative transfer. It is necessary to include an estimate of drift velocities to reproduce high yields of dust and low wind velocities.


INTRODUCTION
Winds of AGB stars are believed to be driven by radiation pressure on dust grains, which create an outflow when they in turn drag the atmospheric gas along. These winds are relatively slow (∼ 10 km s −1 ), but mass-loss rates can be high (∼ 10 −5 M ⊙ yr −1 ) owing to high densities. The type of dust forming in AGB-star atmospheres depends on the chemical composition of the gas: oxygen-rich stars (C/O < 1) form mostly silicate-type grains (but also iron dust can form in significant quantities, see Marini et al. 2019), whilst carbonrich stars (C/O > 1) form mainly amorphous carbon (amC) grains and smaller amounts of grains of SiC and polycyclic aromatic hydrocarbons (PAHs). The latter type of stars is usually referred to as "carbon stars" and represents evolved stars with initial masses in the range 1.5-4 M ⊙ that undergo so-called thermal pulses. That is, after the helium shell runs ⋆ Dedicated to Agda.
† E-mail: ChristerSandin@yahoo.se out of fuel, the star derives its energy from hydrogen burning in a thin shell; eventually, accumulated helium from the hydrogen burning ignites, causing a helium shell flash. During the thermal pulses, which last a few hundred years, material from the inner regions is mixed into the outer layers. This process is referred to as dredge-up and changes the surface composition of the star; in particular, this is how an oxygen-rich AGB star evolves into a carbon star. The amount of carbon expelled by carbon stars is significant and they may thus play role for the evolution of carbon (and carbonaceous dust) in the universe, although it cannot be ruled out that massive stars may be equally important (e.g., Gustafsson et al. 1999;Mattsson 2010). The carbon production of carbon stars is important and understanding the wind-formation mechanisms is essential to the full picture.
Radiatively accelerated dust grains exert a drag force on the gas; this drag force depends on how well grains couple to the gas, which in turn depends on the radiation flux, the density and temperature of the gas, as well as the extinction and cross section of dust grains. Grains that are massive enough (typically grains with radii > 0.1 µm) will have significant inertia and behave as so-called "inertial particles", which can decouple from their carrier fluid (the gas). Tiny grains, in comparison, can be regarded as coupled to a fluid element of its carrier fluid. In case gas and dust are perfectly coupled, often referred to as complete momentum coupling, the momentum gained by dust grains from the radiation pressure is immediately transferred (or "shared") with the gas. But if grains behave as inertial particles, only a fraction of the momentum gained from radiation is effectively driving the wind. As a result, dust and gas will have different flow velocities v and u, respectively, a phenomenon known as drift where v D = v − u. Not only is drift ignored when dust and gas are assumed to move with the same velocity [v D = 0, position coupling (PC)], but also the existence of inertia.
Most attempts at modelling dust-driven winds accounting for gas-to-dust drift use a stationary wind approach. Gilman (1972) checks the validity of complete momentum coupling in dust-driven winds through an analytical approximation, and finds that momentum coupling is complete for both O-rich and C-rich stars as the ratio of the drag force to the radiative pressure on dust grains ≫ 1/2. Kwok (1975) formulates a quantitative model of a dust-driven stellar wind, where he accounts for the drift velocity in grain growth, assumes complete momentum coupling, and he finds that grains are ablated when v D > 20 km s −1 . Goldreich & Scoville (1976) model mass loss and accounts for heating of the gas owing to drift. Berruyer & Frisch (1983) relax the assumption of complete momentum coupling and study the influence of different terms in three different regions; the authors find high drift velocities in their model where v D ≈ 100 km s −1 . MacGregor & Stencel (1992) scrutinize the momentum coupling completeness and find that this assumption is invalid for small grains (r d < ∼ 0.08 µm). Mastrodemos et al. (1996) use simplified models to justify the neglect of inertia of dust grains. Additional studies of dust-driven winds account for effects of drift in some respect (Salpeter 1974;Lucy 1976;Kwan & Hill 1977;Menietti & Fix 1978;Phillips 1979;Deguchi 1980;Tielens 1983;Papoular & Pégourié 1986;Dominik et al. 1989;Netzer & Elitzur 1993;Habing et al. 1994;Krüger et al. 1994;Ivezić & Elitzur 1995;Krüger & Sedlmayr 1997;Crosas & Menten 1997;Liberatore et al. 2001).
Amongst studies that consider drift, Elitzur & Ivezić (2001, hereafter EI01) and Ivezić & Elitzur (2010, hereafter IE10) deserve an honourable mention, as the authors prove the importance of drift analytically. Although they make several simplifying assumptions, their conclusions are robust: drift and reddening play crucial roles in shaping the velocity structure of dusty winds and the mass-loss rate must be strongly correlated with drift velocities.
Despite the demonstrated importance of drift, the number of time-dependent models is still sparse. Simis et al. (2001) calculate a dynamic wind model, which, however, does not include stellar pulsations or a detailed treatment of radiative transfer (RT); their model presents a too high density profile. The full system of radiation hydrodynamics with drift using one mean dust velocity is then, for the first time, modelled by Sandin & Höfner (2003a, hereafter Paper I), Sandin & Höfner (2003b, hereafter Paper II), and Sandin & Höfner (2004, hereafter Paper III) -these works are summarised with Sandin (2003) -who use grey RT and find that models including drift form more dust in the form of larger grains. The numerical approach of these models is improved with the work presented in Sandin (2008, hereafter Paper IV).
Somewhat surprisingly, gas-to-dust drift has in recent years been disregarded, despite the rich literature on stationary winds that include drift and the finding in Paper III that amounts of formed dust increase significantly when drift is included in the grain formation. However, the exclusion of drift is perhaps understandable in light of the limited effects observed on outflow properties in those pioneering works by Sandin & Höfner. On top of that, detailed models without drift appear to be in good agreement with observations (see, e.g., Mattsson et al. 2010, hereafter M10;Mattsson & Höfner 2011, hereafter MH11;Eriksson et al. 2014, hereafter E14;and Bladh et al. 2019, hereafter B19), which also ignore drift. The inclusion of drift is also a great computational challenge and the computational expense is another plausible reason why drift has been neglected. But from the point of view of understanding dust-driven winds, these justifications are unsatisfactory. Effects of inertia are undeniable and the timescale associated with drift -generally known as the stopping time -is strongly dependent on the gas-density gradient. To this date, drift must have a greater impact in models with frequency-dependent RT, knowing that such models have steeper density gradientsin particular in the regime of critical wind formation. There are therefore solid reasons for a scrutiny of the role of drift in more realistic models, in particular since the computational power required for this is less of an inhibiting factor than it was a mere decade ago.
With the simulation code we present here, T-800, we extend our gas-dust drift models with frequency-dependent gas and dust opacities, where we can choose between the small particle limit (SPL) and Mie scattering (following our approach in MH11), and we also calculate RT using a Feautrier-based solver (using the readily available description of Hubeny & Mihalas 2015, hereafter MH15) that allows the calculation of models with as many gridpoints as when grey opacities are used; this solver replaces the one of Balluch (1988), which is based on the approach of Yorke (1980). Amongst several numerical changes, we also use the high-accuracy advection scheme introduced in Paper IV. Our models are unprecedented in the respect that there are no similar stellar-wind models that include as much physics and can operate with the exceptional numerical accuracy needed to deal with drift. Our results vitiate the current consensus that drift plays a minor role (Höfner & Olofsson 2018, footnote 6).
Here, we re-evaluate the wind formation for a set of the model parameters in M10. Our study thereby also applies to E14. Our aim is to scrutinise basic differences between drift and PC models. We first describe the physics and numerical features of our models in Section 2. The modelling procedure and results are then presented in Section 3. We discuss the influence of drift on our results in Section 4 and summarise the paper with our conclusions in Section 5.

MODEL FEATURES AND IMPROVEMENTS
As in the four previous papers in this series, we consider three interacting physical components in the outer atmosphere and wind of the star: gas, dust, and radiation field. The three components are described by a system of equations that conserve and describe the interchange of mass, energy, and momentum. We call our simulation code T-800 and our simulation code for calculating initial structure John Connor (our analysis tool is similarly named Sarah Connor. Whilst the physical system is described in part in Paper I-Paper II and references therein, a number of adjustments merit a more complete description of the current capabilities. In comparison to other AGB star wind models, T-800 share most features with the darwin PEDDRO-type 1 models that originate with Höfner et al. (2003).
We describe the hydrodynamic equations and the physical terms of T-800 in the following three subsections: the gas component, Section 2.1; the dust component, Section 2.2; the radiation field, Section 2.3; and the numerical method, Section 2.4. All used abbreviations and symbols are collected in Tables A1-A3.

The gas component
Matter is present in either gas or dust phase. Five equations describe the gas phase: the equation of integrated mass, the equation of continuity, the equation of motion, the equation of inner energy, and the equation of number density of the condensible material: ∂ ∂t ρ g + ∇ · (ρ g u) = −m 1 S (2) ∂ ∂t (ρ g u) + ∇ · (ρ g u u) = −∇P g − Gm r r 2 ρ g + ∂ ∂t (ρ g e) + ∇ · (ρ g e u) = −P g ∇ · u + 4π ρ g (κ J J − κ S S g ) + where m r is the integrated mass at radius r, ρ g the gas density, t the time, u the gas velocity, m 1 the dust monomer mass, S the amount of condensible material in the gas phase (see below), P g the pressure, G the gravitational constant, c the light speed, κ H (κ J , κ S ) the gas opacity weighted with the first radiative moment H (the zeroth radiative moment J, the source function S g ), v the dust velocity, f drag the drag force, e the specific internal energy of the gas, ε the fraction of specular collisions between gas and dust particles, v D = v − u the drift velocity, and n C the number density of the condensible material. As before, an ideal gas law is used for the equation of state 1 Pulsation-enhanced dust-driven outflow (PEDDRO), see for example Höfner & Olofsson (2018) where T g is the gas temperature, γ = 5/3 is the ratio of specific heats, µ = 1.26 is the mean molecular weight, N A is Avogadros constant, and m p is the proton mass. Assuming local thermal equilibrium, the source function equals the Planck function, S g = B. Abundances of carbon and oxygen are specified as (log) number fractions relative to hydrogen, i.e. ǫ X = log 10 (n X /n H ) + 12, and the carbon-to-oxygen-ratio The number densities of the molecules in the gas phase that are part of the grain formation are calculated in an equilibrium chemistry of H, H 2 , C, C 2 , C 2 H, and C 2 H 2 . Partial pressures of single H and C atoms are calculated according to the description in (Gail & Sedlmayr 2014, chapter 10.3, hereafter GS14), and we use dissociation constants (K) of Sharp & Huebner (1990).

The dust component
Five equations describe the dust phase: four moment equations K 0 -K 3 of the grain-size distribution function (GS14), and the equation of motion of the dust.
where 0 ≤ j ≤ 3, ρ d = m 1 K 3 the dust density, τ −1 the net grain growth rate, N l is the lower size-limit of macroscopic grains (we use N l = 1000 carbon atoms), J ⋆ is the nucleation rate, and χ H the extinction coefficient weighted with the first radiative moment H. The net grain growth is gr is (homogeneous and heterogeneous) grain growth, τ −1 dc grain decay (by evaporation and chemical sputtering), and τ −1 ns ablation (by non-thermal sputtering), see Paper III for details. Dust particles are assumed to be spherical grains of amC. Effects of non-zero drift velocities are included in the grain growth and destruction terms according to the description in Paper III.
The moments of the size distribution allow calculation of average properties of the dust, including the total number density of dust grains n d = K 0 , the mean grain radius r d = r 0 K 1 /K 0 (r 0 is the monomer radius), the mean grain cross section σ = πr 2 0 K 2 /K 0 , 2 and the mean grain size N = K 3 /K 0 . The monomer radius is defined as (Gail et al. 1984, equation 2.2) where A is the atomic weight of the dust-forming species, and ρ m the intrinsic density of dust grains. The amount of condensible material S equals the right-hand side of equation (7) when j = 3, S = τ −1 gr K 2 + N l J ⋆ . The (average dust velocity) drag force is the same as in Paper I, and we only consider specular collisions (ε = 1): where C LA D is the limits approximation of the drag coefficient (Paper I, equation 23), T d the dust temperature (equation 23), and v ζ the thermal velocity v ζ = ζT g , and ζ = 128k B /(9π µm H ), where m H is the mass of a hydrogen atom.
In models that use position coupling (PC) instead of gas-dust drift, all terms in the dust equation of motion (equation 8) are added to the gas equation of motion (equation 3), and the equation is replaced with the relation v = u.

The radiation field component
The radiation field is described by the zeroth and first moment equations of the radiative transfer equation where J represents the radiative energy density, H the radiative energy flux, K is the second moment of the specific intensity, and q the sphericality (equation 16). The two moment equations depend on three moments of the radiation field, which is why it is necessary to solve the RT equation to calculate the Eddington factor f Edd , which gives K = f Edd J. We describe our approach to solve the RT equation next in Section 2.3.1. The interactions between the radiation field and the gas and dust are described separately in Sections 2.3.2 and 2.3.3. Regarding the frequency-dependent problem of RT, we follow the guidelines of Mihalas & Weibel-Mihalas (1984, §82).

Radiative transfer
We solve the equation of RT in spherical geometry, without scattering. Up to now, all our models have used the timeindependent spherical-geometry method of Yorke (1980), Balluch (1988), andBodenheimer et al. (2007); the method considers rays that are individually first integrated inwards and then outwards. Whilst this approach has worked well in models using grey RT and reasonably well in models using frequency-dependent RT, we found that it sometimes is unstable.
Following the literature (Mihalas & Weibel-Mihalas 1984;Annamaneni 2002, MH15), we chose to write a new solver based on a differential-equation technique that uses a Feautrier-type solution along individual impact parameters. We used the description of (MH15, chapter 19.1); because of numerical inaccuracies in optical shells that are more thin, we found it is necessary to either use quadruple precision when solving the resulting tridiagonal system of equations or rewrite the equations according to the mention in Nordlund (1982, section 3.9.1) and description in Rybicki & Hummer (1991, appendix A).
The inputs to the solver are the total absorption coefficient κ ν,g+d = κ ν + χ ν /ρ g , where κ ν and χ ν are the frequencydependent gas opacity and the dust extinction efficiency, respectively; the source function where B ν is the Planck function; the radiative temperature at the outer boundary, T ext r = 0 K; and a frequencydependent expression for the radiative flux at the inner boundary H int , which is placed where the diffusion limit applies (MH15, equation 11.176), where κ R is the Rosseland mean opacity (equation 17). The RT calculations yield mean-intensity-like and flux-like variables that are integrated over impact parameters according to the description in Yorke (1980) to yield the three radiative moments J ν , H ν , and K ν . And these moments are in turn used in the equations in the next two sections to calculate frequency-integrated properties.
The Eddington factor f Edd and the sphericality factor q are calculated as and where r C is the radius at the inner boundary. The frequencyintegrated zeroth and second radiative moments are overlined here to indicate that these properties are calculated from the output of the RT calculations and not from the hydrodynamic equations. The computing time using N ν frequencies, N D gridpoins, and N C core rays scales as t The radiative-transfer calculations can with advantage be executed in parallel for individual frequencies since no scattering is used.

The interaction between radiation field and gas
For each wavelength, the radiative-transfer equation is solved for the previous time step using the radial structure of the gas density, opacity, and temperature. The first three radiative moments J ν , H ν , and K ν that result from the calculations are integrated over frequency to find the following moment-weighted gas opacities where κ S is the Planck mean opacity. Gas opacities are provided as tabulated values κ ν (ρ g , T g , ν). The tables are created with the coma code (Aringer 2000;Aringer et al. 2009) and include updates regarding abundances (solar composition), frequency interpolation and resolution (B. Aringer, priv.comm.). The frequency integral (0 ≤ ν < ∞) is simply taken as the range of frequencies that is available in the tabulated data.

The interaction betweeen radiation field and dust
The extinction coefficients of the dust that correspond to the weighted gas opacities are calculated as follows using the radial structure of the dust extinction coefficient, the dust temperature, and the mean grain radius ( r d ): where, B d,ν is the Planck function at the dust temperature T d , κ d,ν = χ ν /ρ g , and we also have that (MH11) where a gr is the grain radius, a gr ≡ r d the mean grain radius, and Q ′ abs,ν = Q abs,ν /a gr the absorption efficiency. The absorption efficiency Q abs,ν used to calculate χ J and κ S is and the absorption efficiency Q abs,ν (pr) used to calculate χ H accounting for radiation pressure is where Q ext,ν is the extinction efficiency, Q sca,ν the scattering efficiency, and cos θ ν the average scattering angle. Finally, the Rosseland dust extinction χ R is calculated assuming Q sca,ν = 0, which is why Q abs,ν ( χ R ) = Q ext,ν . The frequency integral is also here taken as the range of frequencies that is available in the table.
Using Mie scattering, we followed the approach of MH11. Whilst calculating Q ext,ν (T d , r d ), Q sca,ν (T d , r d ), and cos θ ν (T d , r d ), we did not as in MH11 use bhmie of Bohren & Huffman (1983) 3 , instead we implemented the theory as described by GS14 (chapter 7.3). The calculation of these extinction coefficients requires tabulated values of the refractive indices n ν (phase velocity) and k ν (extinction coefficient). Assuming SPL, we use (m ν = n ν + ik ν ; e.g. Wickramasinghe 1972, equation 2.30) 3 The code bhmie modified by B. Draine and others is available at https://www.astro.princeton.edu/∼draine/scattering.html. and assume that the average scattering angle vanishes, cos θ = 0. We use the approach of GS14 (chapter 8.3; also see H03) to calculate a unique dust temperature as where σ SB is the Stefan-Boltzmann constant and the radiative temperature T r = (πJ/σ SB ) 1 4 . We also compare our new results with grey calculations where we have assumed that T d,grey = T r and that (see equation 4 in Paper I)

Numerical method
The eleven partial differential equations -equations (2)- (8) and (11)-(12) -are discretised in the volume-integrated conservation form on a staggered mesh, which together with the integrated-mass equation (equation 1) and the adaptive grid equation (see, e.g. Dorfi & Höfner 1991, section 3.2) form a system of thirteen equations. The non-linear system of equations is solved implicitly using a Newton-Raphson algorithm where the Jacobian of the system is inverted by the Henyey method. The flow of mass, energy, and momentum between the gridpoints is described with a second-order volume weighted van Leer advection scheme (Dorfi et al. 2006); details of the implemented advection scheme are described in Paper IV. We give further details of the basics of the numerical method used with T-800 in, for example, Paper I and references therein.

Using the adaptive grid equation or a mostly fixed grid
The fundamental problem with using a single adaptive grid equation for two fluids (gas and dust) or three components (gas, dust, and radiation field) is that it is very difficult to trace features in both the gas and the dust at the same time using a limited number of gridpoints N D . Modelling tests also reveal difficulties when the frequency-dependent carbon-related opacities we have are used together with N D > ∼ 100.
In comparison to our earlier work, we do not use the adaptive-grid equation to resolve shocks or other features. Instead, we fix the outer parts of the grid and set N D = 1024. In this approach, we do not let the outer boundary of the grid expand from the value used in the initial-model calculations to the pre-determined outer boundary of the wind calculations (as we do in Paper I-Paper IV and is done in all other implementations of stellar winds using the adaptive grid equation); the radial range of the initial model is instead set to the full range of the wind right away where the outer boundary typically is set to 40R ⋆ , but with slower winds we use 20R ⋆ . We calculate a grid using a logarithmic distribution of gridpoints. In some models, we find that it is necessary to increase the spatial resolution in the centre parts to achieve convergence, which is why we doubled the number of gridpoints inside of r 950 of gridpoint 950 and redistribute the remaining 950 gridpoints in the region outside of r 950 .
Owing to a grid that is fixed where there is dust, we use the same amount of artificial viscosity in both PC and drift models; according to the description in Paper I (equations 13 and 14), we set the length scale to l av = 3.5 × 10 −3 r.
We also do not use any artificial mass diffusion. (Notably, in Paper III and Paper IV, we use both artificial mass diffusion and a length scale twice as high in the presented models, where gridpoints move about to some extent throughout the model domain.) With one model setup, we instead tested using l av = r, with negligible differences in resulting physical structures.
We compare some of the new models with models that use the adaptive grid equation, where the gas density and energy are resolved; these models use N D = 100 and extend out to 25R ⋆ (see Section 4.1).

Modelling stellar pulsations at the inner boundary
We model effects of stellar pulsations on the atmosphere and wind region using a piston boundary condition, which is a sinusoidal and radially varying inner boundary that is placed above the region where the, so-called, κ-mechanism supposedly originates. The piston boundary is described with the period P and the amplitude ∆ u p . The adaptive grid equation is used to allow gridpoints where r < 2R ⋆ to stretch with the piston. We do not model any inflow of mass through the inner boundary, as, for example, Simis et al. (2001, section 4.2) and Woitke (2006, section 2.5.1) do. The typical amount of mass in the modelled envelope is 0.14-2.3% of the total mass including the core (see Table 1) before the dynamical modelling of the stellar wind begins. As the model domain is rather quickly depleted of material owing to the stellar wind, long-term modelling, covering thousands of years, is problematic in the current approach.

Additional considerations in the new models
Rosseland and Planck mean opacities and extinction coefficients are calculated when the tabularised frequencydependent opacity and extinction data are loaded; the mean extinction coefficients are calculated for a set of pre-defined temperatures and grain radii. The equation of RT is solved for the previous time step in the first iteration, before the non-linear system of equations is solved for the current iteration. T-800 executes the RT calculations in parallel for individual frequencies using a hybrid approach that makes use of both OpenMP and MPI. 4 The frequency-integrated weighted opacities and extinction coefficients are divided with the Rosseland and Planck mean opacities and extinction coefficients at the first iteration to calculate opacity ratios k it=1 X as follows Only the Rosseland and Planck mean opacities and extinction coefficients are calculated in subsequent iterations, and those values are then multiplied with the opacity ratios k it=1 X of the first iteration to get the respective value in the current iteration.

MODELLING PROCEDURE AND RESULTS
We first describe our modelling procedure in Section 3.1, and the physics setup and choice of model parameter sets in Section 3.2. We present our results in Section 3.3.

Modelling procedure
Our modelling procedure consists of four separate stages. We begin by calculating a hydrostatic dust-free initial model using John Connor (Appendix B). The initial model spans the radial range r int , r ext ≈ [0.9R ⋆ , 1.8R ⋆ ]. The inner radius r int is set as small as possible, but the model will not converge if it is lower than a model-dependent threshold radius. Simultaneously, the external radius r ext cannot be too large in the steep hydrostatic structure as that also prevents the model from converging. We have found empirically that the external radius is well selected such that ρ g (r ext ) = 10 −6 ρ g (r int ).
The model domain of the converged initial model is extended to instead use the outer radius r ext final = 40R ⋆ before it is saved to a file. Models with an expected low terminal velocity (u ∞ < ∼ 10 km s −1 ) instead use r ext final = 20R ⋆ to save calculation time. Physical properties are not used in the radial range r ext < r < r ext final , which is used as a "gridpoint reservoir" in the expansion stage.
In the relaxation stage, the initial model of John Connor is relaxed using the system of equations of T-800, still without simulating stellar pulsations or allowing dust to form. A model has relaxed to be hydrostatic (using the dynamic code) when the time step becomes higher than 10 14 s. In the expansion stage, the piston is switched on from zero to full amplitude in 2 pulsation periods. Gridpoints i where r i < 2R ⋆ thereby move along with the inner boundary, whilst remaining gridpoints are held fixed. Simultaneously with the piston activation, all dust equations and terms are switched on. Dust begins to form in the cooler outer layers of the initial model, i.e. at r ≈ 2R ⋆ ; dust grains absorb radiative momentum and are accelerated outwards whereby they drag the gas along, thus initiating a wind. The expansion of the physical region proceeds until the outer region has moved from r ext and reached r ext final of the fixed grid. In the final wind stage, the wind model is evolved for a time interval of about 12-200 P. The modelled time interval depends on the time that is needed for transients in the expansion phase to leave the model domain through the outer boundary. Temporally averaged properties are measured only after this time. The interval is longer with models that do not develop a periodic or a close to periodic structure.

Physics setup and selection of model parameters
Our approach is to use a physical setup that mostly is identical to that of M10. We used 319 frequencies in the RT calculations, whilst M10 use 64 frequencies. The abundances are set to solar values, with the exception of the carbonto-oxygen ratio that is an input parameter. We assume ∆ u p = 4 km s −1 and use M ⋆ = 1.0 M ⊙ . We use the same pulsation periods as before, but note that these are somewhat different from the period relation of Wood (1990, see equation 2.44 in Lamers & Cassinelli 1999, where the form below was achieved using equation B5 and assuming T eff /T eff,⊙ = 0.572 and M ⋆ = 1.0M ⊙ ) and the period-luminosity (P-L ⋆ ) relations of Feast et al. (1989) and Groenewegen et al. (1998) m bol (Feast 1989 m bol (Wood 1990) * = −2.58 log P + 20.33, and (28) m bol (Groenewegen 1998) = −2.59 log P + 20.52.
Here, we used the same distance modulus to the Large Magellanic Cloud as Groenewegen et al. (1998), µ D = 18.50, and the bolometric magnitude of the sun M bol,⊙ = 4.74. The bolometric magnitude m bol is converted to a luminosity using 2.5 log(L ⋆ /L ⊙ ) = M bol,⊙ + µ D − m bol . The resulting P-L ⋆ -relations are shown in Fig. 1 along with the values we use. (For reference, the figure also shows the P-L ⋆ -relation for O-rich Miras of Whitelock et al. 2009.) The values we use match the relation of Wood (1990) the best. The dust consists of amC, where the intrinsic dust density ρ m = 1.85 g cm −3 , the surface tension σ grain = 1400 erg cm −2 , and the tabulated values (n ν and k ν ) used to calculate the extinction are taken from Rouleau & Martin (1991). We assume SPL, but calculate three models using Mie scattering as comparison. Furthermore, amongst other parameters, Andersen et al. (2003) study the role of the sticking coefficients to properties of stellar winds that use frequency-dependent RT and find that mass-loss rates, outflow velocities, and the degree of condensation change. More recent models (beginning with Mattsson et al. 2007) use sticking coefficients that are set to unity as these form a stellar wind easier owing to larger grains than when the original lower values are used. We use the original values that are used in all models in this series of articles Ξ 0.34 Table 1. Model parameters, see Section 3.2 for further details. The model name is given in Column 1. The following five columns specify: stellar luminosity L ⋆ , effective temperature T eff , chemistry (C/O ratio), and pulsation period P. The last column shows the fraction of the stellar mass contained in the radial domain of the initial model. The stellar mass M ⋆ is nearly 1.0 M ⊙ in all models, and the pulsation amplitude ∆ u p = 4 km s −1 . parison of results of three models that use both Ξ 0.34 and Ξ 1.00 (ξ C = 1.00, ξ C 2 = 1.00, ξ C 2 H = 1.00, ξ C 2 H 2 = 1.00) shows non-trivial differences (Appendix C). The model calculations are computationally demanding, and we do not calculate a model grid at this point. The model setups are taken from M10 (table 2). We show our complete set of model parameters in Table 1. We collect all physical and numerical assumptions in Table D1 for easy reference.

Results
We selected the models in this paper with the intention of studying relative changes compared to previously calculated wind models. This approach allows a quantitative and qualitative estimate of the importance of the adopted differences in modelling and physics.
Our main objective here is to introduce effects of drift, using one mean dust velocity. The drift models are compared to (non-drift) PC models that in all other aspects are the same as the drift models.
In addition to a modelled drift velocity v D in drift mod-els, we calculate the equilibrium drift velocityv D for both PC and drift models assuming complete momentum coupling. That is, we assume that the drag force is equal to the radiative pressure on dust grains after subtracting the gravitational pull on the same dust grains (this term is negligible, but we keep it for completeness of the argument). Assuming specular collisions, equations (8) and (10) give If we ignore the gravitational term as well as the thermal velocity term v ζ and use M = 4πr 2 ρ g u, we get a simpler expression valid at supersonic velocities Notably, equation 31 can also be applied to results of PC models, but the resulting values ofv D are calculated with models that disregard dilution of the dust and are therefore not comparable with v D of drift models. Wind models are characterised with properties temporally averaged at the outer boundary. Two properties characterize the gas: the gas mass-loss rate M and the gas terminal velocity u ∞ . The dust is characterized with four properties: the degree of condensation f cond , the dust-togas mass-loss ratio M d / M , the mean grain radius r d , and the terminal drift velocity v D,∞ (only for drift models). The dust-to-mass mass loss ratio is where Moreover, since the dust component is diluted by the drift factor, we define the degree of condensation as where ρ tot c is the total density of condensable matter (present in both the gas and the dust phases). Assuming PC, the expression reduces to Providing a measure for the variability of the model structure, each outflow property (Q) is accompanied by a relative fluctuation amplituder = σ s /Q, where σ s is the (sample) standard deviation of the property Q in the measured time interval.
We show results of all our model calculations using Ξ 0.34 in Table 2. The average mass loss properties of the drift models are plotted against F D in Fig. 2. Moreover, the average dust mass loss rate is plotted against the drift velocity in Fig. 3. The four remaining averaged properties -outflow velocity, degree of condensation, grain radius, and drift velocity -are plotted against F D in Fig. 4. Finally, the average mass loss ratio is plotted against the average grain radius in Fig. 5.

DISCUSSION
We discuss the following topics in the next six subsections: the role of the spatial resolution (Section 4.1), wind formation when models include drift (4.2), temporal variability and averaged properties (4.3), the importance of using Mie scattering in place of SPL (4.4), a comparison with observations (4.5), and a comparison with the theory of EI01 (4.6). Table 2. Temporally averaged quantities at the outer boundary; see Section 3.3. From the left, the first two columns specify the model name (see Table 1) and if PC or drift is used (P/D). Six column pairs show the averaged: mass loss rate M , terminal velocity u ∞ , degree of condensation f cond , dust-to-gas mass-loss rate δ dg F D , dust radius r d , and drift velocity v D,∞ (only for drift models). A relative fluctuation amplituder is prvided for each property; a subscript m (c, w, µ) indicates that the shown value was multiplied with a factor 10 3 (10 2 , 10 4 , 10 6 ). The final columns show the outflow classification class: irregular (i), periodic (l×p), and quasi-periodic (lq); l indicates the (multi-)periodicity of the gas/dust outflow in the unit of the piston period P, and if the dust extinction is calculated using Mie scattering instead of the default SPL. Rows of drift models are shown in boldface.

How results depend on the used spatial resolution
A strong argument in favour of using an adaptive grid equation is that such models are able to resolve shocks; this is shown by, for example, Dorfi & Feuchtinger (1991) and Feuchtinger & Dorfi (1994) for single shocks and pulsations. We show in Paper IV (where N D = 500 and 700) that winds are better modelled without resolving shocks. The adaptive grid equation is still used to move the inner boundary to simulate stellar pulsations. Models become smoother and also often periodic than when the adaptive grid equation is used to resolve shocks, contrary to what is claimed in the literature. In this paper, we took the approach of Paper IV further and kept the grid fixed at all radii r > 2R ⋆ , which alleviates the numerical advection of dust in drift models; the spikes seen in plots of the drift velocity in Paper I and Paper II are thereby avoided as the relative velocity between grid points and the dust is always greater than zero (Appendix E). We also increased the number of gridpoints to N D = 1024, which is more than a factor 10 higher than what is currently used with extant results of the darwin stellar wind model (e.g., H03, M10, E14, B19), who use N D = 100. We also used the higher-accuracy advection scheme of Paper IV. All shocks of the gas and discontinuities of the dust are resolved with at least one gridpoint across a shock front or a discontinuity with this approach; we illustrate this in Fig. 6 where we show the physical structure of the PC model L3.85T28E88. And they are resolved at all times as the gridpoints are fixed. Shocks are less steep in the outer parts owing to the artificial viscosity, which length scale l av ∝ r, see equation (25). A large fraction of the PC models (mainly) are periodic, mostly with a periodicity 1P (see Section 4.3.1).
In comparison, models that use the adaptive-grid equation and N D = 100 track one or two shocks in the outer parts of the wind, and show much higher variability; such models appear to be unable to adequately resolve multiple shocks. The reason is that most gridpoints gather about one or two shocks leaving remaining regions unresolved, see Fig. 6. Thereby, it is hardly possible to follow shocks that develop with each pulsation period; there are no gridpoints available to resolve them. Our conclusion could perhaps have been different if N D was higher in these models. Notably, our comparison of benchmark test results shows reasonable agreement between results of T-800 and darwin that use the same physics setup (Appendix F). Our results suggest that wind models are more accurate overall when the adaptive grid equation is not used and all shocks and all regions are resolved. This is particularly important in drift models where dust is not affixed to gas shocks.

Wind formation in our new models
4.2.1 The radial structure of model L3.70T28E88 We illustrate the physical structure for the PC and drift models of setup L3.70T28E88 in Fig. 7; this drift model shows the smallest dust mass-loss rate of all our models. (We show a similar figure for the PC and drift models L3.85T28E88 using Mie scattering instead of the SPL in Fig. 12). The PC model shows a periodic structure with small-amplitude variations in both the gas and the dust. Dust formation takes place throughout the envelope (Fig. 7e), although at much reduced rates in the outer parts and the rate of grain growth at, say, radii r > ∼ 10 R ⋆ is usually negligible. Neither the degree of condensation (Fig. 7d) nor the grain radius (Fig. 7f) increase by much for r > ∼ 7R ⋆ . The degree of condensation is rather low compared to denser grey opacity models (cf. figs. 2 and 4 in Paper IV). We calculated the equilibrium drift velocityv D (equation 30), which reaches about 20 km s −1 already at small radii (Fig. 7c); this velocity is supersonic for r The drift model also shows a periodic structure. The mass loss rate is less than half of the corresponding PC model. The drift velocity (Fig. 7c) attains values of v D = 20 km s −1 already at r ≃ 2.5R ⋆ , and increases to about 30-35 km s −1 at larger radii; this is in sharp contrast to earlier grey models where drift velocities were always small (v D < 16 km s −1 in Paper I-Paper IV; cf. figs. 2c and 4b in Paper IV). The equilibrium drift velocityv D of the drift model is very similar to v D for r < ∼ 10 R ⋆ and shows larger deviations of up to 5 km s −1 in less dense regions at larger radii. Our own tests show that the momentum coupling is nearly complete (Appendix G). There is -in this casesome similarity betweenv D of the drift model and that of the PC model (see Section 4.2.2). The drift velocity is supersonic at large Mach numbers throughout the radial domain The variations in the drift velocity and the dust density ( Fig. 7b) show discontinuities that are moving outwards at speeds more than twice as high as the gas. The dust appears to accumulate in separated shells, where there is little dust between the shells. The high variability of the drift model makes a direct comparison of amounts of dust between the PC and drift models difficult (Fig. 7d).
In the shown snapshot, dust formation mostly occurs in a narrow region, where 1.5 < ∼ r < ∼ 2.2R ⋆ (Fig. 7e). At larger radii, gas particles do not stick to dust particles any more if drift velocities are too high (cf., e.g., fig. 1 in Paper IV). In such locations, grain decay dominates, and in particular grains are ablated by gas particles (τ −1 ns ). However, the grain decay is inefficient.
It is easier to scrutinise the wind formation in a temporal average of the radial structure in models of pulsating atmospheres. We show such an average plot for the same model in Fig. 8. The figure reveals smooth structures in all shown properties of both models. In particular, the gas velocity structure is -in this case -very similar to the selfsimilar structure presented by IE10 (equation 1, at small optical depth) where we set the dust condensation radius r c = 2R ⋆ . The drift model shows structures that are as smooth as   Figure 7. Radial structure of a snapshot of setup L3.70T28E88 for the full modelled region, using the sticking-coefficients setup Ξ 0.34 . The drift (PC) model is shown with purple (orange) lines. From the top left, the 12 panels show: (a) gas velocity u, sound speed c s , and self-similar gas velocity u ss ; (b) gas density ρ g , dust density 10 4 × ρ d (log); (c) drift velocity v D , equilibrium drift velocityv D , and thermal velocity v ζ ; (d) degree of condensation f cond ; (e) net growth rate τ −1 gr , net decay rate τ −1 dc , nucleation rate J ⋆ , and non-thermal sputtering rate τ −1 ns (log); (f) average grain radius r d ; (g) Eddington factor f Edd ; (h) extinction coefficient χ H , grey extinction coefficient χ; (i) gas temperature T g , radiative temperature T r , and dust temperature T d ; (j) opacity κ H and Rosseland mean opacity κ R (log); (k) temperature ratios T d /T r and (T g /T r ) eq ; and (l) extinction coefficient ratio χ H /χ. All properties are drawn as function of the stellar radius R ⋆ (lower axis) and astronomical units (AU; upper axis). Grey horizontal lines are guides.
those of the PC model. The mass loss rates of both the gas and dust are constant for r > ∼ 5R ⋆ (Figs. 8d and 8f). Meanwhile, both the degree of condensation (Fig. 8c) and average grain radius (Fig. 8e) show a significant increase throughout the model domain; this cannot be explained by dust formation as both grain growth and decay are negligible (Fig. 8g). Instead, the slope appears when dust leaks (cf. "free streams", Hopkins & Lee 2016) into more dust free regions between dust fronts (such dust fronts are seen in Figs. 7b and 7d); more dust has moved into regions between fronts as the wind reaches larger radii. The larger the separation between dust fronts, the steeper the slope. Finally, also the drift velocity increases with the radius (Fig. 8h).

Radial structure properties of the full model sample
In our analysis of radial structures of the full model sample, we consider four properties that show some change: outflow velocity, drift velocity, dust formation, and amount of formed dust. We first discuss our PC model results and thereafter our drift models.
Outflow velocities of PC models typically increase by 10-70 per cent in the radial interval 10-40 R ⋆ . Values are usually higher with higher model luminosity (L ⋆ ) and carbon-to-oxygen ratio (C/O). The increase is lower in the radial interval 20-40 R ⋆ , up to 20 per cent. The final outflow velocity is reached already at 20R ⋆ in model L3.85T30E88. The self-similar solution for the outflow velocity u ss (equation 35) provides a good description of the velocity structure  Figure 8. Temporally averaged radial structure of setup L3.70T28E88 for the full modelled region. The drift (PC) model is shown with purple (orange) lines. From the top left, the eight panels show the: (a) gas velocity u, sound speed c s , and self-similar velocity u ss ; (b) gas density ρ g , dust density 10 4 × ρ d (log); (c) degree of condensation f cond ; (d) dust mass loss rate M d (log); (e) average grain radius r d ; (f) dust-to-mass loss rate ratio M d / M, dust-to-gas density ratio δ dg (log); (g) net growth rate τ −1 gr , net decay rate τ −1 dc , nucleation rate J ⋆ , and total non-thermal sputtering rate τ −1 ns (log); and (h) drift velocity v D . All properties are drawn versus the stellar radius R ⋆ (lower axis) and astronomical units (AU; upper axis). Grey horizontal lines are guides.
in five models: L3.70T26E85, L3.70T28E88, L3.85T24E85, L3.85T28E85, and L3.85T28E85, where all but one model have a low C/O ratio; the velocity structure is somewhat to much less steep than u ss in the remaining models.
Grain growth occurs at some rate throughout the model domain, but is balanced by grain decay through evaporation and chemical sputtering for r < ∼ 2R ⋆ . In our examination of the radial dust mass loss rate structure, we see that dust formation is complete (within a few per cent) at r ≈ 5-10 R ⋆ , typically, but in a few models it appears that a larger radial interval is needed, r ≈ 20-30 R ⋆ . In models with a higher mass loss rate, the dust mass loss rate structure increases more or less monotonically with the radius, but in all models with a lower mass loss rate, a peak is reached at about r = 2 R ⋆ , where grain decay causes decreased values out to, say, r ≈ 3 R ⋆ (L3.8T24E85, L3.70T26E85, L3.85T26E85, L4.00T26E82, L3.85T28E85, and L4.00T32E88).
In drift models, the terminal velocity is typically reached at shorter radii. The self-similar solution describes average outflow velocities well in three models: L3.70T24E88, L3.70T28E88, and L3.85T28E88; the velocity structure is somewhat to much less steep than u ss in the remaining models. There is one exception, the velocity structure of M3.85T30E88 increases more than u ss ; also here, the optically thin exponent of u ss (2/3) provides a better fit than when using the optically thick exponent (2/5, see IE10). The outflow velocity increases by up to 40 per cent in the radial interval 10-40 R ⋆ (up to 12 per cent for 20-40 R ⋆ ).
The drift velocity varies with the model setup and the radius. Dust grains accelerate fast from, say v D ≃ 2-10 km s −1 when they form at r < ∼ 2 R ⋆ , to higher velocities v D ≃ 20-100 km s −1 at 3R ⋆ < ∼ r < ∼ 5-10 R ⋆ . Values are typically higher in regions between more dense shells of dust where 5 < ∼ v D < ∼ 30 km s −1 . Such differences are not seen in a temporal average of the radial structure. The equilibrium drift velocity overlaps the drift velocity well at lower radii r < ∼ 10 km s −1 , and mostly shows a somewhat larger deviation at larger radii where the difference is 5-50 km s −1 . Differences are larger in regions between dust shells. The equilibrium drift velocity of the PC models is sometimes similar to the drift velocity throughout the model domain (L3.85T2400E91 and L3.70T3000E91), but is more often similar in denser dust shells. For most models, the equilibrium drift velocity bears little resemblence to the actual drift velocity, which is not strange considering that the two values are calculated using different physical structures! Grain formation would be interrupted by ablation when drift velocities are high enough. Typically, such higher values are reached where r > ∼ 4R ⋆ , and then in regions where there is less dust between dust shells. We find that ablation has a minor effect on the dust formation and is unimportant in models where v D,∞ < ∼ 30 km s −1 . Our scrutiny shows that it is necessary to model a larger region that extends out to, say, 40 R ⋆ to calculate a more accurate outflow velocity. The full region of dust formation should mostly be covered in a model that extends to 10R ⋆ . Currently, our models use one average dust velocity for grains of all sizes. It is possible that our results would be different if the models would include grain size-dependent dust velocities, which could affect the drift velocity and thereby the rate of ablation for grains of different size.

Reasons for higher drift velocities than in grey models
Our results using grey and constant opacities show minor effects of drift (Paper I-Paper IV), where the average drift velocity is typically v D ≈ 5 km s −1 ; for the full sample of grey models, v D < ∼ 16 km s −1 . The results we present here show larger effects. In our analysis of the reason behind this discrepancy, we here examine the terms that are different in the two approaches.
The dust temperature is calculated differently in the grey and frequency-dependent approach. In grey models, we have used T d,grey = T r , whilst in frequency-dependent models, the radiative temperature is weighted with the extinction coefficient ratio ( χ J / χ S ) 1/4 (equation 23). Furthermore, in radiative equilibrium, Temperature ratios T g /T r eq that deviate from 1 indicate an importance of non-grey RT (cf. section 3.1 in H03). We show the gas, dust, and radiative temperatures as well as the temperature ratios T d /T r and T g /T r eq in Figs. 7i and 7k (cf. figs. 3c, 5b, and 5c in H03). In this context, our PC model ratio T d /T r shows a good agreement with what H03 present for a model with different parameters. Meanwhile, our model shows a ratio that is lower than theirs; the ratio decreases towards T g /T r eq ≃ 0.5 for larger radii whilst their value is T g /T r eq > 0.8 for r > ∼ 4R ⋆ . The drift model ratios are similar, with a somewhat steeper temperature gradient of the gas, which might indicate an even stronger importance of non-grey RT when drift is included.
A fundamental difference between grey and frequencydependent models is how the extinction is calculated. The ratio between χ H (equation 18) and χ (equation 24) for L3.70T28E88 is about 1.8 when dust first forms and increases to larger values with the radius, see Fig. 7l (cf. fig. 5d in H03). The ratio increases somewhat faster with radius in the drift model.
To examine consequences owing to this discrepancy more closely, we calculated a drift model of setup L3.70T28E88 where we used the grey extinction χ instead of the frequency-dependent weighted average extinction χ H . We show the resulting radial structure of the two drift models (that use χ H and χ) in Fig. 9; the figure shows all variables in the equilibrium drift velocityv D (equation 31) except the radiative flux H, which is nearly identical in the selected snapshots of the two models.
The drift velocity v D is significantly lower in the model that uses the grey extinction χ, 10 < ∼ Fig. 9a. Dust grains collect in shells with nearly dust free regions between shells (Figs. 9c and 9d); the shells are smeared out when moving outwards. The gas density ρ g is similar throughout the radial domain, Fig. 9b. The dust extinction χ (the dust moment K 2 ) is, moreover, a factor 10 8 (10 7 -10 9 ) higher at the front of the dust shells than between them; the ratio is smaller towards the outer boundary. The ratio between the dust extinction and the dust moment K 2 (Fig. 9e) illustrates more clearly that these two properties are responsible for the higher drift velocity when using χ H instead of χ (see equation 30). However, because the physical structure changes with the different physical conditions (whence all three variables are modulated), we have not proven that the lower grey extinction is the only reason behind the differences.

Parameter combinations that fail to form a wind
Models fail to form a wind when too little dust forms. Such conditions are characterized by smaller amounts of carbon C/O and low luminosity L ⋆ , as well as high effective temperature T eff . In five PC models with low outflow velocities (L3.85T24E85, L4.00T26E82, L3.85T28E85, L4.00T28E85, and L4.00T30E85), the dust stays in place and accumulates until the combined radiative pressure on larger amounts of dust is able to form a wind. In the corresponding drift models, the drift velocity reaches high values near where dust is first formed. The high drift velocities ablate dust grains, and instead of accumulating, dust grains move outwards through the gas and leave the gas without dragging it along. Our models are not setup to handle such cases where the outer boundary falls back towards the photosphere. Instead, the increased amounts of dust around the star heat up and affect the radiative transfer and overall physical structure in the enclosed star.
Six models show a situation where the mass loss rate and outflow velocity in the PC model become higher in the drift model (L3.85T24E91, L3.85T30E88, L4.00T30E91, L4.00T32E88, L4.00T32E91, and L4.00T28E88 that is calculated using Mie scattering). The more efficient dust formation in the drift model is able to form more dust that is able to drive the wind more efficiently than in the PC case.
The aim of our study is not to set the limits of wind formation, but we find that the dust mass loss rate M d < ∼ 2.0-5 × 10 −10 M ⊙ yr −1 in setups that fail to form a wind. It is tricky to make a better determination as the amount of formed dust also depends on if the drift velocity is high enough that grains are ablated by non-thermal sputtering. Also, a smaller value of C/O < ∼ 8.50 appears to make wind formation difficult; in our current set of models, only one such drift model forms a wind, at high luminosity and low effective temperature (L4.00T24E85).

Classification of periodic and irregular variations
The temporally averaged properties in Table 2 were calculated using different time intervals. Some models show periodic variations relatively quickly after calculations begin, whilst other never turn periodic. We use three variability classes: irregular (i), periodic (l×p), and nearly periodic or quasi-periodic (x×q). Variations of quasi-periodic winds are not perfectly periodic in both gas and dust properties or show a multiplicity that deviates from an integer factor of the pulsation period. Winds with a lower outflow velocity show a perfectly periodic radial structure with very lowamplitude variations that do not appear as periodic at the outer boundary -we still classify these winds as periodic.
Out of the 31 PC model winds in Table 2, 26 models are classified as either periodic (15) or quasi-periodic (11). The remaining five models are classified as irregular. The classification of the 22 drift models are different with 8 irregular structures, and 14 structures that are either periodic (10) or quasi-periodic (4). We show four temporal structures of both drift and PC models to illustrate the three classifications in Fig. 10, the two model setups are L4.00T28E88 and L3.70T24E88. The PC model of setup L3.70T24E88 is only evolved for a couple of periods after the structure turns quasi-periodic, of period 1.6P. Both PC models show a variability in the terminal velocity and mass-loss rate that change in amplitude and are also not integer multiples of the respective pulsation period, which is why both model structures are classified as quasi-periodic instead of periodic. The drift model L3.70T24E88 shows a clear irregular structure, and finally the drift model L4.00T28E88 shows a clear periodic structure.
The classification is occasionally a bit uncertain between the periodic and quasi-periodic classes. The classification might also change with longer modelling intervals. However, we believe that it is more important that the models are further developed with more physics, which could change the structure completely, before the assessments on this level of detail are attempted anew. We note, however, that the classification of all more numerically accurate models in Paper IV are stationary or periodic; this result is not reproduced here with our new models, currently.

Outflow velocity versus radiative acceleration
The outflow velocity (and mass-loss rate) is often related to the ratio of radiative to gravitational acceleration (of the dust; e.g. Lamers & Cassinelli 1999, chapter In M10 (equation 7), we define the wind-formation efficiency parameter α = δ dg L ⋆ M −1 ⋆ that is, in principle, proportional to Γ. Here, we adjust this parameter to account for the dilution of the dust component owing to drift (see equation 32) The wind-driving mechanism is stronger at higher outflow velocities. Our PC model values, as well as the associated values of M10, overlap the values of E14 well. Some of our PC model α values are higher than what both M10 and E14 find. It is more interesting to see that the drift-model values are shifted towards higher to drastically higher values of α. Whilst the relation shows a high correlation for PC models, the correlation is lower with drift models. The figure illustrates the drastically higher amounts of dust formed in drift models at increasing values of α In comparison to our study, E14 define four different pulsation classes for PC models that form winds: steady winds with small temporal variations (ws), winds with periodic variations in properties (wp), winds with more irregular non-periodic variations (wn), and winds that show an intermittent episodic outflow (we). Most models are nonperiodic, slightly fewer models are steady or episodic (at lower outflow velocities (u ∞ < ∼ 13 km s −1 ) and the remaining models are periodic (see their figure 5). The authors show a fact sheet for one model that they classify as periodic (wp; fig. C.1), which reveals a temporally variable structure that we would here classify as irregular. We believe our models show more accurate and periodic structures than darwin (Appendix E2). We find periodic variations across a larger region of the plot than E14; including all models with u ∞ < 17 km s −1 and all drift models with u ∞ < 25 km s −1 . Notably, our PC models that were calculated using Mie scattering are found at higher terminal velocities than all our other models (at the same α value). Figures 2 and 4 illustrate well defined ranges of physical values that result with our current set of wind models. Additional sets of models could likely extend the relations further towards both less and more massive winds.

Characterizing mean wind properties with F D
The figures reveal an exponentially decreasing dependence with F D in the mass loss rate M (Fig. 2a), dust mass loss rate M d (Fig. 2b), dust-to-gas mass loss ratio M d / M = δ dg F D (Fig. 2c), terminal velocity u ∞ (Fig. 4a),  and degree of condensation f cond (Fig. 4b); fits are shown in the respective figure panel. The mean drift velocity v D instead increases with F D (Fig. 4d). There is little correlation between the luminosity of each fit with the central star luminosity, except that the models using log L ⋆ = 3.85 seem to result in the best exponential fits. It is unclear that any similar relation exists for the mean grain radius r d (Fig. 4c).
Our current models reveal a maximum mass loss rate of M (F D = 1) ≈ 10 −5 M ⊙ yr −1 , where v D = 0 km s −1 . It appears that drift models do not form higher mass loss rates; optically dense models might play a role to understand the occurrence of such high mass-loss rates (see the discussion in IE10). And the mass loss rate decreases with increasing drift factor, down to M(F D ≃ 4.1)) ≈ 2×10 −7 M ⊙ yr −1 . Our models show a lack of lower mass loss rates, which is likely a result of the model parameters we have chosen for our calculations. But the results are also a function of the assumptions in form of grain properties and the interaction between the gas and dust.
The fit to the dust mass loss rate is most strongly correlated and it is also shows the steepest decrease with the drift factor. The range of values is 2.4 × 10 −10 (F D = 4.1) < ∼ M d < ∼ 2.5 × 10 −7 (F D = 1.0) M ⊙ yr −1 . Whilst dust formation in the form of grain growth increases with lower values of the drift velocity (cf. fig. 1 in Paper III), the grain growth quickly becomes less efficient at higher drift velocities where grains are also ablated by non-thermal sputtering; compare Fig. 3, which shows low dust mass loss rates in five out of six models where v D > 30 km s −1 . The exception is the high carbon-content model L4.00T32E91. At some point where v D > ∼ 40 km s −1 , there is a cutoff where dust is unable to drive a wind (see Section 4.2.4).
All values of the dust-to-gass mass loss ratio lie in the range 0.66 < ∼ 10 3 δ dg F D < ∼ 2.6, which is higher than most of the non-drift dust-to-gas density ratios δ dg that are reported in earlier studies (see figure 4 in M10, figure 5 in E14, and figure 8 in B19), also see Section 4.5. Non-drift models show too low values except with the highest mass loss rates when F D = 1.
The terminal velocity decreases with the drift factor, which indicates increasing difficulties at forming high outflow velocities when the drift velocity increases. The maximum outflow velocity in our models is about 65 km s −1 at F D = 1; notably, observations do not show such high values (see Section 4.5). All terminal velocities u ∞ > 39 km s −1 are found in the high carbon-to-oxygen-models, which are also not seen in observations. Moreover, the degree of condensation decreases with the power 1.3 of the drift factor and all values are found in the range 0.11 < ∼ f cond < ∼ 0.60. The mean grain radius shows no evident relation with the drift factor. Instead, we show the dust-to-gass mass loss ratio M d / M versus the mean grain radius r d in Fig. 5. The figure shows that the mean grain radius increases with the carbon-to-oxygen ratio.
The drift velocity reveals a trend of increasing values with the drift factor. The minimum value of our fit is v D (F D = 1) ≃ 12 km s −1 , whilst the models show v D > 10 km s −1 . Drift velocities are with few exceptions higher than in earlier higher accuracy Planck-mean models where 3 ≤ v D ≤ 5 km s −1 (Paper IV), and lower-accuracy Planckmean models where 4 ≤ v D ≤ 16 km s −1 and constantopacity models where 2 ≤ v D ≤ 13 km s −1 (Paper II). The upper limit on the drift velocity in these older papers is indicated with a horizontal guide in fig. 4d. Notably, the older results are calculated for a different region in the parameter space.

Replacing SPL with Mie scattering in the dust extinction
We show in Section 4.2.3 that the higher dust extinction of frequency-dependent models makes a difference compared to when a lower grey extinction is used. For the models presented in this paper, which have a mean grain radius in the range 0.070 < ∼ r d < ∼ 0.75 µm, the radiative pressure efficiency factor Q abs,ν (pr) is up to about a factor 5 higher when Mie scattering is used instead of SPL (see figure 3 in MH11); this implies that effects can be even stronger compared to when using grey extinction values, and drift velocities can be even higher. We calculated three sets of models for setups L3.70T28E88, L3.85T28E88, and L4.00T28E88 using Mie scattering to see how the outcome is affected.
Compared to the results of the respective SPL model, the PC model values are 8.8-38 per cent lower ( M , f cond , M d / M ). Differences are smaller in the mean grain radius r d , 13 per cent lower to 3.5 per cent higher. However, the outflow velocity u ∞ is 32-49 per cent higher. Compared to the values of the model of M10, the grain radius r d is 55 lower to 56 per cent higher, the outflow velocity u ∞ 13 per cent lower to 81 per cent higher, and the remaining properties 17-82 per cent lower.
The changes are high in the two drift models as well where the drift velocity v D is 35 per cent higher in model L3.85T28E88 and the mass loss rate M is 19 per cent higher in model L4.00T28E88 than in the corresponding SPL model. The remaining properties are 24-79 per cent lower; additionally, the periodic structure of the SPL model L4.00T28E88 turns irregular using Mie scattering. Compared to the values of the PC model of M10, the values of L3.85T28E88 are 7-67 per cent lower. The differences are seemingly lower for model L4.00T28E88 where the degree of condensation f cond is 56 per cent lower whilst the other properties are 6.9 per cent lower to 14 per cent higher. No wind forms using Mie scattering with setup L3.70T28E88.
We show the radial structure of the PC and drift models of setup L3.85T28E88 that use Mie scattering in Fig. 12. The drift model shows a more stationary appearing structure in the gas density and velocity structures  than the corresponding SPL model (not shown, but compare with the similarly appearing drift-model structure in Fig. 7). The moderately high and nearly constant value of the drift velocity results in little ablation (Fig. 12c). Whilst the temperature ratios are about the same as for model L3.70T28E88 (Fig. 7k), higher values of the extinction are seen in Fig. 12h; the extinction ratio is about twice as high at the outer boundary and also shows a radial structure.
Considering the amount of formed dust and compar- T d /T r (T g /T r ) eq T ratio Figure 12. Radial structure of a snapshot of setup L3.85T28E88 using Mie scattering. The drift (PC) model is shown with purple (orange) lines. From the top left, the eight panels show: (a) gas velocity u, sound speed c s , and self-similar gas velocity u ss ; (b) gas density ρ g , dust density 10 4 × ρ d (log); (c) drift velocity v D , equilibrium drift velocityv D , and thermal velocity v ζ ; (d) degree of condensation f cond ; (e) net growth rate τ −1 gr , net decay rate τ −1 dc , nucleation rate J ⋆ , and non-thermal sputtering rate τ −1 ns (log); (f) average grain radius r d ; (g) temperature ratios T d /T r and (T g /T r ) eq ; and (h) extinction coefficient ratio χ H /χ. All properties are drawn as function of the stellar radius R ⋆ (lower axis) and astronomical units (AU; upper axis). Grey horizontal lines are guides.
ing with the values of the SPL models (Fig. 11), the three PC models form less dust and attain drastically increased outflow velocities, which places the models above and away from all other model values. The drift models also form less dust than when SPL is used, but in combination with the lower outflow velocities the new values are closer to the other model values.
In our comparison of results of Mie and SPL models in MH11, we find for two sets of models that effects are rather small; model values are both lower and higher than when SPL is used, although we see a tendency towards smaller amounts of formed dust in all models that use Mie scattering. Here, we have modelled three setups using Mie scattering -which is why it is too early to generalise the results to all other models. However, the lower outflow velocity achieved with Mie scattering places the two drift models in a region that is closer to values of observations (Fig. 13). The large changes, which are larger than we find in the PC models of MH11, indicate that, in models that use a high spatial resolution, Mie scattering is an important process that needs to be used in place of SPL in all models of carbon-rich mass loss rates.
Models that are calculated using Mie scattering result in lower outflow velocities, when models include drift, which also implies higher drift velocities. The models using Mie scattering follow the same trends in properties as the SPL models when the properties are related to the drift factor.

Radio and infrared observations
One approach is to measure mass loss rates (of the gas) using radio observations of emission lines of CO and apply the stationary wind approach of Morris (1980) and Knapp & Morris (1985), see, for example, Olofsson et al. (1993), Knapp et al. (1998), Schöier & Olofsson (2001 derive their values using an improved RT model), Schöier et al. (2002), and Groenewegen et al. (2002a,b).
A more common approach, currently, is to measure spectra in the infrared wavelength range for both C-rich and O-rich stars at a known distance [typically the Magellanic Clouds (MC)], fit a spectral energy distribution (SED) to the spectrum and use the dust optical depth to extract mass loss rates of the gas (e.g., van Loon et al. 1999;Groenewegen et al. 2007Groenewegen et al. , 2009 or dust (e.g., Jura 1986;Zijlstra et al. 1996;Srinivasan et al. 2009;Sargent et al. 2010Sargent et al. , 2011Srinivasan et al. 2011;Jones et al. 2012Jones et al. , 2014Riebel et al. 2012;Nanni et al. 2019). The chosen approach is to -amongst other assumptions -use a fixed outflow velocity of the dust v ∞ = 10 km s −1 and fix the dust-to-gasratio δ dg = 0.005; the drift velocity is also assumed to be zero, so u ∞ = v ∞ .
The parameter selection is important. The results of our drift models -that we calculated using solar metallicity opacities -indicate that it is very hard to reach both the highest mass loss rates and low outflow velocities (cf. Fig. 2a and 4a). For example, M = 1.0 × 10 −5 M ⊙ yr −1 requires u ∞ = 75 km s −1 , which is a high value that is not observed. When we instead use u ∞ = 45 km s −1 , v D,∞ = 17 km s −1 (Fig. 4c), M = 3.9 × 10 −6 M ⊙ yr −1 (Fig. 2a), and δ dg = 6.8 × 10 −3 (Fig. 2c). Please note, however, that the error bars are significant! Considering the conditions where u ∞ = 10 km s −1 -as in the lower metallicity MC-specific parameters mentioned above -the drift velocity v D,∞ = 55 km s −1 , the mass loss ratio M = 4.9 × 10 −6 and the dust-to-gas ratio δ dg = 1.1 × 10 −4 (Fig. 2c). These values are naturally not directly applicable to the MCs; instead one should use relations based on models that are calculated using MC-specific metallicities). Differences are large, the dust-to-gas ratio assumed by the authors is 50 times higher than this value; in good agreement with the range of values B19 present for the MCs. Drift velocities are, moreover, never negligible, and it appears plausible that they are even higher in low-metallicity environments.
The gas mass loss rate can be derived in the SED approach without knowledge of the drift velocity, assuming it is possible to estimate the dust density ρ d and dust-togas density ratio δ dg . However, to measure the dust mass loss rate M d , it is necessary to estimate the drift velocity v D,∞ . Nearly all mass-loss estimates based on SED data ignore drift altogether. Nanni et al. (2019) mention drift, but ignore it based on the results of the grey and stationary model of Krüger & Sedlmayr (1997), who include grain growth and drift and find that drift velocities are very small ( v D,∞ < ∼ 5 km s −1 ). The systematically higher mass loss rates found by Groenewegen et al. (2002b) in their CObased study imply lower drift velocities, v D,∞ ≃ 3 km s −1 ; Schöier et al. (2002) find the same value. Groenewegen et al. also point out that their observations of gas-to-dust ratios are in good agreement with theory, but note that the cited theoretical studies are all based on grey gas opacities and as their own study, none considers effects of drift. There appear to be two exceptions to low drift velocities: based on stationary models, Papoular & Pégourié (1986) and Olofsson et al. (1993) estimate higher drift velocities that are more similar to our model values.
Average drift velocities in stellar winds can be estimated with time-dependent models such as those that are calculated here. Alternatively, as a first approach, one can assume equilibrium between the radiative pressure and the drag force, as in equation (30). Physical assumptions regarding both the gas and the dust influence the result. The dust velocity is required for accurate estimates of yields of dust, which are often provided as dust-to-gas (equation 32; or gas-to-dust) ratios.

Model values of M10, E14, and observations
In their comparison of results of C-rich models using darwin, E14 plot observed values along with model values of the respective mass-loss rate and outflow velocity. We re-create this figure with Fig. 13, and consequently include the C-star observations of Netzer & Elitzur (1993, who gather values of 18 references), Knapp et al. (1998), Schöier & Olofsson (2001), and Groenewegen et al. (2002b). We also overplot our new values as well as the delimited region of lower optical depth drift-dominated mass loss according to IE10 (equa-tion 10); we use their fig. 3 to set the delimiting values 10 ≤ A ≤ 30 and they take the data from Olofsson et al. (1993). Reddening dominates in the region with higher mass-loss rates (above the upper line), and here drift is of less importance. In this plot, we omit our models using log (C − O) + 12 = 9.1, with one exception, as they achieve very high outflow velocities.
Observations of Knapp et al. (1998) and Schöier & Olofsson (2001) mostly lie within the limits of drift-dominated models, whilst more values of Groenewegen et al. (2002b) lie in the region that is classified as reddening-dominated based on the data of Olofsson et al. (1993); mass-loss rates of six objects that overlap in these two studies agreee to within a factor two when using the same object distance (M. Groenewegen, priv.comm. 2020) and three (two) of these six objects are shifted out to the redshift-dominated (in to the driftdominated) region when using the distances of the more recent study. The mixed-origin values of Netzer & Elitzur (1993) are spread out over most of the plot where there are values of observations.
A majority of the PC models of E14 are found in the drift-dominated region, with the exception of models with very low outflow velocities, which are also in the reddeningdominated region. Except low-velocity "episodic" model values, remaining model values are found in the strip where −6.5 < ∼ log 10 M < ∼ − 5.0 M ⊙ yr −1 . Few values are found in the region −7.5 < ∼ log 10 M < ∼ − 6.3 M ⊙ yr −1 , 3 < ∼ u ∞ < ∼ 14 km s −1 , which combines lower outflow velocities and mass loss rates.
Our new PC model values, which use more gridpoints, more frequency points in the RT, and improved numerical features, differ from the original values of M10 throughout the plot. Differences are the largest in the drift dominated region, where outflow velocities change by −87 to +81 per cent and mass loss rates by −97 to +5.3 per cent. Few of our new values change so much that they enter the region of values of Knapp et al. (1998) and Schöier & Olofsson (2001). All our models that use log (C − O) + 12 ≤ 8.5 are found in the reddening-dominated region at velocities higher than in the drift-dominated region.
Values of our drift models are both higher and lower compared to the respective PC model; changes in the outflow velocities and mass loss rates are with one exception −67 to +48 and −58 to +38 per cent, respectively. The changes are the largest in model L4.00T32E88, +380 and +1100 per cent, respectively. The correspoinding changes relative to the PC model values of M10 are −69 to +56 per cent and −83 to −3.2 per cent. The PC models with the lowest terminal velocities ( u ∞ < ∼ 8 km s −1 ) and where log(C − O) + 12 = 8.2 do not form a wind. The drift models L3.70T26E88 and L4.00T24E85 lie in the reddening-dominated region. Model setups L3.85T30E88 and L4.00T32E88 are two cases near the lower limit of the drift dominated region where both the terminal velocity and the mass-loss rate of the drift models are higher than in the corresponding PC model. Most of the observations of Groenewegen et al. (2002b) show higher values that are not covered in our parameter space.
There is a deficiency of drift models with a lower outflow velocity. Except the three models in the reddening dominated region and model L3.70T28E88, there are no (SPL) drift models with v D < 17 km s −1 . Notably, when model L3.85T28E88 uses Mie scattering instead of SPL, the outflow velocity decreases by 39 per cent to u ∞ = 12.7 km s −1 , and the mass loss rate decreases by 46 per cent (also see Section 4.4). It appears that a study is valuable where all models use Mie scattering to get model values that agree better with observations. Effects of drift dominate in the optically thin region, as predicted by EI01, and drift will be an important component when finding model setups that match the observations in this region. (More difficult to model the high mass loss rate values in the optically thick region.)

Comparison with the theory of EI01
Instead of solving the formidable wind formation problem described here, EI01 present a solution that is much less numerically demanding and yet yields general properties of dust driven stationary winds and grain drift. A comparison with our results is justified.
The main idea of the work of EI01 is that properties of the wind at large distances are expressed using attributes in the wind formation region. The assumptions are the following. All dust forms promptly at the grain condensation radius r c and grains neither grow nor are destroyed outside of this radius. Grains are described with their type, assuming a constant size, condensation temperature, and absorption and scattering efficiencies. The pressure gradient is ignored as winds of interest are assumed to be highly supersonic. The star is losing mass at specified mass-loss rate at the condensation radius r c = 2R ⋆ . The authors present results in the form of radial density and velocity structures, structures plotted versus a drift-effect parameter that is the ratio of radiation pressure to drift effects (P), and reddening correction factors plotted versus the optical depth τ V (at the visual wavelength 550nm; EI01, equation 17).
The winds of our study with the lowest (L3.70T28E88) and highest (L4.00T24E91) dust mass loss rate yield optical depths in the range 0.65 < ∼ τ V < ∼ 15. For the discussion here, we show temporally averaged radial structures of the drift model L3.70T28E88 in Fig. 14. In comparison to Fig. 7, the properties show larger variations here owing to a shorter time interval (5P) and low numbers of models (54) used to calculate the averages. Moreover, we only consider values where r ≥ r c .
The velocity structure (Fig. 14a) agrees well with the self-similar velocity structure u ss of IE10 (their equation 1 and our equation 35) when we use the gas velocity at the outer boundary as u ∞ . The agreement is less optimal when  Figure 14. Temporally averaged radial structure of a set of snapshots of the drift model L3.70T28E88 showing calculated properties of EI01. From the top, the five panels show: (a) gas velocity u, sound speed c s , self-similar gas velocities u ss and u ss,∞ ; (b) average grain radius r d ; (c) drift profile EI01 ζ ; (d) calculated (EI01, equation 18) and parameterised (EI01, equation C7) dust density profiles EI01 ηy 2 = EI01 η(r/r c ) 2 ; and (e) reddening correction factors EI01 K 1 and EI01 K 2 that are used in the circled region at the condensation radius r c . See Fig. 7 for more details.
we instead use equation 35 in EI01 to calculate u ss,∞ and their expressions for u ∞ 5 . The resulting velocity structure u ss,∞ is about 2.5 times higher than u ss . Furthermore, the drift profile EI01 ζ (u/v = F −1 D , equation 12 and figure 2 in EI01) shows the same trend as in EI01 (Fig. 14c). Our model value lies closer to the τ V = 0.08 line of EI01. This property is always less than 1 when the model accounts for gas-to-dust drift.
The density profile EI01 η (see Fig. 14d; compare equation 18 and figure 2 in EI01) is very similar to what EI01 show. Admittedly, our line is shows rather large variations owing to the small number of models used in the shown average (see above). The trend of the parameterised density ( EI01 η calculated using equation C7 in EI01) is very similar to EI01 η. The agreement of the two profiles supports the statement that the density profile does not depend on drift (EI01).
We also show the reddening correction factors EI01 K 1 and EI01 K 2 (equations 48 and 49 and the top two panels in figure 4 in EI01) in Fig. 14e; these corrections are scalar properties that apply to the wind formation point r c . We show radial profiles of these structures to emphasize that as the wind formation point is not a fixed spatial coordinate in our models. According to its definition, the reddening correction EI01 K 1 (τ V = 0.7) ≃ 1.1 and EI01 K 1 ≥ 1 ∀ τ V . We find a somewhat lower value, EI01 K 1 ≃ 0.3. Finally, we find EI01 K 2 ≃ 0.3, which is also lower than the value according to its definition, EI01 K 2 (τ V = 0.7) ≃ 0.95 and EI01 K 2 < 1 ∀τ V . Both values are, however, in qualitative agreement with EI01, and in particular EI01 K 2 indicates a stronger need for significant reddening corrections for the lower optical depth value of our model.
The authors show that the mass-loss rate depends on the optical depth as τ 3/4 V (equation 55) owing to drift. Considering how influential we find that the inclusion of drift is in this study, it would be interesting to see how mass-loss rates derived in infrared observations are affected if they would use this relation instead of a usual relation that is linear with τ V . Moreover, according to the theory of EI01 and accounting for effects of drift, the mass-loss rate is proportional to the outflow velocity as M ∝ u ∞ 3 and it is not related to the stellar luminosity. When drift is ignored, the outflow velocity is instead proportional to the luminosity as u ∞ 4 ∝ L ⋆ . With a limited number of exceptions, the mass-loss rate of the PC models of M10 (figure 1), E14 (figure 4), B19 (figure 6) show a non-existant to very weak dependence on the outflow velocity (all model values appear in a wide band), seemingly in agreement with the finding of EI01. Our drift models follow the expected behavior in the sense that they, with two exceptions, lie inside the driftdominated region indicated in Fig. 13. Notably, we find that models calculated using Mie scattering show a trend of both lower outflow velocities and mass-loss rates. It seems plausible that all models calculated using such scattering lie inside this region. Our three PC models calculated using Mie scattering show increased drift velocities and appear to follow the relation u ∞ ∝ L ⋆ .
Our models are time-dependent and not stationary, and drifting dust initiates the stellar wind and dust grains form wherever conditions are suitable in a pulsating atmosphere. The two approaches show a qualitative agreement, but the physical differences result in quantitative differences as shown.

CONCLUSIONS
We have extended our grey dust-driven wind models of Papers I-IV to include frequency-dependent RT in both the gas and dust components. We have also rewritten the RT solver and included our improvements to the numerical description of advection terms and the discretisation scheme. With our new model code T-800, we are able to model more realistic configurations of stellar winds of carbon-rich stars than was possible in our earlier papers; we can use the same physics as the other stellar wind model code that is currently used, darwin (Höfner et al. 2016) -with the added advantage that we can also model effects of gas-to-dust drift (two-fluid flow) and use high spatial resolution. To our knowledge, this is the first time a cool stellar wind is modelled at this high level of physical detail. And it appears that effects of drift become stronger with the level of detail.
Based on the set of models of M10, we have calculated both PC (non-drift) and drift models to clarify differences owing to drift. As in Paper IV and in contrast to the approach of darwin, we have skipped the adaptive grid equation to resolve shocks in the gas. Instead of the commonly used 100 gridpoints, our models use 1024 gridpoints that are fixed in space where they resolve shock fronts in the gas and fronts in the dust simultaneously, at all times. Here, we have mostly calculated dust extinction rates using the SPL, but we used Mie scattering in a few cases for comparison. Benchmark results of T-800 and darwin using the same PM model setups shows reasonable agreement.
We have found periodic variations in a majority of the calculated structures of both PC and drift models, whilst earlier studies have reported irregular variations; we attribute this difference to the higher spatial resolution and improved numerical accuracy in our models as we also achieve irregular structures when we use the same numerical and physical setup as those studies. The results reveal intermediate to large changes of 50-1000 per cent in properties such as outflow velocities and mass-loss rates when we compare with our PC models in M10. Outflow velocities, in particular, appear to increase greatly when PC models use Mie scattering instead of SPL; changes are greater than we find in MH11.
Calculating dust yields, the dust-to-gas density ratio is multiplied with the drift factor. Consequently, our drift models yield 1.2-4.2 times as much dust as a corresponding PC model that ignores drift, and where the dust-to-gas density ratio is unchanged; this result corroborates our find in Paper III and Sandin (2003) that dust formation is more efficient owing to drift. It is impossible to get correct yields when drift is ignored. Furthermore, our results show that ablation is unimportant, grain destruction rates are too small at the large radii where drift velocities are high enough to activate ablation by non-thermal sputtering.
With two exceptions, our drift models lie in the interval of drift-dominated outflows as discussed by EI01 assuming M ∝ u ∞ 3 (Fig. 13); this interval is in turn based on the CO-based observational data of Olofsson et al. (1993). We could not calculate corresponding drift models for many PC models that lie outside of this interval. And consequently, most of the high-mass loss observations of, for example, Groenewegen et al. (2002b) also lie in a different parameter space than our models. Nearly all observational studies disregard drift and in view of our results therefore achieve too small yields of dust, in particular for lower mass loss rates where M < ∼ 1.0 × 10 −5 M ⊙ yr −1 . A comparison between the results of our timedependent models and the simplified theory developed by EI01 shows a qualitative agreement, but the simplifications prevent a quantitative agreement. We agree with EI01 that drift is always present and an important factor of the wind when calculating models that result in realistic outflow velocities that match observations.
Our current dust-driven wind models do not reproduce the combination of higher mass loss rates -M > ∼ 1.0 × 10 −5 M ⊙ yr −1 -at low expansion velocities. The models can, of course, be improved further with more physical detail. We could add a set of additional equations of motion of the dust to describe binned size-dependent drift velocities. Considering the already high complexity of our models, it is difficult to predict what the effects of such a treatment would be. Plausibly, particles of some sizes would move faster through the gas, which might increase ablation rates. Moreover, the current set of models need to be calculated using Mie scattering only, SPL is too inaccurate in the grain size interval we model. It would also be of interest to calculate models at lower metallicities to see how drift affects results in such circumstances (Mattsson et al. 2015;B19). In a longer perspective, it would be highly valuable to extend T-800 with oxygen-rich chemistry including silicates to include drift in models of M star winds.
son Lindblom at PDC for assistance concerning technical aspects in making T-800 run on the PDC resources.

DATA AVAILABILITY
The data underlying this article will be shared on reasonable request to the corresponding author.  Avogadros constant N C number of rays inside the model n C cm −3 gas phase total number density of condensible material N D number of gridpoints N ν number of frequencies used in the RT n ν refractive index, phase velocity n d cm −3 dust number density, n d ≡ K 0 N l lower size limit of macroscopic grains ν s −1 frequency P g dyn cm −2 gas pressure q sphericality Q abs,ν χ J , κ S : Q abs,ν = Q ext,ν − Q sca,ν , χ H : Q abs,ν = Q ext,ν − cos θ ν Q sca,ν χ R : Q abs,ν = Q ext,ν , as Q sca,ν = 0 Q ext,ν extinction efficiency Q sca,ν scattering efficiency Q ′ X,ν cm −1 absorption / extinction efficiency, Q ′ X,ν = Q X,ν /a gr r cm radius r 0 cm monomer radius R ⋆ cm stellar photosphere radius r c cm dust condensation radius r d cm mean grain radius, gas density ρ m g cm −3 grain intrinsic density S cm −3 s −1 amount of condensed material S g erg cm −2 s −1 gas source function σ cm 2 mean grain cross section, σ = πr 2 0 K 2 /K 0 σ grain erg cm −2 grain surface tension K grey dust temperature, T d,grey = T r T g K gas temperature T r K radiation temperature, T 4 r = J π/σ SB T r,ext K radiative temp. at the outer boundary τ −1 s −1 net grain growth rate, total grain decay rate τ −1 ns s −1 total non-thermal sputtering rate τ V optical depth of the dust, at 550nm u cm s −1 gas velocity u ss cm s −1 self-similar gas velocity using sticking coefficient setup ξ atom and molecule-specific sticking coefficient χ cm −1 grey extinction coeffient, χ = πr 3 0 K 3 × 4.4T d,grey χ ν cm −1 Hz −1 dust extinction coefficient χ J,H,S cm −1 frequency average of χ ν , weighted by J, H, and S g , respectively χ R cm −1 Rosseland mean extinction coefficient relative fluctuation ampl.,r = σ s /Q Q temporal mean of the property Q σ s standard deviation u ∞ km s −1 terminal gas velocity v ∞ km s −1 terminal dust velocity ordinary differential equations (ODEs), where the four primary variables are: m r , P g , H, and T g .
The initial model code of T-800, John Connor, solves the four ODEs in equations (B1)-(B4) together with four conditions that specify the stellar mass M ⋆ , the effective temperature T eff , the luminosity L ⋆ , and the oxygen and carbon abundances ǫ O and ǫ C [C/O = log (C − O) + 12 = log (10 ǫ C − 10 ǫ O )], as well as the mean molecular weight µ. Additional input is required in form of the radial extent of the model domain ([r int , r ext ]), frequency-dependent opacities, required accuracy of the calculations, etc.
John Connor calculates an initial model in the form of a hydrostatic stellar atmosphere of an AGB star in four steps: (i) Calculate a model using grey RT. The calculations begin at the photosphere radius R ⋆ with preset values of the mass, effective temperature, and luminosity. The photosphere radius is taken from the relation (B5) Using an initial guess of the pressure, the equations are integrated outwards to the outer boundary. This integration is iterated by adjusting the photosphere pressure until the temperature at the outer boundary is lower than the temperature T g in the equation where a value at the outer boundary is denoted with the subscript "ext", and the angular intensity distribution of the radiation fieldμ = (H/J) ext is taken from the RT calculations; initially,μ = 0.5. Thereafter, the ODEs are integrated inwards from the photosphere to the inner boundary, without solving the RT equation whilst assuming f Edd = 1/3; this step provides initial values for all four primary variables at the inner boundary. So far, this procedure is very similar to the description given in Dorfi (1998, chapter 4.10). Finally, the ODEs are integrated from the inner boundary to the outer boundary whilst calculating RT. The pressure at the inner boundary is again iterated to achieve the estimated temperature at the outer boundary. This step makes use of an ODE solver, such as the implicit predictor-corrector ODE solver slga (Raith et al. 1979), which handles stiff equations.
(ii) Calculate a model using frequency-dependent RT based on the primary variables at the inner boundary of the first step as a starting point. The ODE integration is anew iterated to achieve the estimated temperature at the outer boundary, which will be different compared to the grey calculations.
(iii) Calculate temperature corrections to achieve radiative equilibrium using the inner boundary values of the second step as starting point. The temperature corrections are calculated for the full radial extent using the Unsöld-Lucy procedure and the description in MH15 (chapter 17.3). Consequently, three ODE:s are solved in this step, equations (B1)-(B3), making use of the ODE solver dlsode 6 . Temperature corrections are applied iteratively until the maximum temperature and luminosity corrections are both below a preset threshold. The pressure at the inner boundary is anew iterated until the temperature at the photosphere radius is close enough to the specified value T eff .
(iv) The resulting physical structure is relaxed on the adaptive grid, also when all grid weights are set to zero. Regardless of how many gridpoints the ODE solver has used, the solution is interpolated to use all or a part of the N D required gridpoints, the latter option is the case when using a fixed grid with a larger extent than the initially modeled hydrostatic region. The gridpoint concentration in the centremost region is doubled relative to the remaining domain according to the description in Section 2.4.1. The physical structure is stored in a binary file for use with T-800.

APPENDIX C: THE ROLE OF STICKING COEFFICIENTS
We wanted to compare our new results with those of M10 who use higher sticking coefficients. Here, we also calculated a set of PC models using the adaptive grid equation and N D = 100; these models use the RT of Yorke (1980) and 64 frequencies, as well as the same second order van Leer advection.
We present our results of the three models that we calculated using the two sets of sticking coefficients Ξ 0.34 and Ξ 1.00 in Table C1. We show the results in Fig. C1 drift velocity of model L3.85T28E88, and the mass loss ratio of model L4.00T28E88.
Except that smaller amounts of dust forms in models that use Ξ 0.34 , there is no simple relation between property values and sticking coefficient. In this paper, we use the set of sticking coefficients suggested by Gail et al. (1984), Ξ 0.34 .

APPENDIX D: PHYSICAL AND NUMERICAL SETUP
We collect values and references of most physical and numerical parameters and assumptions in Table D1. Parameter values of the benchmark test models in Appendix F are shown in column 3, separate from the parameter values used in the other parts of this paper in column 2.

APPENDIX E: NUMERICAL ISSUES IN DRIFT MODELS E1 Troughs in the drift velocity
Our first generation of drift models often show high-value spikes in the drift velocity at the front of dust "shocks", see for example figure 2 in Paper I. Such spikes occur when the relative speed between a moving gridpoint (of the adaptive grid equation) and the dust is very close to zero. When this happens, the numerical diffusion may become too small to prevent spurious results to appear (see Paper I, section 3.2). The problem is mitigated in the Planck-mean models presented in Paper II-Paper IV, where both densities and drift velocities are lower than in the constant-opacity models of Paper I.
Densities and drift velocities of our new frequencydependent RT models presented here are higher than in the Planck-mean models. And if gridpoints are allowed to move about using the adaptive grid equation, conditions are anew favourable for the appearance of spurious spikes in the drift velocity. We decided to attempt to avoid these spurious spikes by keeping the grid fixed, which ought to work for as long as the dust velocity satisfies v 0 km s −1 . As a compromise, we kept all gridpoints fixed where dust is present, r > 2R ⋆ , while gridpoints at smaller radii were still allowed to stretch with the movements of the piston at the inner boundary.
Our new drift models are free of spikes in the dust component owing to vanishing numerical diffusion. However, most models show "troughs" of lower -and even negativedrift velocities where there are strong negative radial gradients of the dust density (in dust fronts), for example see r ≃ 15, 21, 27, 33, and 38R ⋆ in Figs. 7b and 7c. A term that is susceptible for the origin of these features is the advection in the equation of motion of the dust (equation 8).
VanderHeyden & Kashiwa (1998) present an improved van-Leer-type ("compatible") advection scheme that unlike previous schemes is designed to preserve a monotone character of both the momentum and the velocity in an advected property such as the dust momentum (ρ d v). We performed tests where we replaced our regular advection scheme with this one, but we could see no improvement where the troughs disappear. We do not currently understand the origin of these troughs, numerical or physical, but they occur where there are minute amounts of dust and their influence on the dynamics appears to be insignificant.

E2
The role of numerics to properties of the model variability A majority of both our PC and drift models presented here as well as in Paper IV result in periodically varying structures. A vast majority of the models of darwin meanwhile result in irregularly varying structures (S. Höfner, priv.comm. 2020). As we do in Paper IV, we attribute the occurrence of periodic variations, as opposed to irregular variations, to our use of a high numerical accuracy in our models owing to the volume-weighted advection scheme and volume-weighted averaging on the staggared mesh. Additionally, we use a higher number of gridpoints (N D = 1024 instead of N D = 100) and the adaptive grid equation is kept fixed where r > 2R ⋆ and is not used to resolve shocks anywhere. We calculated a few additional PC setups of model L3.85T26E88 for comparison where we used the same approach as in Paper IV and allowed all gridpoints to stretch with the movements of the piston at the inner boundary. The model calculated using the PPM advection scheme and arithmetic averages instead of the volume-weighted van Leer scheme and volume-weighted averages does not become periodic in the first 77 periods; the models are in all other aspects identical. Likewise, the model where the adaptive grid is not fixed at r = 2R ⋆ also fails to develop periodic variations, as does a model that uses the regular van Leer advection scheme. Periodic variations put high demands on the numerical method, which is why we advocate using as many gridpoints as possible, volume-weighted averages, volumeweighted (van Leer) advection, and fixing the adaptive grid wherever possible to avoid the degradation of the solution when all regions are not equally resolved and gridpoints restructure.

APPENDIX F: BENCHMARK TEST: darwin VS. T-800
We have compared the outcome of PC models of T-800 and darwin in a benchmark test where we used the same physical and numerical setup -as far as we know -as the models calculated using darwin; the parameter-value setup is shown in Table D1 on the right-hand side of the values used in the remaining parts of this study. We show the model parameters of this sample in Table F1 together with the corresponding model values of M10 and E14; additionally, we used M = 1.0M ⊙ and ∆ u p = 4 km s −1 . In these models, which use the adaptive grid equation to resolve shock fronts, the outer radius is at first allowed to move outwards from the location of the outer boundary of the hydrostatic inital model to r = 25R ⋆ , before the wind is calculated. We evolved the models for a period of 200-1000P. We follow the approach of E14 and calculated temporally averaged values for the interval that begins when the initial transient has left the model domain and ends after 1000 periods or when 20 per cent of the mass in the envelope remains, and we use the same properties as with the other models in this study. We show the results of our comparison in Table F2.
Our model values mostly agree with the values of both Table C1. Properties temporally averaged at the outer boundary varying the sticking coefficients. From the left, the first four columns specify: the name; the model type, PC (P) or drift (D); sticking-coefficient setup (Ξ); and the number of gridpoints N D . Remaining columns are the same as in Table 2. All models have been calculated with the outer boundary fixed at r ext final = 40 R ⋆ (N D = 1024) and r ext final = 25 R ⋆ (N D = 100). Rows of drift models are shown in boldface. The four last lines show values of Andersen et al. (2003,  E14 and M10 within 30 per cent, and the agreement is typically better than that. Our terminal velocities for models L3.70T26E85 and L3.85T26E85 are 90 and 120 per cent higher than the other two studies, and our mass-loss rate of model L3.70T26E88 is about 60 per cent higher. Whilst both T-800 and darwin use the same modelling approach as well as the same modelling parameters, as far as we can fathom, there are still differences in the discretization. For example, for one model, L3.85T28E85, we got values that were off by 40-60 per cent from our other set of values when we changed the discretization to subtract decay rates from growth rates before integrating the values over the volume of each grid cell (not shown in Table F2). Additionally, the amount of mass in the model domain of these models is smaller than in our other models, compare the last column in Tables F1 and 1. The implication of the smaller amount of mass is that the model domain is emptied of mass faster -the selection of the temporal interval used to calculate the mean values is meanwhile crucial. Add to that, in some cases we find it difficult to estimate the when the initial transient has left the model domain. The larger dif-ferences in the outflow velocities plausibly originate in some other difference.
The agreement overall is decent to good, and we conclude that T-800 is able to reproduce results of darwin. Notably, T-800 does not [yet] include a description of oxygenrich chemistry, which is why we cannot reproduce such models. Gilman (1972) is the first to study the degree of momentum coupling between dust and gas. He finds that the momentum coupling is instantaneous and complete. Later, Berruyer & Frisch (1983) show that the dust decouples from the gas at large radii (r > ∼ 1200R ⋆ ). MacGregor & Stencel (1992) find that complete momentum coupling holds when grains are large (a gr > ∼ 0.1 µm) and that the gas and dust phases decouple with small grains (a gr < ∼ 0.05 µm). Our models are dynamic instead of stationary, include both grain growth and ablation with dust grains that show different (mean) sizes, as well as radiative transfer. The models are more complex than before, and we think a comparison is justified.

APPENDIX G: ON COMPLETE MOMENTUM COUPLING
We use the approach of MacGregor & Stencel (1992) and compare the magnitudes of the following four terms in the dust equation of motion (equation 8): the inertial and temporal term f in , the gravitational force f grav,d , the radiative pressure on the dust f rad,d , and the drag force f drag . The inertial term is where we discretized the advection term using the simpler arithmetic mean expression of van Leer. A comparison of the force terms for our set of models shows a mostly complete momentum coupling, where f drag ≃ f rad,d and f in < f drag . However, we see some support for a relaxed coupling in the outer parts. We plot the radial structure of the force terms for the model with the lowest dust mass loss rate L3.70T2800E88 (Section 4.2.1) in Fig G1; we show all other properties of this model considered here in Fig. 7. The figure shows a nearly full coupling where f drag ≃ f rad,d − f grav,d and f in < f drag at all radii, except at dust fronts (e.g. r ≈ 10 R ⋆ ) where the difference is smaller. There is no clear decoupling as MacGregor & Stencel (1992) find (cf. fig. 4), plausibly because grains are larger. The influence of decoupled phases ought to be small, as forces are weak at the large radii where decoupling occurs. Effects might be stronger at much larger radii, as Berruyer & Frisch 1983 find.
We also show the radial structure of model L3.85T30E88 where the decoupling is somewhat stronger, see Fig. G2. Here, f drag ≃ f rad,d − f grav,d and f in is also similar to the other terms for r > ∼ 24R ⋆ , indicating more decoupled phases. As expectedv D ≃ v D for all radii, whilst there are some differences in the same outer parts.
Complete momentum coupling appears to be a suitable approximation in our models. It seems justified, however, to check this condition anew in future models where grains of different size move at separate drift velocities. Table D1. Physical and numerical assumptions Table F2. Temporally averaged quantities at the outer boundary; see Appendix F. From the left, the first two columns specify the model name (see Table F1) and a reference for the values (when other than our own). Five columns show the averaged: mass loss rate M , terminal velocity u ∞ , degree of condensation f cond , dust-to-gas ratio δ dg , and dust radius r d . The final column shows our outflow classification class (irregular i) and the one of E14.  Figure G2. Radial structure of a snapshot of the drift model L3.85T30E88 for the full modelled region. Panels a, b, and c correspond to panels a, c, and b in Fig. 7. Panel d corresponds to Fig. G1.