Pulsation-induced atmospheric dynamics in M-type AGB stars. Effects on wind properties, photometric variations and near-IR CO line profiles

TL;DR

This study investigates how different pulsation-driven boundary conditions in dynamical models of M-type AGB stars affect atmospheric structure, wind properties, and observable features, aligning models more closely with observations.

Contribution

It introduces varied boundary conditions in dynamical models to better replicate observed stellar wind velocities, mass-loss rates, and photometric and radial velocity variations in Mira variables.

Findings

Models with adjusted boundary conditions produce realistic wind velocities.

Photometric variations in models match observed Mira star variations.

A phase shift of 0.2 between luminosity and radial velocity maxima is necessary for accurate radial velocity curves.

Abstract

Wind-driving in asymptotic giant branch (AGB) stars is commonly attributed to a two-step process. First, matter in the stellar atmosphere is levitated by shock waves, induced by stellar pulsation, and second, this matter is accelerated by radiation pressure on dust, resulting in a wind. In dynamical atmosphere and wind models the effects of the stellar pulsation are often simulated by a simplistic prescription at the inner boundary. We test a sample of dynamical models for M-type AGB stars, for which we kept the stellar parameters fixed to values characteristic of a typical Mira variable but varied the inner boundary condition. The aim was to evaluate the effect on the resulting atmosphere structure and wind properties. The results of the models are compared to observed mass-loss rates and wind velocities, photometry, and radial velocity curves, and to results from 1D radial pulsation models. Dynamical atmosphere models are calculated, using the DARWIN code for different combinations of photospheric velocities and luminosity variations. The inner boundary is changed by introducing an offset between maximum expansion of the stellar surface and the luminosity and/or by using an asymmetric shape for the luminosity variation. Models that resulted in realistic wind velocities and mass-loss rates, when compared to observations, also produced realistic photometric variations. For the models to also reproduce the characteristic radial velocity curve present in Mira stars (derived from CO $\Delta v = 3$ lines), an overall phase shift of 0.2 between the maxima of the luminosity and radial variation had to be introduced. We find that a group of models with different boundary conditions (29 models, including the model with standard boundary conditions) results in realistic velocities and mass-loss rates, and in photometric variations.

Tables

Table 1. Table 1: Parameters of the dynamical model (adopted from Nowotny et al., 2010 ) used in this study.

Model:	M2
$L_{\star }$ [L⊙]	7000
$M_{\star }$ [M⊙]	1.5
$T_{\star }$ [K]	2600
$R_{\star }$ [R⊙]	412
[Fe/H]	0
C/O [by number]	0.48
Period [days]	490
$\Delta u_p$ [kms-1]	6
$f_L$	1.5

Table 2. Table 2: Polynomial constants, for Eq. 8 .

m	$a_m$	$\Delta {\varphi }_s^m$
0	-5.3309e+04	$\Delta {\varphi }_s^0$
1	3.6367e+04	$\Delta {\varphi }_s^1$
2	-9.7661e+03	$\Delta {\varphi }_s^2$
3	1.3385e+03	$\Delta {\varphi }_s^3$
4	-1.0460e+02	$\Delta {\varphi }_s^4$
5	5.5368e+00	$\Delta {\varphi }_s^5$

Table 3. Table 3: Number of terms used in this work, for different values of Δ ​ ϕ s Δ subscript italic-ϕ 𝑠 \Delta\phi_{s} .

$\Delta {\varphi }_s$	$w$	N where
0.05	2.66	13
0.1	1.81	37
0.15	1.36	115
0.2	1.16	260

Equations

$$\beta =\frac{\dot{M}{v}_{\infty }c}{L_{\star }}$$

$$\Delta {\varphi }_{\mathrm{v}}=\frac{\Delta max+\Delta min+\Delta zero}{3}.$$

$$R_{in}(t)=R_0+\frac{\Delta u_{p}P}{2\pi }\sin \left(\frac{2\pi }{P}t\right)$$

$$u_{in}(t)=\Delta u_p\cos \left(\frac{2\pi }{P}t\right)$$

$$\Delta L_{in}(t)$$

$$=f_L\left({\left[1+\frac{\Delta u_{p}P}{R_{0}2\pi }\sin \left(\frac{2\pi }{P}t\right)\right]}^2-1\right)\times L_0$$

$$\Delta L_{in}(t,\Delta {\varphi }_p)=f_L\left({\left[1+\frac{\Delta u_{p}P}{2\pi R_0}\sin \left(2\pi \left(\frac{t}{P}+\Delta {\varphi }_p\right)\right)\right]}^2-1\right)\times L_0$$

$$\Delta L_{in}(t,w(\Delta {\varphi }_s),\Delta {\varphi }_p)=k\sum _{n=1}^N\frac{1}{n^{w(\Delta {\varphi }_s)}}\sin \left(\frac{2\pi n}{P}\left(t+\Delta {\varphi }_p\right)\right)$$

$$w(\Delta {\varphi }_s)\approx \sum _{m=0}^M{a}_m\times \Delta {\varphi }_s^m$$

$$L2=\sqrt{\frac{1}{x_{steps}}\sum _i^{x_{steps}}{\left(f_i-g_i\right)}^2}$$

Full text

*Aims. *We test a sample of dynamical models for M-type AGB stars, for which we kept the stellar parameters fixed to values characteristic of a typical Mira variable but varied the inner boundary condition. The aim was to evaluate the effect on the resulting atmosphere structure and wind properties. The results of the models are compared to observed mass-loss rates and wind velocities, photometry, and radial velocity curves, and to results from 1D radial pulsation models. The goal is to find boundary conditions which give realistic atmosphere and wind properties.

*Methods. *Dynamical atmosphere models are calculated, using the DARWIN code for different combinations of photospheric velocities and luminosity variations. The inner boundary is changed by introducing an offset between maximum expansion of the stellar surface and the luminosity and/or by using an asymmetric shape for the luminosity variation. Ninety-nine different combinations of theses two changes are tested.

*Results. *The model atmospheres are very sensitive to the inner boundary. Models that resulted in realistic wind velocities and mass-loss rates, when compared to observations, also produced realistic photometric variations. For the models to also reproduce the characteristic radial velocity curve present in Mira stars (derived from CO $\Delta v=3$ lines), an overall phase shift of 0.2 between the maxima of the luminosity and radial variation had to be introduced. This is a larger phase shift than is found by 1D radial pulsation models.

*Conclusions. *We find that a group of models with different boundary conditions (29 models, including the model with standard boundary conditions) results in realistic velocities and mass-loss rates, and in photometric variations. To achieve the correct line splitting time variation a phase shift is needed.

1 Introduction

The mass loss of asymptotic giant branch (AGB) stars through a slow stellar wind is presumably caused by a two-step process. First, stellar pulsations create shock waves in the surface layers of the star that levitate matter in the atmosphere. This compressed, levitated material then reaches temperatures that are cool enough for dust condensation to occur. In this second stage the dust is accelerated outwards by radiation, through scattering or absorption, depending on the chemical composition and size of the dust grains. A stellar wind is triggered, as momentum is transferred from the dust to the gas through collisions.

