Alfv\'en wave-driven wind from RGB and AGB stars

We develop a magnetohydrodynamical model of Alfv\'en wave-driven wind in open magnetic flux tubes piercing the stellar surface of Red Giant Branch (RGB) and Asymptotic Giant Branch (AGB) stars, and investigate the physical properties of the winds. The model simulations are carried out along the evolutionary tracks of stars with initial mass in the range of 1.5 to 3.0 $M_{\odot}$ and initial metallicity $Z_{\rm ini}$=0.02. The surface magnetic field strength being set to be 1G, we find that the wind during the evolution of star can be classified into the following four types; the first is the wind with the velocity higher than 80 km s$^{-1}$ in the RGB and early AGB (E-AGB) phases; the second is the wind with outflow velocity less than 10 km s$^{-1}$ seen around the tip of RGB or in the E-AGB phase; the third is the unstable wind in the E-AGB and thermally pulsing AGB (TP-AGB) phases; the fourth is the stable massive and slow wind with the mass-loss rate higher than 10$^{-7} M_{\odot}$ yr$^{-1}$ and the outflow velocity lower than 20 km s$^{-1}$ in the TP-AGB phase. The mass-loss rates in the first and second types of wind are two or three orders of magnitude lower than the values evaluated by an empirical formula. The presence of massive and slow wind of the fourth type suggests the possibility that the massive outflow observed in TP-AGB stars could be attributed to the Alfv\'en wave-driven wind.