This scenario has been investigated extensively, and is supported by various observations. The gas dynamics in the region where shock waves develop are studied through high-resolution spectroscopy, e.g. vibration-rotation lines of CO. The CO molecule is stable throughout the atmosphere and the wind acceleration region, making CO lines suitable to probe the inner regions, with regular pulsation, the shock waves and the steady outflow (see e.g. Hinkle et al., 1982; Nowotny et al., 2005a, b, 2010). Interferometry and high angular resolution imagining are also becoming important for studying both the structure and the dynamics of AGB star atmospheres, giving new insights into the shock propagation and dust condensation distances (see e.g. Chandler et al., 2007; Karovicova et al., 2013; Ohnaka et al., 2012, 2016).

The wind driving process is usually studied with dynamical models (see Höfner, 2015, for a recent review of these models) which simulate the stellar atmosphere where radiative effects dominate. This is a vastly different region from the stellar interior, where convection is the main energy transport mechanism and where the stellar pulsations originate. Dynamical wind models typically have an inner boundary situated just below the stellar photosphere, and therefore do not include any description of either convection or pulsation driving. However, the variation of the upper layers of the star is vital to wind driving as this induces the shock waves that facilitate dust formation. A parameterised prescription is typically used in the dynamical wind models, to simulate the temporal variation of the stellar luminosity and the radial velocity of the gas layers at the inner boundary. This prescription has historically had a simple sinusoidal form, based on the initial efforts to describe radial variation in AGB stars by Bowen (1988), and has the advantages of few free parameters. However, self-excited interior pulsation models (e.g. Ireland et al., 2011) and observations (e.g. Nowotny et al., 2010; Lebzelter, 2011) both indicate that this approach should be improved.

It has further been shown that assumptions made about the inner boundary conditions may have strong implications for the structure of the resulting model atmosphere, and consequently for the mass-loss rate and wind velocity, at least in the case of C-type AGB stars (Liljegren et al., 2016).

In this work we focus instead on self-consistent time-dependent models for dynamical atmospheres and dust-driven winds of M-type AGB stars. The models by Höfner (2008) and Bladh et al. (2013, 2015) are based on a wind acceleration scenario by photon scattering on large, near-transparent silicate grains, which is supported by recent observational studies (see e.g. Norris et al., 2012; Ohnaka et al., 2012, 2016). Here we test a range of different inner boundary conditions for M-star atmospheric models. As the silicate dust present in the atmosphere of M-stars are more transparent than the amorphous carbon in C-stars, M-stars are generally less obscured by dust. There is therefore more observational material available for M-stars, e.g. high-resolution spectra, which facilitates comparisons between model results and observations.

To approach realistic pulsation properties and to pin down free parameters in the DARWIN models, we investigate the systematic effects of different pulsation properties (phase shifts between radial and luminosity variations, and different shapes of luminosity variations) on dynamical wind models. The resulting models are also compared to various observables. A set of 99 models with the same stellar parameters, but with a combination of different inner boundaries (see Sect. 2.1.1), is calculated and then evaluated using three criteria:

• •

Wind velocity and mass-loss rate - A direct output from dynamical wind models (DARWIN models, Sect. 2.1), and a measure of the general dynamics. The models are evaluated by comparing the velocity and mass-loss rate combinations with observations from Olofsson et al. (2002) and González Delgado et al. (2003) in Sect. 4.1.

• •

Photometric variations - The brightness in individual filters will vary during a pulsation cycle for AGB stars, due to various absorbers occurring throughout the atmosphere. This creates observable loops in colour-colour diagrams. The shape and positions of these loops will depend on the overall atmospheric structure and dynamics. Synthetic $(J-K)$ vs. $(V-K)$ loops are calculated for the DARWIN models, using the COMA code (Sect. 2.2.1), and are compared to observations in Sect. 4.2.

• •

High-resolution line profiles - Vibration-rotation CO lines are formed in different layers of the atmospheres and wind acceleration regions. The Doppler-shifted CO second overtone line profiles and derived radial velocity (RV) curves reflect the velocity field of the pulsating photospheric layers of the star. Synthetic spectra for this line are calculated using the model atmospheres with COMA (Sect. 2.2.2) and resulting lines and RV curves are compared to observations in Sect. 4.3.

These tests are indicators of different aspects of the AGB atmosphere, and realistic models should reproduce the observations of all the three aspects at once. To the authors’ knowledge this is the first study of its kind on M-type AGB star wind models.

The possibility of using the tests as diagnostics for pulsation models is also explored. The results from the three tests are compared to pulsation properties derived from 1D self-excited radial pulsation models by Bessell et al. (1996), who provided a mathematical description in the form of Fourier series of the luminosity variation and the sub-photospheric dynamics from their models. This can be directly compared to the boundary condition used in the DARWIN models. Furthermore, these pulsation models and follow-up simulations have been used extensively for the interpretation of various observations over the years, despite the lack of dust-driven winds in the models of Bessell et al. (1996).

2 Modelling methods

2.1 Dynamical wind models - DARWIN

Pulsation, shock, and dust formation processes occurring in the atmospheres of AGB stars are inherently time dependent. To investigate the influence of different pulsation properties, atmosphere models are calculated using the Dynamical Atmosphere and Radiation-driven Wind models with Implicit Numerics (DARWIN) code. The DARWIN models simulate the time-dependent structures of an AGB star atmosphere using 1D frequency-dependent radiation-hydrodynamics, including a detailed treatment of dust growth and evaporation (for description see Höfner et al. 2016 and references therein).

For M-type AGB stars the winds are driven by photon scattering on large silicate (Mg2SiO4) dust grains, which grow through reactions with Mg, SiO and H2O. The treatment of dust growth is time dependent, following the method described by Gail & Sedlmayr (1999), with the growth rate limited by the available Mg atoms. As the nucleation of dust grains in an oxygen-rich environment is not well understood (see Gail et al., 2016; Gobrecht et al., 2016, and references therein), pre-existing seed particles consisting of 1000 monomers are introduced. The seed particle abundance is a free parameter and is defined as the ratio between number density of initial grains and number density hydrogen atoms. This is set to $3\times 10^{-15}$ for all models calculated in this paper, following the findings in Bladh et al. (2015).

The starting point for the simulations is a hydrostatic dust-free atmospheric structure whose fundamental parameters are effective temperature, luminosity, and chemical composition. The variations of the luminosity and the gas velocity at the inner boundary are gradually ramped up to simulate the pulsation of the star, transforming the hydrostatic initial atmosphere into a dynamical model. Using an adaptive spatial grid, the full system of non-linear partial differential equations describing gas, dust, and radiation is solved simultaneously with a Newton-Raphson scheme.

2.1.1 Inner boundary conditions

The spatial range of the DARWIN models, for models that develop a wind, reaches from just below the photosphere ($\sim 0.9R_{\star}$) out to where the wind has reached its terminal velocity $v_{\infty}$ ($\sim 20-30R_{\star}$). Since the inner boundary is located above the pulsation driving region, the variability of the stellar radius and the luminosity present in AGB stars is simulated with an ad hoc temporal variation of the physical quantities at the inner boundary. The form of these temporal variations has historically been very simple; $\Delta R_{in}\propto\sin(2\pi t/P)$ for the radial variation and $\Delta L_{in}\propto\Delta R_{in}^{2}$ for the luminosity variation. Such a sinusoidal description, where the expansion and contraction of the stellar surface is locked in phase with luminosity variation, has been used in several different variations and generations of dynamical atmosphere and wind models (e.g. Bowen, 1988; Winters et al., 2000; Höfner et al., 2003, 2016).

This simplified way of describing the time evolution of the luminosity at the inner boundary is not, however, representative of what is known about AGB stars pulsation properties. The shape of the luminosity variation is known to differ: Lebzelter (2011) found that around 30$\%$ of Mira light curves deviated significantly from a sinusoidal curve and presented different shapes and secondary maxima. Interior pulsation models, by e.g. Bessell et al. (1996) or Ireland et al. (2008, 2011), predict phase shifts between the variation of the surface layers and the luminosity. In Nowotny et al. (2010) atmosphere models of C-stars with the standard boundary condition were calculated; however, a synthetic phase shift had to be added to match the models to observations. To arrive at a better approximation of what a realistic boundary condition for the DARWIN models should be, we perform a systematic investigation of the influence of the inner boundary condition on various model properties, and compare the results with available observations.

The inner luminosity boundary is modified in two ways: i) by introducing a phase offset $\Delta\phi_{p}$ between the radial variation and luminosity variation and ii) by changing the shape of the luminosity variation, making it increasingly asymmetric (measured by $\Delta\phi_{s}$), as seen in panels 2 and 3 in Fig. 1. The result of both these changes is a phase shift between the maximum expansion and maximum luminosity of the star. The notation to describe these phase shifts introduced in Liljegren et al. (2016) is used here: $\Delta\phi_{p}$ for case i) (panel 2 in Fig. 1) and $\Delta\phi_{s}$ for the asymmetric case ii) (panel 3 in Fig. 1). The parameter $\Delta\phi_{s}$ is analogous to $\Delta\phi_{p}$, and measures the difference between the maximum stellar expansion and the maximum of the luminosity, which in case ii) is shifted because the luminosity variation is asymmetric. A larger $\Delta\phi_{s}$ then indicates an increasingly asymmetric light variation. The combination of these two changes will lead to a total phase shift between the maximum expansion and maximum luminosity as $\Delta\phi_{tot}=\Delta\phi_{p}+\Delta\phi_{s}$ (see panel 4 in Fig. 1).

For a full mathematical description of the boundary condition, see Appendix A.

2.1.2 Model set

The model parameters used for the starting model (see Table 1) are based on stellar parameters of Model M in Nowotny et al. (2010); the only difference is in the chemical composition (in previous work C/O$=1.4$, here C/O$=0.48$). The atmospheric model here therefore has an O-rich chemistry, in contrast to the C-rich chemistry in the atmospheric model used in Nowotny et al. (2010). Continuing the naming tradition we call our model M2. This combination of parameters was found to reproduce very well the discontinuous S-shaped radial velocity (RV) curve, characteristic of Mira stars. The RV curve is derived from high-resolution time series spectra containing lines formed in regions with strong shocks (around $R\approx 1R_{\star}$), and can be used to deduce the dynamics of the inner layers of Mira atmospheres (e.g. Sect. 5.2 of Nowotny et al., 2010).

Lebzelter & Hinkle (2002) found that RV curves of a sample of Miras, covering a range of periods between 145 and 470 days, showed a universal behaviour with very similar velocity amplitudes (difference between minimum and maximum velocities) and discontinuities in the RV curve around maximum bolometric phases. They also showed that while properties such as mass-loss rate and wind velocity differ between different Mira stars, the atmospheric dynamics of the inner layers seem to be similar, independent of metallicity, spectral type, period and chemistry. It is therefore probably a good assumption that the result for model M2 can be generalised to other combinations of stellar parameters for Mira stars.

A set of 99 models with different inner boundary conditions for Model M2 is investigated. An overview of the set of models can be seen in Fig. 2. As observations (see e.g. Nowotny et al., 2010) and models (see e.g. Ireland et al., 2011) both predict a positive phase shift, i.e. the stellar surface variation lags the luminosity variation, we constrain the investigation to cases where $\Delta\phi_{tot}\in[0,0.5]$.

Of the 99 models, around 20$\%$ did not converge, in most cases due to the time step becoming very small. Such crashes are caused by extreme conditions (strong temporal changes) in models where the radial variation and luminosity variation are very much out of phase. Fortunately, the subsequent analysis showed that these models are far from the realistic region in parameter space, and will therefore not be discussed further.

2.2 Spectral synthesis

Photometry and high-resolution spectra are calculated for comparison with observational data. These calculations are done by applying an a posteriori approach, using resulting atmospheric structures from the DARWIN runs. For every model twenty snapshots per pulsation cycle, equidistant in time, are chosen. For each of these snapshots opacities are calculated using the COMA code, with the assumptions of LTE and a microturbulence of $\xi_{t}=$2.5 kms*-1*. Radiative transfer is then solved for the derived opacities. For a detailed description of the COMA code see Aringer et al. (2016) and for a description of how the line profiles are calculated see Nowotny et al. (2005a, b, 2010).

2.2.1 Photometry

Comparisons using photometric colours are of interest as they trace variations of spectral energy distributions during a pulsation cycle. Large variations in the luminosity create strong variations in gas temperature, which in turn influence the molecular abundances in the atmosphere. These variations are especially prominent in the V-band, which is dominated by TiO, while H2O affects the near-infrared wavelength region. Here the approach of Bladh et al. (2013, 2015) is used, investigating the colour-colour diagram $(J-K)$ vs. $(V-K)$. Typically, observed M-stars show large variations in $(V-K)$ during a cycle, caused by changes in the molecule abundances (mainly TiO), while the variations in $(J-K)$ are small.

Observed photometric variations are represented by sine fits to light curves for the M-type Mira stars R Car, R Hya, R Oct, R Vir, RR Sco, T Col and T Hor (for more details see Bladh et al., 2013). The difference between the observed loops are most likely due to differences in the fundamental properties of the stars, which is further explored in Bladh et al. (2013), Bladh et al. (2015), and Höfner et al. (2016).

Spectra for selected snapshots of every model are calculated, using a resolution of R $=10,000$ and covering a wavelength region between 0.3$$\rm\upmu m and 25$$\rm\upmu m. We follow the same method here as described in Nowotny et al. (2011) and Aringer et al. (2016). Photometric filter magnitudes are calculated from these spectra, following the Bessell system (Bessell & Brett, 1988; Bessell, 1990). To be consistent with the observations, sine curves are fitted to the synthetic photometry variations of the models to produce loops in the $(J-K)$ vs. $(V-K)$ plane.

The observed loops can be compared to loops calculated from the models, as the dynamical model should be able to reproduce the general characteristics of these variations.

2.2.2 Line profiles in high-resolution spectra

The velocity field in the line formation region will affect the line profiles of molecular lines originating in the dynamical atmosphere. While the light curves computed in the previous section reflect the luminosity variations, the line profiles of molecular lines are affected by the velocity fields in the respective line-forming region. By calculating photometric variations and synthetic line profiles for the models and then comparing them with observations, it should be possible to deduce information about a potential phase shift between the luminosity and radial variation.