INTRODUCTION
In the post main-sequence phase, the envelope of low and intermediate star expands during the contraction of central core; in the Red Giant Branch (RGB) and the early Asymptotic Giant Branch (E-AGB) phases, the stellar radius reaches typically several tens of the solar radius, and several hundreds of the solar radius in the thermally pulsing AGB (TP-AGB) phase. As stars ascend along the RGB and the AGB, the mass-loss rate becomes larger, which has been revealed by various observations (e.g., Mészáros et al. 2009;Schöier et al. 2013;Ramstedt & Olofsson 2014), and the mass-loss rate in these phases increases with decreasing the surface gravity (Judge & Stencel 1991). The mass loss influences the observable properties such as the spread of the horizontal branch (HB) in the Hertzsprung-Russell (HR) diagram of globular clusters (e.g., Catelan 2000) and the number ratio of TP-AGB to RGB stars (N TP−AGB /N RGB ) on the color magnitude diagram (Rosenfield et al. 2014(Rosenfield et al. , 2016.
Furthermore, Third Dredge-Up (TDU) and Hot Bottom Burning (HBB) experiencing on TP-AGB, whose efficiencies depend on the mass-loss rate, affect the surface abundances of C, N, and O elements, the effective temperature, and the surface gas density (e.g., Marigo 2002;Ventura & Marigo 2010;Tashibu et al. 2017). Consequently, the mass-loss rate influences the chemical composition and amount of dust formed around TP-AGB stars (e.g., Ventura et al. 2012Ventura et al. , 2018 as well as the spectral appearance of the stars. In the stellar evolution calculations, the empirical massloss formulae of Reimers (1975) and Schröder & Cuntz (2005) have been often employed on the RGB and the AGB (e.g., Weiss & Ferguson 2009;Cristallo et al. 2009Cristallo et al. , 2011Marigo et al. 2013). On the other hand, recently Rosenfield et al. (2014Rosenfield et al. ( , 2016 have claimed that an en-hanced mass loss in comparison with the rates evaluated by the empirical formulae is favorable in the pre-dust phase on the TP-AGB to reproduce N TP−AGB /N RGB in nearby galaxies with metallicity ranging from [Fe/H] = -1.59 to -0.56. So far, as the promising physical processes related to mass loss on the RGB or/and AGB, the following three processes have been considered; the dissipation of Alfvén wave which has been believed to act in both phases as long as the chromospheric structure is present; levitation of the atmosphere induced by stellar pulsation; radiation pressure acting on dust grains and the resulting momentum transfer from dust to gas through collisions. The second one alone cannot induce the outflow observed in AGB stars (Schirrmacher et al. 2003). However the hydrodynamical models of pulsation-enhanced dust-driven wind that combines the latter two can induce massive outflow from AGB stars. In the early models, the third one was treated as a simple function of gas temperature (Wood 1979;Bowen 1988). Then the models have been further developed by considering the formation process of dust concretely for carbon-rich (C-rich) AGB stars (e.g., Fleischer et al. 1992;Höfner & Dorfi 1997;Yasuda & Kozasa 2012). Although recent high-resolution observation shows that the pulsation-enhanced dust-driven wind mechanism actually works in evolved AGB stars, the applicability for oxygen-rich (O-rich) AGB stars has yet to be explored (e.g., Ohnaka et al. 2016). Thus, it has been desirable to investigate how the first process affects the mass-loss event on the RGB and the AGB using magnetohydrodynamical (MHD) models (e.g., Vidotto & Jatenco-Pereira 2006).
The mass loss in the RGB phase has been considered to be induced by an outward-directed flux of Alfvén waves in the stellar atmosphere (e.g., Belcher & Olbert 1975;Haisch et al. 1980;Hartmann & MacGregor 1980). As mentioned by Holzer et al. (1983), in the steady-state models, the massive slow winds can be reproduced by adding the energy not only in the subsonic region but also in the supersonic region; the energetics of the Alfvén wave-driven wind needs to be treated as accurately as possible. Furthermore it is unclear whether the stellar winds continuously stream out in the framework of the steady-state models. Then the dynamical MHD simulations for the RGB phase were performed by Suzuki (2007) applying the wind model originally developed for the Sun that reproduces the physical properties of the observed solar corona and wind (Suzuki & Inutsuka 2005. Suzuki (2007) showed that the model for the RGB stars can generate the massive wind whose mass-loss rateṀ varies in time from 10 −10 to 5 × 10 −7 M ⊙ yr −1 However, as pointed out by Airapetian et al. (2010), the wind model does not include magnetic diffusion terms in the induction equation. If the magnetic diffusion is taken into account, the Alfvén waves dissipate more rapidly near the stellar surface, which reduces the energy supply in the outer region. As a result, the wind could be expected to be less dense in realistic situations. One of the main purposes in this paper is to extend the model of Suzuki (2007) by including effects of the magnetic diffusion.
In addition, the evaluation of the radiation field is crucial for the energetics of the stellar wind. In particular, the radiative cooling and heating processes in the region with the temperature lower than ten thousands directly affects the formation of chromosphere in the Alfvén wave-driven wind. However, in the above mentioned dynamical models for the Alfvén wave-driven wind (Suzuki 2007;Airapetian et al. 2010), the radiative heating process has not been treated.
In order to investigate Alfvén wave-driven winds in open magnetic flux tubes piercing the stellar surface of RGB and AGB stars, we extend the MHD model of Suzuki (2007) as follows; first we include the magnetic diffusion terms via Joule dissipation and ambipolar diffusion. Second, evaluating the radiation field based on a modified Unno-Kondo method by Winters et al. (1997), we take into account radiative heating and cooling processes. Then, simulating the stellar winds for model stars whose stellar parameters are taken from the stellar evolution calculations, we examine the dynamical properties of winds along the evolution tracks.
This paper is organized as follows. In Section 2, we describe our wind model, focusing on the parts extended from the model of Suzuki (2007). Then, after the wind types are classified in Section 3.1, Section 3.2 presents the results of the simulated winds on the RGB and the AGB of stars with initial mass M ini = 1.5, 2.0, and 3.0 M ⊙ and initial metallicity Z ini = 0.02 as a representative of solar metallicity. In Section 3.3, the transition of wind type is investigated as a function of stellar parameters. In Section 4.1, we discuss the results of simulations, analyzing the change of wind properties along the evolution tracks of stars. Then we compare the results of the numerical simulations with those based on observations for RGB and E-AGB stars (Section 4.2.1) and for TP-AGB stars (Section 4.2.2). The summary is presented in Section 5.

Magnetohydrodynamical model
In order to construct the Alfvén wave-driven wind model for evolved stars with extended atmospheres, we extend the one-dimensional (1D) MHD simulation code developed for the RGB stars in Suzuki (2007). Here we describe the extensions of the code.
We set an open magnetic flux tube by a filling factor f , and the functional form, following Kopp & Holzer (1976) and Suzuki et al. (2013), is given by where r is the radial distance from the center, f 0 is a filling factor at the stellar surface of the radius r 0 , and h 1 is a typical height of closed magnetic loops (e.g., see Figure 1 of Suzuki et al. 2013). In the simulations, f 0 is set to be 1/1000 as a standard value (e.g., Suzuki et al. 2013), and we adopt h 1 = 0.6 h p with the pressure scale height h p at the stellar surface. Then the radial component of the magnetic field B r is given as by the conservation of magnetic flux, and the radial magnetic field at the stellar surface B r,0 is set to be 1000G as a standard value so that the averaged field strength contributed from the open field regions f 0 B r,0 = 1G.
Using the filling factor f , we solve the equations for the conservation laws of mass, momentum, and energy, and the induction equation that are expressed as where symbols are almost the same as those used in Suzuki (2007); t, ρ, v, p, and B are time, mass density, velocity, thermal pressure, and magnetic field strength, respectively, and subscripts r and ⊥, i denote radial and ith tangential component, respectively; d/dt and ∂/∂t denote Lagrangian and Eulerian derivatives, respectively; ǫ is the specific internal energy and is given by ǫ = p/[(γ − 1)ρ] for an ideal gas with the specific heat ratio γ; G is the gravitational constant; F c = κ 0 T 5/2 (dT /dr) is the thermal conductive flux by Coulomb collisions, which is effective only in high temperature regions with temperature T ≥ 10 4 K, where κ 0 = 10 −6 g cm s −3 K −7/2 ; q R represents the net radiative cooling rate described in Section 2.2; the magnetic diffusivity η is the sum of the Joule diffusion η J and the ambipolar diffusion η AD , each of which is given by and respectively, where c is the the speed of light, e is the elementary charge, and subscripts e, i, and n denote electron 1 , ion, and neutral species, respectively; m e is the mass of electron; ν en is the collision frequency between electron and neutral species and is given as ν en = n n < σ en v en > with the number density of neutral species n n , the electron-ion collision cross section σ en , and the relative velocity between ion and neutral species v en , and where the angle bracket means average over the velocity space. χ is given as where σ in , v in , and m i (m n ) are the ion-neutral collision cross section, the relative velocity between ion and neutral species, and the mass of ion (neutral) species, respectively; x e = n e /n H is the ionization degree that is defined as the ratio of n e to the number density of hydrogen nuclei n H ; We evaluate the numerical values of right-handed sides of equations (8) and (9) by setting < σ en v en > = 8.3 × 10 −10 T /K cm 3 s −1 and < σ in v in > = 1.9 × 10 −9 cm 3 s −1 , referring to equations 14 and 12, respectively, of Draine et al. (1983). Equation (9) is valid if and only if x e ≪ 1. As x e approaches to unity, ρ n becomes lower than ρ, and equation (9) overestimates η AD . Furthermore, in the circumstellar envelope of the model star on the RGB and on the E-AGB, the ionneutral collision time typically exceeds the maximum period of the perturbation waves (set at the stellar surface) in the temperature region of T > 20000 K. Thus we turn off η AD for x e > 0.8. Here we note that the Hall effect does not work in this 1D simulations with-out dependence on the two tangential components. To estimate x e , we disregard the second and higher ionized metal species for simplicity. Then x e is approximately expressed as where n p , n He + , and n He ++ are the number densities of hydrogen ions, first ionized helium ions, second ionized helium ions, respectively; A j is the abundance of jth metal species relative to hydrogen; R j 1c and R j c1 are the photoionization rate and the radiative recombination rate, respectively, between the ground and continuum states for the jth metal species. We derive the ratios n p /n H , n He + /n H , and n He ++ /n H referring to the scheme of Hartmann & Avrett (1984) and Harper (2014, priv. comm.). The ratio R j 1c /R j c1 is set as where k B and ν ı,0 are the Boltzmann constant and the frequency of the photoionization edge. Here we assume that the radiation field of star is approximated by optically thin blackbody radiation characterized by the effective temperature T eff and a geometric dilution factor W defined as where R * is the stellar radius. In equation (11), we also include a galactic ionization term due to the interstellar radiation field (ISRF) particularly in the ultraviolet (UV) region. Referring to the work of Mathis et al. (1983), we set the temperature of the UV component T gal as 7500 K and the dilution factor of the component W gal as 10 −14 .

Radiative cooling and heating
The net radiative cooling term q R appearing in equation (6) is expressed by the difference of the radiative cooling rate q R,cool and the radiative heating rate q R,heat as follows; In this paper we take the following simple treatment for the cooling term q R . We introduce the following three temperature regimes; T < T low , T low ≤ T < T high , and T high ≤ T , where T low = 6000 K and T high = 10000 K. In the high temperature regime with T ≥ T high , we basically adopt optically thin cooling rate given by q R,cool = n p n e Λ, where Λ is the cooling function for the optically thin plasma with solar metallicity (Sutherland & Dopita 1993) 2 . However, this treatment sometimes overestimates the cooling in opticallythick region. In order to avoid this overestimation, we adopt a "cap" taken from the empirical radiative rate = 4.5 × 10 9 ρ erg cm −3 s −1 of the solar chromosphere (Anderson & Athay 1989;Moriyasu et al. 2004) for q R,cool ; q R,cool,high = min(n p n e Λ, 4.5× 10 9 ρ erg cm −2 s −1 ). (14) In the low-temperature range, T < T low , the radiative cooling rate q R,cool is the same as in the constant opacity model for pulsation-enhanced dust-driven wind from TP-AGB stars (e.g., Tashibu et al. 2017). The value in this case q R,cool,low is given by using the Rosseland mean gas opacity at the stellar surface κ surf derived by the stellar evolution calculations, and is simply expressed as where B(T ) and σ are the frequency integrated Planck function and the Stefan-Boltzmann constant, respectively. In the intermediate range, T low ≤ T < T high , we interpolate between q R,cool,low and q R,cool,high in a logarithmic manner: where α = log(T /T low )/ log(T high /T low ). Although equation (16) is a quite simplified prescription, we suppose that it can give a reasonable estimate because the timescale of the radiative cooling is almost logarithmically monotonic in the temperature range of 6000 K < T < 10000 K ( Figure 11 of Woitke et al. 1996). The radiative heating rate q R,heat is simply estimated by using the frequency integrated mean intensity J as follows; 2 On the TP-AGB where the surface abundances of C, N, and O considerably change due to TDU and HBB, this estimate of the cooling rate based on (scaled) solar abundances is problematic in the hot component of chromosphere with T higher than 15000 K. However the chromospheric structure is not formed on this evolutionary phase as shown in Section 3. Therefore we expect that the disregard of the dependency of the cooling rate on the chemical components in the high temperature regime does not affect the simulation results.
This equation is a reasonable approximation in optically thick region unless the scattering is dominant. It should be cautioned that using the Rosseland mean opacity κ surf could underestimate the heating rate in optically thin region, which will be investigated in our future work. We estimate the frequency integrated mean intensity J in equation (17) by using the Unno-Kondo method in the gray approximation (Winters et al. 1997): where R out is the radial position of the outer boundary, and µ r is the cosine of the separation angle, which is derived by integrating from the inner boundary R in with the boundary condition The above framework to estimate radiation field has been adopted in our previous wind models from TP-AGB stars (e.g., Tashibu et al. 2017) with the addition of the dust opacity.

Input parameters
The input parameters of our MHD simulations are the followings: the stellar mass M * , the stellar luminosity L * , the effective temperature T eff as well as the gas opacity κ surf , the average magnetic field f 0 B r,0 , and the fluctuation amplitude δv 0 at the stellar surface.
Among the input parameters, the stellar parameters M * , L * , T eff , and κ surf along the evolutionary track for a model star with a given set of initial mass M ini and initial metallicity Z ini are calculated by the Modules for Experiments in Stellar Astrophysics (MESA) code (version r6208) (Paxton et al. 2011(Paxton et al. , 2013(Paxton et al. , 2015. The method of calculation is the same as in Tashibu et al. (2017), thus we briefly outline the method here; The standard mixing length theory (Cox & Giuli 1968) is applied and convection is approximated as a diffusive process within the convective region defined by Schwarzshild criterion, and the mixing length parameter α MLT is set to be 2.0. The overshooting parameter f is set to be 0.014 at all convective boundaries, except for the bottom of the He-shell flash region at which f is set to be 0.008 throughout the evolution after the first thermal pulse (TP) (Paxton et al. 2011).
In this paper we adopt the CNO-enhanced low temperature opacity constructed by using AESOPUS tool (Marigo & Aringer 2009), and apply the mass-loss formula of Schröder & Cuntz (2005) in the post mainsequence phase. Along the evolutionary track of a model star, the stellar parameters are sampled with the stellar luminosity spacing of 10 0.05 -10 0.2 L * on the RGB and on the E-AGB, and with the stellar mass spacing δM * ≃ 0.05 M ⊙ on the TP-AGB, excluding the periods of thermal pulses.
Following Suzuki (2007), we set f 0 B r,0 = 1G; the value is not so deviated from that estimated for S-type Mira star χ Cyg by Lèbre et al. (2014) (2-3G) and those estimated for three C-rich AGB stars by Duthu et al. (2017) (1.1-9.5G). We note that, among the three AGB stars, the average value of the surface magnetic field strength for IRC+10216 is 3.8 G while those of RW LMi and RY Dra are lower than 4.4 and 4.8G, respectively.
The other parameter, δv 0 , is determined as follows. First we refer to Renzini et al. (1977) for the dependence of the excited acoustic flux F a,0 on the surface gravity g and T eff : where ρ 0 and c S are the mass density and the sound velocity at the stellar surface, respectively. Then, we get the scaling relation of the surface fluctuation, Here the symbol ⊙ represents the Sun, ρ 0,⊙ = 2.07 × 10 −8 g cm −3 , 3 and T eff,⊙ = 5.77 × 10 3 K at the stellar surface where the diluted optical depth is 2/3 in the initial static atmosphere of our MHD model. The value of δv 0,⊙ is set to 1.0 km s −1 (Suzuki 2007). Note that δv 0,scale for the stars on the AGB with R * > 60 R ⊙ could exceed the sound velocity c S except for the case of T eff ≤ 3000 K. Supersonic fluctuation at the surface simply results in the increase of the radiative loss at low altitudes though the structures of the upper atmosphere and envelope are not strongly affected (Suzuki et al. 2013). Therefore, we adopt The wave spectrum of the fluctuation at the stellar atmosphere is set so as to depend on the inverse of frequency according to the recipe described in Section 2.2.3 of Suzuki (2007). Both the longitudinal and transverse fluctuations are included, and those amplitudes are set to be δv 0 .

RESULTS
We present in this section the results of the MHD simulations of the Alfvén wave-driven wind from the RGB and AGB stars whose stellar parameters derived from the stellar evolution calculations for the model stars of M ini = 1.5, 2.0, and 3.0 M ⊙ with Z ini = 0.02, focusing on the wind properties; note that we do not simulate the stellar winds in the HB phase (gray colored regions in Figures 1, 7, and 8) as well as those of stars descending from the tip of RGB (in particular the region between M * =2.9983 and 2.9964 M ⊙ in Figure 8).
The MHD simulations are carried out as follows; Given a set of the input parameters (see Section 2.3), the simulation is continued until a continuous gas-outflow is observed at the outer boundary located at 25 R * and the continuous gas-outflow is referred as the stable wind hereafter for convenience; otherwise we stop the simulation when the simulation time reaches 80 R * /c S and we call the wind observed in this case as the unstable wind hereafter for simplicity. The mass-loss rateṀ and the radial gas velocity v gas that characterize the stable wind are calculated by averaging the values at the outer boundary over the latter half of the simulation time.
In the following subsections, first we summarize the types of winds observed in the simulations (Section 3.1), then the properties of wind observed in the simulations are described for each model star (Section 3.2). The dependence of the transition of wind properties during the evolution of star on the effective temperature T eff as well as the surface gravity g is clarified in Section 3.3.

Classification of wind types
The properties of Alfvén wave-driven wind change during the evolution of star. Roughly speaking, the wind properties are classified into the following four types ( Table 1). The first type is characterized by warm 4 and fast wind with the radial outflow velocity (radial component of the gas velocity) v gas > 80 km s −1 . The wind of the second type exhibits almost time-steady, slow and cool wind with the mass-loss rateṀ < 10 −10 M ⊙ yr −1 . In the third type, steady wind does not stream out but intermittent winds are driven in a sporadic manner. The fourth type is characterized by cool and massive wind withṀ > 10 −7 M ⊙ yr −1 . The first and second types are found in the RGB and E-AGB phases. The third type is seen in the tip of RGB and the AGB phases while the fourth type appears in the late stage of the TP-AGB phase. The detail is provided for M ini = 1.5, 2.0, and 3.0 M ⊙ model stars in the following subsections.

3.2.
Wind properties of the model stars 3.2.1. Mini = 1.5 M⊙ model Figure 1 shows the time evolution of the effective temperature T eff , the stellar radius R * in units of R ⊙ , the mass-loss rateṀ in units of M ⊙ yr −1 , and the radial outflow velocity v gas in units of km s −1 from the RGB to the E-AGB phase; note that the value ofṀ and v gas are time-averaged values at the outer boundary. The blue line in the middle panel shows the mass-loss rate used in the stellar evolution calculation.
Immediately after entering the RGB phase, R * increases from 4.55 to 11.0 R ⊙ ,Ṁ increases from 2.34 × 10 −14 to 1.38 × 10 −13 M ⊙ yr −1 , and v gas somewhat decreases from 539 to 421 km s −1 . Then, as the star evolves, bothṀ and v gas are not so changed until R * becomes 24.0 R ⊙ .
After then, the wind properties drastically change twice on the RGB-phase; first at M * = 1.4609 M ⊙ where R * = 58.49 R ⊙ and second M * = 1.448 M ⊙ where R * = 70.00 R ⊙ . In particular, this is clearly seen in the temporal variation of the outflow velocity shown in the bottom panel. Just before the first change, the winds are tenuous and fast;Ṁ is lower than 1 × 10 −11 M ⊙ yr −1 and v gas is higher than 80 km s −1 , and the wind is classified as the first type (Section 3.1). After the first transition, the chromospheric structures are not always seen, and the wind becomes more massive but quite slow; v gas is lower than 10 km s −1 . This massive wind is classified as the second type (Section 3.1). Then, after the second transition, the winds become further slower, and steady wind does not stream out. We call this unstable wind as the third type. These two transitions between three wind types are also found in E-AGB phase as shown in the bottom panel in Figure 1 (M * = 1.401 and 1.392 M ⊙ ). In the following, we show the radial structures of the three types of the winds. Figure 2 depicts the radial structure of a typical example of the first type wind. We can see from the middle panel that the first wind is a warm wind with the gas temperature T exceeding 10 5 K in the region outside 1.46 R * . The radial gas velocity v r exceeds 200 km s −1 at 6.56 R * , but is still sub-Alfvénic, as can be seen from the bottom panel. As shown in the top panel, ρ is quite low, and does not exceed 10 −20 g cm −3 in the region outside 20 R * . Roughly speaking, the first type wind for the star with R * < 41 R ⊙ is almost in steady state. RGB and E-AGBṀ ≤ 10 −11 M⊙ yr −1 and vgas > 80 km s −1 > 1.1 second type RGB and E-AGB 5 × 10 −11 M⊙ yr −1 ≤Ṁ ≤ 2 × 10 −10 M⊙ yr −1 and vgas < 10 km s −1 0.91 to 1.1 third type tip of RGB and AGB unstable (or sporadic) -1.0 to 1.0 fourth type TP-AGBṀ > 10 −7 M⊙ yr −1 and 2.4 km s −1 < vgas < 15 km s −1 < -0.4 A typical structure of the wind in the second type is shown in Figure 3. As shown in the middle panel, the wind is cool and T is lower than T eff . The wind in the region of r > 9.61 R * moves outward as shown in the bottom panel. The second type wind is not quasi-steady, and it is not clear whether the wind could be sustainable over longer simulation time or not. The stability of the wind of this type will be discussed in Section 4.1.
A structure of the third type wind is presented in Figure 4. The middle panel shows that the chromospheric structure disappears. The outer atmosphere is levitated due to the injection of the Alfvén wave, which can be seen in the top panel. However, as shown in the bottom panel, inflows occur in intershock regions. Therefore this type of wind is considered not to cause any significant mass loss. The behavior of flowing gas in the E-AGB  Figure 5 shows the time evolution of the stellar parameters and the wind properties on the TP-AGB. The Alfvén wave-driven wind again starts to blow continuously after M * ≤ 0.8998 M ⊙ (see the middle and bottom panels), firstly because of the decrease of the surface gravity and secondly because of the suppression of the damping of the magnetic energy associated with the Alfvén waves (see Section 4.1 for the details). We classify this wind as the four type which is characterized as a dense and slow stable wind appearing in the late stage of TP-AGB. The calculated values ofṀ and v gas are 1.7-2.8 × 10 −6 M ⊙ yr −1 and 10-15 km s −1 , respectively.  A snapshot of the radial structure at M * = 0.8998 M ⊙ on the TP-AGB is shown in Figure 6. We can see from the middle panel that the chromospheric structure does not appear. Nevertheless the stable outflow streams out due to the low surface gravity of log g = −0.455. The temperature profile is almost smooth except for the re- gions at r = 2.0, 2.6, and 3.0 R * because the gas temperature is primarily determined by the radiation on the TP-AGB; the temperature of gas heated in the high density flow rapidly falls down to the radiative equilibrium one. Furthermore, in the region outside 4.0 R * , the wind is super-Alfvénic, v r > v A (see the bottom panel), which is in sharp contrast to the first and second types of wind shown in Figures 2 and 3. Thus the physical properties of the fourth type wind is quite different from those of the first or second type winds.
The presence of the fourth type wind clearly demonstrates that Alfvén wave-driven mechanism is promising even for TP-AGB stars not emitting chromospherically active lines. Also, it should be pointed out that the M ini = 1.5 M ⊙ model star is oxygen-rich before the onset of final TP at M * = 0.8539 M ⊙ . Thus the Alfvén wavedriven mechanism should be considered to be one of possible mechanism to trigger the mass loss from O-rich AGB stars for which no self-consistent hydrodynamical model of pulsation-enhanced dust-driven wind is available (see Section 1).  The fourth type wind onsets at M * = 1.000 M ⊙ . In the evolutionary stage after 1.000 M ⊙ , the star is in the C-rich AGB phase. During this phaseṀ = 1.8-2.0 × 10 −6 M ⊙ yr −1 and v gas = 7.9-12 km s −1 , respectively. Although v gas increases with decreasing the stellar mass M * , the mass-loss rateṀ does not increase, in contrast to the model of M ini = 1.5 M ⊙ .

Mini = 3.0 M⊙ model
The evolution of the stellar parameters and the wind properties from the RGB to AGB is depicted in Figure 8. The transition from the first type to the second type occurs at M * = 2.973 M ⊙ (R * = 101.1 R ⊙ ). The wind types changes to the third type at M * = 2.950 M ⊙ (R * = 164.6 R ⊙ ) when the model star enters into the TP-AGB phase, as in the model of M ini = 2.0 M ⊙ .
Then the model star turns to C-rich AGB stars at M * = 2.186 M ⊙ . During the evolution, the effective temperature T eff in the interpulse phase decreases and then turn up. After reaching to the minimum (∼ 2620 K), the fourth type appears at M * = 1.050 M ⊙ withṀ = 5.5 × 10 −7 M ⊙ yr −1 and v gas = 3.4 km s −1 (Figure 8). However, in contrast to the previous two models,Ṁ and v gas gradually decrease as the star evolves;Ṁ = 3.6 (0.95) × 10 −7 M ⊙ yr −1 and v gas = 3.1 (2.4) km s −1 at We can see from the figure that the transition from the first type to second type seems to depend only on the effective temperature; the warm and fast wind (the first type wind) can be seen only in the region of log T eff ≥ 3.66. Also, as can be seen Figure 9, the star with the second type wind characterized by v gas < 10 km s −1 , M < 10 −10 M ⊙ yr −1 , and relatively cool wind (Figures 2  and 3) are located in the very narrow range of 3.646 < log T eff < 3.656. On the other hand, the fourth type (slow and massive) wind can be seen in the late stage of TP-AGB, and the position of transition to the fourth type depends on not only the effective temperature but also the stellar luminosity. The behavior of the fourth type wind seems to depend on the surface gravity as is mentioned in Section 3.2.1. In this section, we shall investigate how the wind properties and the transition of wind type depend on the effective temperature and the surface gravity of stars. Figure 10 shows the dependence of the mass-loss ratė M (top) and outflow velocity v gas (bottom) on the effective temperature T eff for the the first type wind. Roughly speaking, there is a tendency thatṀ increases and v gas decreases with decreasing T eff except for the region where 4800 K <T eff < 5200 K. Figure 11 plotsṀ (top) and v gas (bottom) for all the simulated stars with the stable wind on the RGB and the AGB. The figure clearly shows that the transition from the first to second type occurs at T eff = 4530 K, and the the second type wind disappears at T eff = 4420 K. The stable wind is never driven from the simulated stars whose effective temperature are in the range of 3450 K < T eff < 4420 K, although the atmospheres are lifted up by the Alfvén waves (Figure 4). In T eff < 3450 K , dense, slow, and cool stable winds develop by the dissipation of long wavelength Alfvén waves in the wind acceleration region (see Section 4.1).  Linsky & Haisch (1979). Note that we delete the tracks during the periods of thermal pulse to avoid overlapping the tracks.    Figure 12 plotsṀ (top) and v gas (bottom) versus log g for the first type wind.Ṁ increases and v gas decreases with decreasing log g except for the region of 1.9 < log g < 2.5.

Dependence on surface gravity
As is presented in Figure 13, the wind turns to the second type at log g = 1.098 (1.087, 0.901) for the simulated star of M ini = 1.5 (2.0, 3.0) M ⊙ . No stable wind is driven in the range of -0.456 < log g <0.475 and the value of log g at which the transition to the unstable (the third type) wind occurs slightly depend on M ini ; log g between 0.779 and 1.068, (0.759 and 0.951, 0.475 and 0.761) for M ini = 1.5, (2.0, 3.0) M ⊙ , respectively. Although AGB stars can develop the fourth type wind below log g = -0.456 in the model stars considered in this paper, the value of log g at which the fourth type wind starts to blow significantly depends on M ini . Also the wind properties during the evolution are influenced by not only the surface gravity but also the surface opacity κ surf , as mentioned in Section 3.2.3. The condition for the apperance of the fourth type wind derived from the examination of the results of MHD simulations is expressed as 4. DISCUSSIONS

Analysis of wind properties
Here we examine how the wind types found in the MHD simulations relate with the dissipation of Alfvén waves, using the time averaged radial profiles of physical quantities characterizing the Alfvén wave-driven wind. The radial variation of magnetic perturbation being associated with the dissipation of Alfvén waves directly affect the properties of stellar winds, which is expressed in equation (6). We introduce the relative variation X at r defined as where B ⊥, * ,1 (B ⊥, * ,2 ) is the first (second) tangential component of the magnetic field strength at the stellar surface. The actual tangential components of these quantities at the stellar surface are set to be zero as the boundary condition while the tangential components of gas velocity at the radial position are set as the input parameters as described in Section 2.3. We substitute the values at the stellar surface with the actual values at the second innermost radial computational grid. The value of the magnetic energy density derived from the substituted quantities (B 2 ⊥, * ,1 + B 2 ⊥, * ,2 )/(8π) is in the range of 0.2 to 0.8 × ρ δv 0 2 . Hereafter we investigate the relation between the wind types and the time averaged value of X for the model star of M ini = 1.5 M ⊙ . The radial structures of the time averaged values of X are shown in the top panel in Figure 14 together with the time-averaged radial structures of the gas temperature, the radial velocity, v r , and the Alfvén velocity, v A , in the panels from the second to bottom; Thin black, gray, and thick black lines represent the typical radial structures for the first type wind at M * = 1.4918 M ⊙ , the second type wind at M * = 1.4009 M ⊙ , and the fourth type wind at M * = 0.8998 M ⊙ , respectively.
In the first type wind, X decreases in a exponential manner according to the density structure in the inner nearly static region, and is not so changed in the outer wind region. The value of X decreases by more than five orders of magnitude from the stellar surface to the outer wind region. Even at the radial position where the distance from the star r -R * is several × 10 −2 R * , X decreases by two orders of magnitude and the gas temperature begins to increase due to the dissipation of Alfvén waves. Then the gas velocity starts to increase steeply around r -R * is 0.1-0.2 R * . The increase of the gas temperature decreases the relative abundances of neutral species, which suppresses the ambipolar diffusion and the further decrease of X in the outer wind part.
After the effective temperature T eff drops down below 4900 K, as the star evolves, the wind speed of the first type wind decreases the with stellar evolution. This seems to be ascribed to the decrease of the thermal enthalpy due to the more efficient radiative cooling in denser wind. Then the winds become denser and slower (second type) as long as a part of the dissipated energy is converted to the kinetic energy of gas necessary to maintain the stable outflow. Considerably later the stable wind in the fourth type is formed. As shown in the second panel in Figure 14, the temperature at the wind base in this fourth type is quite lower than that for the first type. The determination of the stability of the various types of winds requires the analysis of the energetics in the outer wind region.
We analyze the energetics in the outer wind region by introducing the quantities Y i (i=1,2,3) defined as where Y 1 , Y 2 , and Y 3 are kinetic energy density, thermal enthalpy, and the magnetic energy density associated with the Alfvén waves normalized by the absolute value of the product of the gas density and the gravitational potential energy, respectively. We derive these quantities from the time averages of numerators and denominators on the right-hand sides of equations (26), , and (28). The radial profiles of Y -s are given in Figure 15. Here we note that unnormalized energy densities themselves decrease with r.
The top panel of Figure 15 shows the results at M * = 1.4918 M ⊙ . In this first type wind, the three kinds of energy are predominant over the gravitational energy in the outer wind part, and the wind is warm and fast. The predominance of thermal enthalpy leads to the suppression of the ambipolar diffusion, which results in the predominance of the perturbation energy.
In the second type wind at M * = 1.4009 M ⊙ (the middle panel of Figure 15), Y 3 is much less than that in the first type wind, and is 0.4-0.6 in the region of r > 2.2 R * . This is due to the damping of the magnetic perturbation in the outer wind region as shown in the top panel of Figure 14. Y 2 is less than 0.1 through the outer wind region because of the more efficient radiative cooling in denser wind. Also, in this type, the gas temperature is not strongly affected by the dissipation of Alfvén waves; as shown in Figure 3, the gas temperature does not deviate from the radiative equilibrium one so much. Therefore, although the perturbation energy is not so quite lower than the gravitational energy, the wind is slow due to the inefficient energy conversion from the magnetic energy to the kinetic energy. The value of Y 1 never exceeds unity even at the outer boundary (25R * ); This casts doubt on the stability of the second type wind. It is considered that the sporadic event such as the propagation of shock barely enables the second type wind to stream out continuously.
The results at M * = 0.8998 M ⊙ for the fourth type wind is presented in the bottom panel of Figure 15. We can see that Y 3 is 0.07-0.4 in the region of r > 2.0 R * ; it reaches the maximum value (0.4) at r = 2.3 R * and decrease with r. This attributes to the further suppression of the damping of the magnetic perturbation in the inner atmosphere (see top panel of Figure 14), and the efficient energy conversion in the wind part. In the region of r > 2.3 R * , Y 2 is larger than 0.1, and is larger than that in the second type wind, which reflects the considerable decrease of the surface gravity arising from the increase of R * and the decrease of M * although T eff substantially decreases. Also, as in the second type, the gas temperature is not strongly affected by the dissipation of Alfvén waves as shown in Figure 6. However the ratio of the thermal enthalpy to the perturbation energy (Y 2 /Y 3 ) is higher, and the energy conversion from the perturbation to kinetic energy is more efficient than in the second type. The kinetic energy increases with increasing the distance, while the perturbation energy decreases. At r = 6.7 R * , Y 1 reaches 1.0, and the stable super-Alfvénic wind is formed.
The results of analysis are summarized as follows: The transition from the first to second type occurs mainly by the suppression of the thermal enthalpy at the base of the wind. The second type wind could not be always stable. The appearance of the fourth type wind seems to be related with the degree of the suppression of the damping of the magnetic perturbation in the inner region; a larger fraction of the input Alfvénic Poynting flux remains at a higher altitude to derive dense wind.

Comparison with observations 4.2.1. Comparison with RGB and E-AGB stars
The mass-loss rateṀ of the first type wind seen in the simulated RGB and E-AGB stars with T eff (log g) higher than 4530 K (1.1) is at least two orders of magnitude smaller than the rates from the empirical formula of Schröder & Cuntz (2005) (see middle panels of Figures 1, 7 and 8); The calculatedṀ is also smaller than the rates of metal poor red giants with high T eff evaluated from the optical/near infrared observations by Dupree et al. (2009) andMészáros et al. (2009), despite that the mass-loss rate is considered to be higher in metal-rich star than metal-poor star Mészáros et al. (2009) or little dependent of metallicity McDonald & Zijlstra (2015).
The outflow velocity v gas of the simulated stars with high T eff ( 4530 K) exceeds 200 km s −1 . The value of v gas is hard to estimate from observations of chromospheric lines while the observations of He I 10830Å line by Dupree et al. (2009) may provide the radial velocity of the wind with gas temperature higher than 20000 K. However, for 13 giants with high T eff in their samples, v gas is at most 104 km s −1 for HD121135. The difference of v gas may reflect, at least in part, that their estimated values do not necessarily represent the wind velocities in outer low density region; As shown in the bottom panel of Figure 2, the outflowing gas increases with increasing the distance from the center, even at ∼ 10 R * , and it is not clear whether the higher velocity component contributes to the blue wing of observed line profiles or not.
The outflow velocity v gas drops down lower than 10 km s −1 at T eff (log g) around 4500 K (1.1) in our simulations, being accompanied with the transition from the first to second type. However, the UV observations of chromospheric lines clearly show that the wind velocities slow down lower than 20 km s −1 at T eff = 3500 K more gradually, as is depicted in Figure 11 by dotted line (see equation (2) of Wood et al. 2016).
Thus, the problematic discrepancies between the results of MHD simulations and the observations are the following two. (1) the mass-loss rate for the stars with T eff 4530 K. (2) the outflow velocity for the stars with T eff in the range of 3450 to 4420 K.
The first discrepancy may arise from low values of the input parameters f 0 B r,0 and δv 0 ; f 0 B r,0 being increased from 1 to 5G,Ṁ increases by an order of magnitude or more in some cases. Some active giants with high T eff at the base of the RGB (EK Eri, V390 Aur, EI CnC, and ι Cap) have further strong surface magnetic field (Aurière et al. 2015). We will investigate systematically the dependence of the surface magnetic field in future work. For the weakly-active magnetic star at the base of RGB such as Pollux whose surface averaged longitudinal magnitude is sub-Gauss (Aurière et al. 2009), larger δv 0 may be required; recently O'Gorman et al. (2017) estimates an upper limit ofṀ of Pollux;Ṁ to be 3.7 × 10 −11 M ⊙ yr −1 , which is higher than those of the first type wind in our model with f 0 B r,0 = 1G.
On the other hand, the second discrepancy is not improved by adjusting the input parameters; When f 0 B r,0 increases from 1 to 5G, the second type wind changes to the third type wind due to the more efficient cooling in denser wind. In order to maintain the chromospheric structure and the continuous gas-outflow for the stars with T eff in the range of 3450 to 4420 K the formation of the local low density region(s) with high tempeature could be required, which will be able to be investigated in higher dimensional simulations. Also in this context, the applicability of the Rosseland mean opacity at the stellar surface should be examined since the radiative cooling rate evaluated in the present paper may be overestimated.
However, it is difficult to clarify whether the observed wind properties can be reproduced by the MHD simulation, since the evolutionary stage as well as the stellar parameters is not well constrained from the observations, in particular for dust-enshrouded AGB stars. For example, the effective temperature T eff is significantly different, depending on the method of observation and analysis; T eff of R Leo is set to be 1800 (2000) K in the radiative transfer models of Schöier et al. (2013) and Ramstedt & Olofsson (2014) while Perrin et al. (2004) have estimated T eff to be 2980 K based on the interferometric observations.
Here we examine the properties of the fourth type wind from the well-studied two TP-AGB stars; O-rich AGB star o Cet and C-rich AGB star IRC+10216. In Table 2, the input parameters of the MHD simulations for the two TP-AGB stars are provided as well as the wind properties derived from the observations and the simulations. Among the input parameters, M * , L * , and T eff are taken from the literatures (see the caption in Table 2); the data are evaluated from the observations and the models. Since the surface elemental composition is not well-known for both stars, we set κ surf to be 2.0 × 10 −4 cm 2 g −1 which is adopted in the hydrodynamical models of pulsation-enhanced dust-driven wind as a fiducial value (e.g., Bowen 1988;Fleischer et al. 1992;Yasuda & Kozasa 2012).
The observationally unknown value of f 0 B r,0 for o Cet is set to be 4.1G, by referring to the optimized value in the steady-state model of Thirumalai & Heyl (2013). Then, as shown in the second column of Table 2, the cool and dense wind (the fourth type wind) can be generated by the Alfvén wave-driven mechanism. Thus this result Table 2. Observed wind properties (mass-loss rateṀ obs and outflow velocity v gas,obs ), input parameters of our wind models (M * , L * , T eff , κ surf , f0Br,0, and δv0 ) and the wind properties derived from our model (mass-loss rateṀ model and outflow velocity v gas,model ) for two TP-AGB stars. References.
The observed slow and moderate density wind from o Cet may be reproduced by decreasing f 0 B r,0 and δv 0 .
Decreasing the values, we confirm that the wind is still stable as long asṀ is higher than 10 −6 M ⊙ yr −1 . However when f 0 B r,0 ( δv 0 ) is decreased to 1.5G (2.0 km s −1 ), the wind becomes unstable: the observed stable wind withṀ around 10 −7 M ⊙ yr −1 is not reproduced in our model for o Cet.
Here we shall raise the following three possibilities to explain the discrepancy between the observation and our simulation for o Cet. (1) lower surface opacity, since larger opacity makes the wind slower and less dense, as is seen in the model star of M ini = 3.0 M ⊙ (see Section 3.2.3).
(2) the wind from o Cet is actually sporadic rather than stable; slow and less dense wind tends to be unstable not only in our MHD model but also in the pulsation-enhanced dust-driven wind model including the effect of drift of dust (Sandin 2008;Yasuda et al. 2018, in preparation), and this possibility will be explored from the time variation of the radial velocity of H 2 O masers (e.g., Sudou et al. 2017 for R Crt). (3) a hybrid mechanism including dust-driven wind mechanism which may stabilize the unstable wind (the third type wind) in our MHD model with lower f 0 B r,0 or lower δv 0 . The investigations of these possibilities are left for the future works.
The value of f 0 B r,0 for IRC+10216 is set to be 3.8G, by referring to Duthu et al. (2017). Among the stellar parameters, M * is set to 0.8 M ⊙ which is the intermediate value of the current core mass (0.7-0.9 M ⊙ ) estimated by Ladjal et al. (2010) since the core mass is actually lower expected from the adopted lower value of L * . As can be seen from Table 2, the massive quasi-steady wind from IRC+10216 can be driven by the Alfvén wavedriven mechanism, although v gas in our model (11.4km s −1 ) is somewhat lower than that derived from the observations (e.g., 14.6 km s −1 from Knapp et al. (1998)).
The following should be noted; The slow and massive winds from C-rich AGB stars have been investigated in the framework of pulsation-enhanced dust-driven wind models based on the formulation by Gail & Sedlmayr (1988) and Gauger et al. (1990). The hydrodynamical models roughly reproduce the observed dynamical behavior of these stars withṀ ≥ 10 −6 M ⊙ yr −1 (Winters et al. 1994(Winters et al. , 1997Nowotny et al. 2005). The wind properties derived from our MHD simulations should be modified by taking into account the force acting on carbon grains. Furthermore it is a fact that, at present time, no satisfactory hydrodynamical model is available for reproducing the observed physical properties of dusty wind from O-rich AGB stars (e.g., Ohnaka et al. 2016). Thus the wind properties of dusty TP-AGB stars should be investigated by using our MHD model coupled with dust formation calculation.

SUMMARY
The properties and the stability of the Alfvén wavedriven winds from the RGB and the AGB stars are investigated by employing the MHD model extended from the model developed by Suzuki (2007): the Joule resistivity and the ambipolar diffusion are included for the magnetic diffusion terms, the radiative cooling rate at low temperature is modified, and the radiative heating rate is estimated based on the radiation field derived using the Unno-Kondo method of Winters et al. (1997). The model is applied to investigate the stability and the properties of wind from the RGB and the AGB stars whose stellar parameters necessary for the MHD simulations are calculated by the MESA code for the stars with the initial mass M ini in the range of 1.5 to 3 M ⊙ and the initial metallicity of Z ini = 0.02. The results of the simulations are summarized as follows: 1. The stars with T eff > 4530 K in the RGB and E-AGB phases exhibit the chromospheric structures and stable winds, regardless of the initial mass. The stable wind is characterized by the high gas temperature exceeding 10 5 K, the high outflow velocity of v gas > 80 km s −1 , and the low mass-loss rate ofṀ ≤ 10 −11 M ⊙ yr −1 , and is classified as the first type wind. The mass-loss rate is at least two orders of magnitude smaller than that of Schröder & Cuntz (2005), which may indicate that the empirical formulae are derived based on the observations of the magnetically active giants with f 0 B r,0 > 5G.
2. In the stars with 4420 K < T eff ≤ 4530 K on the RGB and E-AGB phases, chromospheric structures are not always seen. The winds blow incessantly, but characterized by the low gas temperature lower than T eff , the quite low velocity of v gas < 10 km s −1 and the mass-loss rateṀ in the range of 5 × 10 −11 -2 × 10 −10 M ⊙ yr −1 . The wind classified as the second type seems not to be stable.
3. In the regime of 3450 K < T eff ≤ 4420 K neither stable wind nor chromospheric structure forms in the M ini = 1.5 M ⊙ , while the atmosphere is lifted due to the injection of Alfvén wave (the third type wind). The levitation of atmosphere is typically seen in relatively swollen stars on the E-AGB and stars on the TP-AGB, as well at M * = 1.4267 M ⊙ around the tip of RGB. In these stars, other mechanisms are required to keep stable wind blowing.
4. On the TP-AGB of all model stars, the formation of stable Alfvén wave-driven referred to as the fourth type is possible because the damping of magnetic perturbation is more suppressed in the inner atmosphere and the gravity in the outer wind part gets lower with evolving the star. The stable stellar winds are slow, dense, and super-Alfvénic without any chromospheric structure.
The dissipation of Alfvén waves injected from the stellar surface can lead to stable winds on the stars in the RGB and AGB phases which are plotted in HR diagram in Figure 9. In particular, it is shown for the first time that the slow and massive wind can be magnetically driven without any chromospheric structure in low-gravity AGB stars. In the future work we will investigate systematically the dependence of the surface magnetic field strength on the wind properties along the stellar evolution. Furthermore we will extend our code to investigate how stellar pulsation and dust formation play a role in driving the wind in AGB stars.