Nowotny et al. (2005a, b, 2010) studied such line formation and velocity effects on various molecular features in C-type AGB star models. It was found that DARWIN models were able to reproduce various spectral features observed, and that there are three system of lines that are particularly interesting: the vibration-rotation CO lines $\Delta v=1,2,3$. These lines probe three different regions in the atmosphere as the temperature where these lines form are roughly 350-500K for CO $\Delta v=1$, 800-1500K for CO $\Delta v=2$, and 2200-3500K for CO $\Delta v=3$. The CO $\Delta v=1$ lines will then probe the outflow, while the CO $\Delta v=2$ lines will probe the dust-forming regions and the CO $\Delta v=3$ lines probe the dust-free layers (e.g. Sect. 5.1 of Nowotny et al., 2010).

The CO $\Delta v=3$ lines are formed in a region where the velocity field is dominated by effects of the shock waves, so the shape of the lines are determined by the propagation of the shock wave, which in turn depends primarily on radial variation due to pulsation. The movement of the mass-shells in such a line-forming region will then be imprinted on the line profile. A schematic of this scenario can be seen in Fig. 3. When material is infalling the line profile will be red-shifted, seen at $\phi\sim\phi_{red}$ to the left in Fig. 3. As an outwards propagating shock wave hits the infalling material, there will be a line splitting due to the gas layers in the line-forming region moving both inwards (red-shifted component) and outwards (blue-shifted component). This is seen at $\phi\sim\phi_{split}$. The line will be blue-shifted at $\phi\sim\phi_{blue}$ when all the layers are moving outwards.

Observations of Doppler-shifted lines in AGB stars will then directly probe the velocity field, and can be compared to model results. In Liljegren et al. (2016) such lines were synthesised for different variations of the inner boundary condition for the case of a C-star atmospheric model, and it was found that line profile variations, particularly the CO $\Delta v=3$ line, combined with information about the visual phase, can be used as a diagnostic tool. This is tested here for the case of M-stars.

Line profiles for the CO line $\Delta v=3$ (CO 5-2 P30 at 1.6573 $\rm\upmu m$) at different times during a pulsation cycle are produced for DARWIN models with a range of boundary conditions. A resolution of R = 300,000 is used, covering a wavelength region between 1.6566$$\rm\upmu m and 1.6578$$\rm\upmu m. Frequency dependent opacities are produced using the COMA code, with assumptions that the line shapes are described by Doppler profiles. Gas velocities are taken into account when performing the radiative transfer because the velocity field, with both infalling and outflowing gas, defines the line profiles. Opacities due to dust sources are accounted for, however, they do not significantly contribute to the result as these dust species are very transparent. For a detailed description of the spectral synthesis see Nowotny et al. (2010).

Synthetic RV curves are also calculated for each model, using the CO line $\Delta v=3$ synthesis, and compared to observed radial velocity curves from Lebzelter & Hinkle (2002).

2.3 Bolometric phase $\phi_{bol}$ and visual phase $\phi_{\rm v}$

Two different measures of time are used throughout this paper: visual phase $\phi_{\rm v}$ and bolometric phase $\phi_{bol}$. Cycle dependent observations are usually linked to the star’s visual phase $\phi_{\rm v}$. For the model calculations we use the bolometric phase $\phi_{bol}$, which is always known, as one of the model outputs is the luminosity lightcurve.

For direct comparison between models and observations, when phases are important (Sect. 4.3), $\phi_{\rm v}$ is calculated for the models. For O-rich Miras the bolometric phase lags behind the visual phase $\phi_{\rm v}$, typically by $\phi_{bol}(max)-\phi_{\rm v}(max)\approx 0.1$ (cf. Appendix A of Nowotny et al., 2010). This is comparable to what is found in the simulations, where the differences between the visual and bolometric phase typically were in the range $0.05-0.15$ for the model atmospheres in this paper.

3 Resulting models

3.1 Overview of the wind properties

Changing the inner boundary may have large effects on efficiency of the wind driving in the resulting model atmosphere. With a mass-loss rate $\dot{M}$, which is due to the radiation pressure, the mass $\dot{M}\delta t$ ejected over a time interval $\delta t$ will acquire the momentum $\dot{M}\delta tv_{\infty}$ for a wind velocity $v_{\infty}$, by absorbing a fraction $\beta$ of the momentum available from the luminosity $L_{\star}\delta t/c$. Therefore $\dot{M}\delta tv_{\infty}=\beta L_{\star}\delta t/c$, which can be solved for $\beta$ as

[TABLE]

Here $\beta$ describes the fraction of momentum of stellar radiation that goes into driving the wind, and can be calculated from the outputs of the DARWIN models. The upper left panel of Fig. 4 shows $\beta$ for all models. Different boundary conditions result in a wide range of values of $\beta$, indicating that the form of the inner boundary is directly correlated with the efficiency of the wind driving in these models. There seem to be systematic effects; a group of models with $\Delta\phi_{s}\in[-0.15,0.2]$ and $\Delta\phi_{p}\in[0,0.2]$ are very efficient in wind driving. A large portion of the models with $\Delta\phi_{p}>0.3$, on the other hand, are inefficient in driving the wind and produce much lower values of $\beta$.

The global properties of mass-loss rate and wind velocity are also highly dependent on the choice of the inner boundary conditions. The lower left and lower right panels of Fig. 4 show mass-loss rate and wind velocity plotted separately. Again some systematic effects can be seen; the models with $\Delta\phi_{s}\in[-0.1,0.1]$ and $\Delta\phi_{p}\in[0,0.2]$ result in the highest wind velocities, with values around $v_{\infty}=12-16$ kms*-1*. Velocities decrease with higher $\Delta\phi_{tot}$, with the lowest velocities being $v_{\infty}\approx 2$ kms*-1*. There are also significant effects of changing the shape of the luminosity variation, specified by $\Delta\phi_{s}$.

Many of the models have a mass-loss rate between $1-3\times 10^{-6}$ $M_{\odot}$ yr*-1*, and do not seem to follow the same trends as wind velocity and $\beta$. There are several models with $\Delta\phi_{s}\in[0.1,0.2]$ that have a relatively high mass-loss rate, but do not reach very high wind velocities. While a significant amount of material is lifted during the wind driving process for these models, the acceleration of this material is not efficient. The resulting wind velocities are too small to be realistic, in combination with the mass-loss rates (shown to the left in Fig. 6). The lowest mass-loss rate is $\sim 10^{-8}$ $M_{\odot}$ yr*-1*, which is two orders of magnitude smaller than the model with highest mass-loss rate.

The upper right panel of Fig. 4, shows the radius of the silicate dust grains for all models, which is consistently between 0.35 and 0.4 $\rm\upmu m$. As the dust radius $a_{gr}$ is related to the condensation degree of Si (denoted by $f_{cond}$) as $f_{cond}\propto a_{gr}^{3}$, this relatively small spread in grain radii still represents a significant difference in the degree of condensation. The models with the smallest grain radii have condensation degrees of $f_{cond}\sim 0.15$, while the models with the largest grains have a $f_{cond}\sim 0.3$ It is the models with the largest grain radii that also reach the highest wind velocities. This sets them apart from the group of models with high mass-loss rate but relatively low wind velocities ($\Delta\phi_{s}\in[0.1,0.2]$), for which the grain radii are comparably small.

3.2 Atmospheric structure and dynamics due to different boundary conditions

The large variety in mass-loss rates and wind velocities are due to the effect the inner boundary has on the structure of the atmosphere. This is discussed in detail for carbon stars in Liljegren et al. (2016), and while the wind-driving dust species in C-stars is different from that of M-stars the mechanisms at play, determining the atmospheric structures and the wind properties, are comparable. A shorter explanation is given here by comparing a model with efficient wind driving (left panel in Fig. 5) to one with poor wind driving (right panel of Fig. 5).

The left plot in Fig. 5 shows the atmospheric dynamics of the original model, where each line tracks the movement of a Lagrangian mass shell with time at different depths. This plot is illustrative of the processes contributing to a pulsation-enhanced dust-driven wind in AGB stars, indicating three distinctly different regions in the dynamical atmosphere. The dust-free region close to the surface at $R<2R_{\star}$ is dominated by the pulsation. When infalling matter collides with the outflowing gas, shock waves develop and propagate outwards. The material would follow ballistic trajectories if not for the onset of the wind by radiation pressure on the dust. Dust condenses at distances of $R\sim 2R_{\star}$. In this dust-forming region the dynamical behaviour changes from strictly periodic to more sporadic, with some matter falling back onto the star and some material being accelerated outwards. A steady outflow is found beyond $R\sim 3R_{\star}$.

Ideally, for effective wind driving, the dust formation should take place in the wake of the propagating shock wave, where the dust can accelerate outward moving dense material further, leading to a strong wind. This is the case for the original model, seen in the left plot of Fig.5, where dust is formed at $\phi_{bol}\sim 0.5$. The dust formation region correlate with the lowest temperatures in the atmosphere, which in turn correlate with the luminosity minimum.

The most important change in atmospheric dynamics, when changing the inner boundary, is related to the phase of dust condensation. When shifting the luminosity variation compared to the radial movement of the upper stellar layers, which is the consequence of introducing new inner boundary conditions, the timing of dust condensation is also shifted. If the luminosity minimum and therefore the temperature minimum occurs such that dust is not formed when gas is moving outwards in a wake of a shock, but rather when material is falling back towards the star (seen in the right panel of Fig. 5), the wind driving becomes highly ineffective. This leads to significantly lower wind velocities and mass-loss rates.

Another important aspect, which also determines the effectiveness of the wind, is when the luminosity maximum occurs with respect to dust formation and shock wave propagation. A luminosity maximum earlier in the pulsation cycle leads to a larger radiative acceleration on the dust, and therefore a higher wind velocity and mass-loss rate. On the contrary, a luminosity maximum that is late with respect to the propagating shock wave means less acceleration and a lower wind velocity and mass-loss rate.

The combination of these two effects leads to the vast variety of different behaviours; there are diverse wind velocities and mass-loss rates for the different inner boundary conditions. This is the reason why some models are much more efficient in driving a wind.

4 Comparison to observables

The models are evaluated via comparison to observations using three criteria: i) a combination of wind velocity and mass-loss rate, ii) colour-colour loops, and iii) RV curves derived from CO $\Delta v=3$ vibration-rotation lines.

4.1 Wind velocity and mass-loss rate

Some of the most studied properties of AGB stars are the wind velocities and mass-loss rates, and DARWIN models should reproduce these properties realistically.

An overview of the wind velocities and mass-loss rates of the models is shown in Fig. 6 (pale green, blue, and dark blue circles and the original model as an orange star). This is compared to a linear fit to observational wind velocities against log of mass-loss rates for M-stars using observations from Olofsson et al. (2002) and González Delgado et al. (2003) (crosses are observations, the dashed line is the fit to the observations). The observational uncertainties on the wind velocity vs. mass-loss rate relationship are dominated by uncertainties on the mass-loss rates, which are determined using CO multi-transitional line observations. The grey area plotted in Fig. 6 indicates the standard deviation from the linear fit in the mass-loss rates of the observations, which are comparable to the uncertainties found in Ramstedt et al. (2008) where the reliability of mass-loss rate estimates for AGB stars are investigated in detail. In the left panel of Fig. 6 the dark grey area marks the region with deviations from the linear fit being below 1$\sigma$, while the light grey area denotes a deviation between 1 and 2$\sigma$.

As previously mentioned changing the inner boundary may have large consequences for both the wind velocities and for mass-loss rates. The wind velocities range between approximately $1$ kms*-1* to almost $16$ kms*-1*, while the mass loss can change over two two orders of magnitude, between $3\times 10^{-6}$ and $10^{-8}$ $M_{\odot}$ yr*-1*.

For higher wind velocities, the wind seems to saturate at some maximum mass-loss rate, in this case around $1.5\times 10^{-6}$ $M_{\odot}$ yr*-1*. This leads to a flat distribution in mass-loss rate above $10$ kms*-1*. The changes to the inner boundary thus have little consequence for the mass-loss rate beyond this point; however, the wind velocity is still sensitive to such changes. Below $10$ kms*-1* there is a wider spread, but the models seem to have mass-loss rates that are too high for such low wind velocities, which is not realistic when compared to the observations.

The right panel in Fig. 6 shows the comparison between the observations and the models for different combinations of boundary conditions. As seen, there is a group of models with $\Delta\phi_{tot}\in[0.0,0.3]$ that produces realistic combinations of wind properties. The models that did not converge (red) are relatively far from this group of models.

There are in total 30 models that simultaneously reproduce a realistic combination of mass-loss rates and wind velocities within $1\sigma$. Most of these models have combinations of $\Delta\phi_{s}\in[-0.2,0.15]$ and $\Delta\phi_{p}\in[0,0.3]$, with $\Delta\phi_{tot}\in[0,0.3]$. It should be noted that the model with $\Delta\phi_{s}=-0.2$ and $\Delta\phi_{p}=0.7$, one of the most extreme models, has a very low mass-loss rate and a very low wind velocity(\dot{M}\approx 10^{-7.7}$$M_{\odot}yr*-1*, $v_{\infty}\approx 1$kms*-1*), and does not reproduce Mira-like properties. This model is therefore excluded, even though it technically falls within the $1\sigma$ range.

The 29 models left are used for further analysis and to investigate the colour-colour diagrams and high-resolution spectra. The model with the standard boundary conditions is among these 30 stars, and it produces a realistic mass-loss rate and wind velocity combination. This is rather reassuring, considering that the stellar parameters of this model emerged from a study based on C-rich dynamical models (with a different wind-driving species). As mentioned in Sect. 2.1.2, the pulsation properties of stars that have different chemical compositions but otherwise similar stellar properties should be comparable.

4.2 Loops in colour-colour diagrams

The photometric comparison with observations is done by calculating synthetic colours for the 29 models that produced satisfying mass-loss rates and wind velocities. Fig. 7 shows the loops formed in the $(J-K)$ vs. $(V-K)$ diagram after sine curves are fitted to the photometric variations, for both observed stars (red and blue loops) and for synthetic colours (grey loops). The model with the original boundary condition is shown in orange.

Models with higher $T_{\rm eff}$ and lower mass-loss rates, which have been studied in previous papers, fell mostly into the area covered by loops shown in red (see e.g. Fig. 8 in Höfner et al., 2016). Model M2, with its lower $T_{\rm eff}$ and higher mass-loss rate, on the other hand, is more comparable to R Oct (shown in blue), which has the highest mean ($J-K$) and ($V-K$) in the observed sample. The ($J-K$) values are commonly considered to be an indicator of $T_{\rm eff}$, which is also reflected by models (see Fig. 7 in Bladh et al., 2013). The ($V-K$) values depend on the strength of molecular features (mostly TiO), which in turn are strongly affected by the variable structure of the dynamical atmosphere and the condensation distance of the wind-driving dust grains (see the discussion in Bladh et al., 2013). Observations and models both loop in an anti-clockwise direction in this diagram, indicating qualitatively similar phase lags between light curves in the different photometric bands.

In an attempt to more quantitatively examine how well the observations and the model results agree and to see if any models significantly deviate, we describe the loops as ellipses in the $(J-K)$ vs. $(V-K)$ plane. An ellipse is defined by five quantities: the position of its centre $(x_{0},y_{0})$, the semi-major axis $a$, the semi-minor axis $b$ and the angle of inclination $\alpha$ . The shape of an ellipse is described by its eccentricity as $e=\sqrt{1-b^{2}/a^{2}}$.

The centre positions of the loops of the 29 models, which were selected based on their wind properties, are within the range of observed ($J-K$) and ($V-K$) colours for larger observational samples (mostly single epoch data, see Bladh et al., 2013), and can therefore not be used to put extra constraints on the boundary conditions. They will not be discussed further. Eccentricity, semi-major axis $a$, the semi-minor axis $b$ are plotted against the angle of inclination $\alpha$ for model loops and for observed loops in Fig. 8, again with observations in red and blue and models in grey with the original model as the orange star.

The left panel in Fig. 8 shows the eccentricities for the observed loops, which are all very close to 1. This very high eccentricity is due to a small variation in the $(J-K)$ plane compared to the variation in the $(V-K)$ plane, and common for all the observed loops. The angle of inclination $\alpha$ is between $0^{\circ}$ and $3^{\circ}$, meaning the observed loops are very narrow and almost flat. The model loops have a spread in both eccentricity and angle of inclination; however, some of the models overlap very well with the observations.

The middle and right panels of Fig. 8 show the semi-major and semi-minor axis respectively. Again, there is a spread for the model loops, but it seems to be realistic as this variation is similar to the spread between the different observations. R Oct seems to be a good fit for both the semi-major and semi-minor axis, and the eccentricity. This indicates that models can reproduce realistic $(J-K)$ (semi-minor) and $(V-K)$ (semi-major) variations simultaneously.

There are no models that obviously deviate from the observations. The spread of the different properties for the models is comparable to the spread of the observations. Therefore, all 29 models are used in the further analysis of the line profile variation.

4.3 Line profiles in high-resolution spectra

High-resolution spectra were calculated for the sample of 29 models that produced realistic mass-loss rates and wind velocity combination as well as realistic colour variations. Synthetic RV-curves from the spectra are compared to observed RV curves for Mira AGB stars.

4.3.1 Investigating the inner atmosphere with CO $\Delta v=3$

The CO $\Delta v=3$ lines originate deep in the atmosphere below the dust-forming layers where the radial expansion and compression of the upper stellar layers induce strong shock waves. The movement of the mass-shells in such a line-forming region will then be imprinted on the line profile, as described in Sect. 2.2.2 and illustrated in Fig. 3.

This behaviour is observed for CO $\Delta v=3$ lines in AGB stars. The first column of Fig. 9 shows time series observation of CO second overtone lines (average of 10-20 unblended lines) of $\chi$ Cyg, a S-type Mira variable. The $\phi_{split}$ of $\chi$ Cyg (described in Fig. 3) occurs at $\phi_{\rm v}\sim 0$. This line splitting at maximum visual light seems to be a ubiquitous feature of for Mira AGB stars.

Such line splitting and the resulting RV curve can also be synthesised using DARWIN models. The line profiles for the standard model, $\Delta\phi_{p}=0$, $\Delta\phi_{s}=0$, are shown in the two rightmost panels of Fig. 9 at different resolutions. The shape of the line changes during one pulsation cycle; when material is outflowing the line has a strong blue-shifted component ($\phi_{\rm v}\approx 0.65-0.0$), while the infalling material results in a red-shifted component ($\phi_{\rm v}\approx 0.2-0.4$). The $\phi_{split}$ in such a model occurs at $\phi_{\rm v}\sim 0.8$, which is too early compared to observations.

There is very little difference between the models with different boundary conditions concerning the actual line shape during a cycle other than the timing of different features, such as the line-splitting, occurring at different $\phi_{\rm v}$. This can be seen by comparing panels four and five (the original model, with $\Delta\phi_{tot}=0.0$) with panels two and three ($\Delta\phi_{p}=0.2$, $\Delta\phi_{s}=0.0$, with $\Delta\phi_{tot}=0.2$). The line profiles of the two models are similar, but different line features occurs at different $\phi_{\rm v}$. This indicates that changing the luminosity boundary has few consequences for the shape of the line, which in turn demonstrates that the material at these radial distances is not subjected to significant radiation pressure.

The main characteristics of the observed line, with a line-splitting around $\phi_{\rm v}=0.0$ and the overall shift of the main component of the line, can be reproduced if a model with $\Delta\phi_{tot}=0.2$ is used (such as shown in panels two and three). The velocity field of this model is then in overall agreement with that of the observed star. There are some differences however; the red component during the line splitting is significantly weaker for the models, probably due to a lower density of the infalling material.

4.3.2 Comparison of radial velocity curves

Radial velocity curves derived from observations of second overtone CO lines in Mira stars show a uniform behaviour: a discontinuous S-shaped curve, with both outflowing and infalling components at $\phi_{\rm v}=0.0$, $RV=0$ at around $\phi_{\rm v}=0.4$ and a velocity amplitude $\Delta RV\approx 25$ kms*-1* (Fig. 10, crosses and grey fit, which is the running mean over $\phi_{\rm v}=0.2$ for the observation). The shape and amplitude of the observed RV curve are very well reproduced in previous works, deriving RV curves from synthetic line spectra of C-star models; however, both the line splitting and $RV=0$ of these models occur earlier than $\phi_{\rm v}=0.0$ and $\phi_{\rm v}=0.4$ respectively (see e.g. Nowotny et al., 2005b, 2010). The same problem is found for the standard M-type model here, when calculating the RV curves from the CO $\Delta v=3$ line synthesis, indicated by orange dots in Fig. 10.

Models with a larger $\Delta\phi_{tot}$ approach observations: a line-splitting occurs closer to $\phi_{\rm v}=0.0$. This can be seen in Fig. 10, where the green dots represent a model with a larger total phase shift, $\Delta\phi_{tot}=0.2$ ($\Delta\phi_{p}=0.2$, $\Delta\phi_{s}=0.0$). Increasing $\Delta\phi_{tot}$ of the model will essentially shift the RV curve to the right. So while the inner luminosity boundary does not affect the line-forming regions directly, it does change the visual phase, which is important when comparing models to observations.

We want to compare the resulting RV curves for the 29 models that have realistic wind velocities and mass-loss rates with the mean RV curve from the CO-line observations of Mira stars. In this case it is not suitable to use a standard L2 norm (a measure of the mean quadratic difference between the points on two curves, for mathematical description see Appendix A.3) as an error estimation. Results from such an error might be misleading, due to the offset between the model curves and observational curve and due to the discontinuous nature of the RV curve. Instead we compare the timing of the synthetic and observed RV curves’ minimum ($\phi_{\rm v}(min)$), maximum ($\phi_{\rm v}(max)$), and the RV zero point ($\phi_{\rm v}(RV=0)$). These points are defined as $\phi_{\rm v}(s,min)$, $\phi_{\rm v}(s,max)$, and $\phi_{\rm v}(s,RV=0)$ for the synthetic RV curves and $\phi_{\rm v}(o,min)$, $\phi_{\rm v}(o,max)$, and $\phi_{\rm v}(o,RV=0)$ for the observed mean curves. The differences between the synthetic and observed RV curve for these points are defined as $\Delta max=|\phi_{\rm v}(o,max)-\phi_{\rm v}(s,max)|$, $\Delta min=|\phi_{\rm v}(o,min)-\phi_{\rm v}(s,min)|$ and $\Delta zero=|\phi_{\rm v}(o,RV=0)-\phi_{\rm v}(s,RV=0)|$. We then define the mean difference as

[TABLE]

The value of $\Delta\phi_{\rm v}$ is calculated for the 29 models. The results are shown in Fig. 11. To have a line splitting occurring close to what is observed, $\phi_{\rm v}\approx 0$, the model atmosphere must have a $\Delta\phi_{tot}\sim 0.2$. For models with lower $\Delta\phi_{tot}$ the line splitting occurs too early in the pulsation cycle.

This is important to note when modelling lines that vary due to shock propagation through the stellar layers. If correct phase information is to be retrieved, either models with $\phi_{tot}\sim 0.2$ need to be used, or a phase shift of $\Delta\phi_{\rm v}\sim 0.2$ (corresponding to need to a difference of $\sim 0.3$ if comparing $\phi_{\rm v}$ with $\phi_{bol}$) needs to be applied. The results also indicate that real Mira stars do have an inherent phase shift between the ascending luminosity and the radial movement of the surface.

5 Comparison to 1D pulsation models

The results found in this work are compared with theoretical predictions for the pulsation properties of AGB stars. Self-excited pulsation models of AGBs can be divided into two categories: 1D models and 3D models. One-dimensional self-excited radial pulsation models have been developed and worked on for decades (see e.g. Wood, 1974; Tuchman et al., 1979; Fox & Wood, 1982; Ostlie & Cox, 1986; Wood, 1990), and contemporary models include non-linear effects and turbulent pressure (e.g. Ireland et al., 2008, 2011; Wood, 2015). The drawback, however, is the treatment of convective energy transport, which is simulated with mixing length theory, and is mentioned as the major shortcoming in this generation of models (see e.g. Barthes, 1998, for a detailed discussion). Three-dimensional interior pulsation models with realistic hydro-dynamical treatment should provide a more realistic description of convection, and therefore the pulsation properties; however, these models are still in their infancy, and only a small part of the relevant stellar parameter space explored so far (see Freytag & Höfner, 2008; Freytag et al., 2017). Furthermore, possible 3D effects on mass loss have not yet been investigated.

Here we compare the different inner boundary conditions explored to results from 1D self-excited models, described in Bessell et al. (1996) where non-linear models representing the prototypical Mira, o Ceti, were calculated and explicit descriptions of the variations in the sub-photospheric layers are given. Two fundamental mode models (model Z and D) and one first overtone model (model E) were examined. However the results, and also later investigations, suggest that Mira variables are fundamental mode pulsators. The first overtone model did not achieve appropriate luminosity amplitude and a factor of six had to be applied in Bessell et al. (1996). So while comparisons to all three models are included here, the first-overtone model (model E) is probably not a realistic representation of Mira stars.

Bessell et al. (1996) provide Fourier coefficients to the time dependence of the luminosity and radius for the three models for the region with mean temperature of $\sim 4000-5000K$. This is close to the location of the inner boundary of the DARWIN models. These variations of the photospheric layers from the 1D pulsation models are therefore directly comparable to the inner boundary condition for the DARWIN models. A comparison is made between different DARWIN inner boundary conditions and the variations given by Bessell et al. (1996), which are rescaled to have the same amplitude and period. This comparison is evaluated by using the a relative least-squares error (LSE, the L2 norm), which is the mean quadratic difference at each point between the curves, as seen in Eq. 9, then normalised to the smallest value. The result for each model presented in Bessell et al. (1996) can be seen in Fig. 12.

The boundary conditions most similar to the two fundamental mode pulsators given by Bessell et al. (1996) are slightly asymmetric (with $\Delta\phi_{s}\sim-0.05$) with an offset (with $\Delta\phi_{p}\sim 0.1$), resulting in a total phase shift $\Delta\phi_{tot}\sim 0.05$. This is a smaller phase shift than found by comparing the DARWIN models with different boundary conditions and RV curves in Sect. 4.3.2. It overlaps with DARWIN models that produce realistic wind velocities and mass-loss rates, however. The first overtone pulsator model from Bessell et al. (1996) predicts asymmetry (with $\Delta\phi_{s}\sim 0.1$) and with an offset (with $\Delta\phi_{p}\sim 0.1$), resulting in $\Delta\phi_{tot}\sim 0.2$. This model is discarded as unrealistic in the Bessell et al. (1996).

6 Discussion

Under the assumption that the results for model M2, representing a prototypical Mira, can be generalised to other Mira stars we can draw conclusions about which inner boundary condition agrees best with observations and about what the consequences are of using a generic boundary condition, like that of the standard DARWIN models, for the resulting atmospheres.

In our original sample of 99 models with a broad range of inner boundary conditions, 29 models reproduced realistic combinations of mass-loss rates and wind velocities. There were some systematic effects; $\Delta\phi_{tot}$ ranged between 0 and 0.3, and most of the models had a slightly asymmetric shape with $\Delta\phi_{s}<0.0$. This sample could be reduced further, by comparing the model RV curves with compilations of CO $\Delta v=3$ RV curves for Mira stars.

The RV comparison indicates that a $\Delta\phi_{tot}\sim 0.2$ is needed for the splitting of this line to occur at correct $\phi_{\rm v}$. As no further information about the shape can be derived from this comparison, any model with $\Delta\phi_{tot}\sim 0.2$ should reproduce line splitting at the correct $\phi_{\rm v}$.

The $(V-K)$ vs $(J-K)$ colour-colour loops were also studied for the 29 models that produced satisfying mass-loss rates and wind velocities. While there was a spread in the ellipse properties, it was comparable to the same size as the spread in observed colour-colour loops. This is not always the case as models with unrealistic temperature structures and resulting molecule abundances in turn result in results in unrealistic loops, a fact that was explored in Bladh et al. (2013, 2015).

The similarities of the model loops here is probably a selection effect as photometry calculations are only performed on models with similar mass-loss rates and velocity properties. These models have similar atmospheric structures, which in turn result in similar colour-colour loops. Because no models were obviously unrealistic, no models were discarded using the colour-colour loop criterion.

The combined conclusion when looking at both wind velocity and mass-loss rate, and that of the RV curve comparison, the ideal inner boundary condition has a $\Delta\phi_{s}$ between -0.2, and 0.0, with $\Delta\phi_{tot}\sim 0.2$. There is thus a degeneracy when it comes to the shape of the luminosity variation.

It is not possible to pin down one ideal inner boundary condition from these tests, only that it should have $\Delta\phi_{tot}\sim 0.2$. However, the results do indicate that there is an inherent offset between the radial movement and the luminosity variation of Mira stars.

The standard boundary condition, which has previously been used extensively in grid calculations such as Eriksson et al. (2014) and Bladh et al. (2015), produces realistic wind velocities and mass-loss rate combinations for the M2 model. If phase information about the ascension of the shock wave is not of interest, the standard boundary condition can therefore be used in the models.

Rescaling the photospheric variations of the luminosity and the gas layers found by 1D radial pulsation models from Bessell et al. (1996), and comparing them with the inner boundary conditions used in this work we find an overlap with slightly asymmetric DARWIN models ($\Delta\phi_{s}\sim-0.05$) with a small $\Delta\phi_{tot}$ of 0.05. The DARWIN models with these inner boundary conditions do reproduce realistic combination of wind velocities and mass-loss rates and good colour-colour loops; however, the line splitting in the RV curves occurs too early in the pulsation cycle when compared to observed RV curves.

7 Summary and conclusions

In this paper we investigate the influence of the inner boundary conditions, used in the DARWIN models to simulate the observable variations of AGB stars, to evaluate the effects on the resulting atmosphere structure and wind properties.

The results can be summarised as follows:

• •

The DARWIN models are sensitive to the inner boundary, and using different boundary conditions can result in significant differences for the mass-loss rates (about three orders of magnitude) and wind velocities (about one order of magnitude). However, not all of the resulting models are realistic;

• •

The mass-loss rates seem to saturate (here at $1.5\times 10^{-6}$ $M_{\odot}$ yr*-1*), but changing the inner boundary will change the wind velocity;

• •

All the models with boundary conditions such that they produced realistic combinations of mass-loss rates and wind velocities, also resulted in synthetic colour-colour loops that agree with observations;

• •

There is an inherent phase difference between $\phi_{bol}$ and $\phi_{\rm v}$ for the DARWIN models of $\sim 0.1$. This has previously been suggested to be the case for observed stars;

• •

Phase information can be derived from RV curves; DARWIN models with a $\Delta\phi_{tot}\sim 0.2$ results in a line splitting phase that corresponds to observations;

• •

The 1D radial pulsation models of Bessell et al. (1996) indicate that the total phase difference $\Delta\phi_{tot}$ is smaller than $\Delta\phi_{tot}\sim 0.2$;

• •

Using the standard boundary condition results in realistic mass-loss rates, wind velocities and colour-colour loops in the case of M-stars; however, the line splitting of the CO dv=3 line will occur at the wrong $\phi_{\rm v}$.

It can therefore be concluded that grid studies of DARWIN models, such as Eriksson et al. (2014) and Bladh et al. (2015), that use the standard inner boundary condition should produce realistic mass-loss rates and wind expansion velocities.

Acknowledgements.

The observational results (RV data as well as FTS spectra of $\chi$ Cyg, S Cep, and W Hya) were kindly provided by Th. Lebzelter and K. Hinkle. S.H. would like to acknowledge the support from the Swedish Research Council (Vetenskapsrådet). The computations were performed on resources provided by the Swedish National Infrastructure for Computing (SNIC) at UPPMAX.

Appendix A Description of the boundary condition

A.1 Inner boundary

Originally in the DARWIN code the effects of stellar pulsation were described by two variable physical quantities at the inner boundary: $R_{in}(t)$ and $L_{in}(t)$. The first $R_{in}(t)$, is the variation in the innermost gas layer. This lower boundary layer is impermeable, meaning no gas flows across it. The second, $L_{in}(t)$ describes the variation in the luminosity. The inner boundary is places below the stellar photosphere, at around $R\approx 0.9R_{\star}$. The $R_{in}(t)$ variation has the form

[TABLE]

where $\Delta u_{p}$ is the velocity amplitude, P is the pulsation period, and $R_{0}$ the average radial distance of the boundary. This corresponds to a gas velocity of

[TABLE]

Assuming a constant flux at the inner boundary the luminosity becomes $L_{in}(t)\propto R^{2}_{in}(t)$. To better match observations of the bolometric flux variation a free parameter $f_{L}$ was introduced in Gautschy-Loidl et al. (2004), so that the amplitude of luminosity can be adjusted separately from the radial amplitude. The original form of the luminosity is

[TABLE]

where $L_{0}$ is the average luminosity of the star. In this paper we alter the luminosity variation, but keep the form of the radial variation.

If a phase shift $\Delta\phi_{p}$ is introduced between $R_{in}(t)$ and $L_{in}(t)$, but the shape is kept unchanged, the luminosity variation becomes

[TABLE]

The asymmetric luminosity variation is achieved by describing the luminosity variation at the inner boundary by a modified, smoothed Fourier saw-tooth wave curve as

[TABLE]

where k is the amplitude set to match the amplitude from Eq. (A.1), P is the pulsation period and $w(\Delta\phi_{s})$ is a smoothing factor. If $w(\Delta\phi_{s})=1$ we retain the saw-tooth shape; however, when $w$ is increased $\Delta L_{in}$ will approach a sinusoidal curve. As seen in Fig. 1, the $\Delta\phi_{s}$ is a measure of how asymmetric the curve is, specifying how far from the original $L_{in}$ the phase of the maximum has been shifted. This measure is analogous to $\Delta\phi_{p}$.

This form also allows for different values of $\Delta\phi_{p}$.

A.2 Estimation of w

The parameter $w$, from Eq. (7), is a function of $\Delta\phi_{s}$ and the relationship between $w$ and $\Delta\phi_{s}$ can be seen in Fig. 13. For a simpler extraction of the $w$-value a polynomial is fitted, of the form

[TABLE]

A good fit was found for M = 5. Table 2 shows the constants for this fit.

A.3 Number of terms used

When using the Fourier series from Eq. (7), the number of terms $N$ needed to avoid significant unwanted overshooting close to the discontinuity (Gibbs phenomenon) depends on the value of $\Delta\phi_{s}$. With a higher $\Delta\phi_{s}$, the luminosity variation will be more asymmetric and more prone to overshooting, which requires a larger $N$.

To estimate the $N$ number of Fourier terms needed for convergence at different $\Delta\phi_{s}$, $\Delta L_{in}$ for one period calculated using $N$ Fourier terms is compared to one period with $N+1$ terms. When the difference between $\Delta L_{in}(N)$ and $\Delta L_{in}(N+1)$ becomes small, it is assumed that $\Delta L_{in}(N)$ has converged.

The difference between the two curves is measured with the least-squares error, defined as

[TABLE]

with $f_{i}=\Delta L_{in}(N+1)_{i}$ and $g_{i}=\Delta L_{in}(N)_{i}$. This is a measure of the the average difference between two points on the curves, and is usually denoted as the least-squares error or the L2 norm. Convergence is assumed when $L2<0.001$.

An overview of the number of Fourier terms $N$ used in this work, for each value of $\Delta\phi_{s}$, is seen in Table 3 and Fig. 14.