EDP Sciences
Free Access
Issue
A&A
Volume 606, October 2017
Article Number A6
Number of page(s) 15
Section Stellar atmospheres
DOI https://doi.org/10.1051/0004-6361/201731137
Published online 22 September 2017

© ESO, 2017

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 fromdynamical wind models (DARWIN models, Sect. 2.1), and ameasure of the general dynamics. The models are evaluated bycomparing the velocity and mass-loss rate combinations withobservations from Olofssonet al. (2002) andGonzález Delgadoet al. (2003) inSect. 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 (JK) vs. (VK) 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 × 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

thumbnail Fig. 1

Different variations of the inner boundary condition used, showing luminosity variation (black) and the radial expansion and contraction (dotted blue). Upper panel: the original boundary condition, where the shape of the variation is sinusoidal for both luminosity and radial variation, and where the phase is locked. Second panel: shape of the luminosity is unchanged, however an offset between the maximum luminosity and largest expansion is introduced, indicated by Δφp ≠ 0. Third panel: shape of the luminosity is made asymmetric, indicated by Δφs ≠ 0. Lower panel: case with both offset (Δφp ≠ 0) and asymmetry (Δφs ≠ 0).

Open with DEXTER

The spatial range of the DARWIN models, for models that develop a wind, reaches from just below the photosphere (~0.9 R) out to where the wind has reached its terminal velocity v (~20 − 30 R). 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; ΔRin ∝ sin(2πt/P) for the radial variation and 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 Δφp between the radial variation and luminosity variation and ii) by changing the shape of the luminosity variation, making it increasingly asymmetric (measured by Δφ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: Δφp for case i) (panel 2 in Fig. 1) and Δφs for the asymmetric case ii) (panel 3 in Fig. 1). The parameter Δφs is analogous to Δφ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 Δφ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 Δφtot = Δφp + Δφs (see panel 4 in Fig. 1).

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

2.1.2. Model set

Table 1

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

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 ≈ 1 R), 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 Δφtot ∈ [ 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.

thumbnail Fig. 2

Overview of the calculated models, showing the Δφs and Δφp used. The points represent each of the models; grey background means the model finished while the red background means the model did not converge (Sect. 2.1.2). The lines show Δφtot = Δφs + Δφp.

Open with DEXTER

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 ξt = 2.5 km s-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 (JK) vs. (VK). Typically, observed M-stars show large variations in (VK) during a cycle, caused by changes in the molecule abundances (mainly TiO), while the variations in (JK) 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, 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μm and 25μ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 (JK) vs. (VK) 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 Δv = 1,2,3. These lines probe three different regions in the atmosphere as the temperature where these lines form are roughly 350–500 K for CO Δv = 1, 800–1500 K for CO Δv = 2, and 2200–3500 K for CO Δv = 3. The CO Δv = 1 lines will then probe the outflow, while the CO Δv = 2 lines will probe the dust-forming regions and the CO Δv = 3 lines probe the dust-free layers (e.g. Sect. 5.1 of Nowotny et al. 2010).

The CO Δ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 φ ~ φ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 φ ~ φsplit. The line will be blue-shifted at φ ~ φblue when all the layers are moving outwards.

thumbnail Fig. 3

Schematic of pulsating atmospheric layers, and lines created in this region at different times. The left panel: shows mass-shells, and the colour indicates whether matter is infalling (red) or outflowing (blue). At times when the gas layers are moving inwards (φ ~ φred) the Doppler-shifted line will have a red component (first panel to the right), while outflowing material (φ ~ φblue) leads to a blue-shifted line (last panel to the right). When a shock waves moves through the line-forming region (φ ~ φsplit) there will be both a red component and a blue component (second panel to the right).

Open with DEXTER

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 Δ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 Δv = 3 (CO 5–2 P30 at 1.6573 μ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μm and 1.6578μ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 Δv = 3 synthesis, and compared to observed radial velocity curves from Lebzelter & Hinkle (2002).

2.3. Bolometric phase φbol and visual phase φv

Two different measures of time are used throughout this paper: visual phase φv and bolometric phase φbol. Cycle dependent observations are usually linked to the star’s visual phase φv. For the model calculations we use the bolometric phase φ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), φv is calculated for the models. For O-rich Miras the bolometric phase lags behind the visual phase φv, typically by φbol(max) − φv(max) ≈ 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

thumbnail Fig. 4

Overview of the calculated models, showing the dependence on Δφs, Δφp and Δφtot. Each square represents one model, and the colours indicate physical properties. The red squares are models that did not converge (see Sect. 2.1.2). Upper left: β, the fraction of momentum available from the luminosity that has gone into driving the wind. Upper right: final radius of the dust grains. Lower left: Log of the mass-loss rate. Lower right: wind velocity of the models.

Open with DEXTER

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 , which is due to the radiation pressure, the mass δt ejected over a time interval δt will acquire the momentum δtv for a wind velocity v, by absorbing a fraction β of the momentum available from the luminosity Lδt/c. Therefore δtv = βLδt/c, which can be solved for β as (1)Here β 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 β for all models. Different boundary conditions result in a wide range of values of β, 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 Δφs ∈ [ − 0.15,0.2 ] and Δφp ∈ [ 0,0.2 ] are very efficient in wind driving. A large portion of the models with Δφp> 0.3, on the other hand, are inefficient in driving the wind and produce much lower values of β.

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 Δφs ∈ [ − 0.1,0.1 ] and Δφp ∈ [ 0,0.2 ] result in the highest wind velocities, with values around v = 12 − 16 km s-1. Velocities decrease with higher Δφtot, with the lowest velocities being v ≈ 2 km s-1. There are also significant effects of changing the shape of the luminosity variation, specified by Δφs.

Many of the models have a mass-loss rate between 1 − 3 × 10-6M yr-1, and do not seem to follow the same trends as wind velocity and β. There are several models with Δφs ∈ [ 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 ~10-8M 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 μm. As the dust radius agr is related to the condensation degree of Si (denoted by fcond) as , 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 fcond ~ 0.15, while the models with the largest grains have a fcond ~ 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 (Δφs ∈ [ 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< 2 R 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 ~ 2 R. 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 ~ 3 R.

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 φbol ~ 0.5. The dust formation region correlate with the lowest temperatures in the atmosphere, which in turn correlate with the luminosity minimum.

thumbnail Fig. 5

Mass shell plots for two models. Pink represents where more than 1% of the available Si has condensed, and indicates the presence of dust. Left: standard model Δφp = 0.0, Δφs = 0.0; dust forms just behind the shock wave, where the density is high, leading to efficient wind driving. Right: model Δφp = 0.4, Δφs = 0.0, with poor wind driving. Dust forms when matter is falling back onto the star, in a low-density region, resulting in significantly less matter condensing into dust, and therefore less material is accelerated in the wind.

Open with DEXTER

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.

thumbnail Fig. 6

Left: mass-loss rates and wind velocities of models with phase shift at the inner boundary (circles), the original model (orange star) and observations (crosses). The dark grey area is the standard deviation σ of the mass-loss rates of the observations, and the light grey is 2σ. Right: each square represent a model; Δφs, Δφp and Δφtot indicated. The colours are the same as in the left panel; models <1σ in the mass-loss rate and velocity plot are in light green.

Open with DEXTER

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 Δ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σ, while the light grey area denotes a deviation between 1 and 2σ.

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 km s-1 to almost 16 km s-1, while the mass loss can change over two two orders of magnitude, between 3 × 10-6 and 10-8M yr-1.

For higher wind velocities, the wind seems to saturate at some maximum mass-loss rate, in this case around 1.5 × 10-6M yr-1. This leads to a flat distribution in mass-loss rate above 10 km s-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 km s-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 Δφtot ∈ [ 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σ. Most of these models have combinations of Δφs ∈ [ − 0.2,0.15 ] and Δφp ∈ [ 0,0.3 ], with Δφtot ∈ [ 0,0.3 ]. It should be noted that the model with Δφs = − 0.2 and Δφp = 0.7, one of the most extreme models, has a very low mass-loss rate and a very low wind velocity ( ≈ 10-7.7M yr-1, v ≈ 1 km s-1), and does not reproduce Mira-like properties. This model is therefore excluded, even though it technically falls within the 1σ 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. Figure 7 shows the loops formed in the (JK) vs. (VK) 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 Teff 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 Teff and higher mass-loss rate, on the other hand, is more comparable to R Oct (shown in blue), which has the highest mean (JK) and (VK) in the observed sample. The (JK) values are commonly considered to be an indicator of Teff, which is also reflected by models (see Fig. 7 in Bladh et al. 2013). The (VK) 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 (JK) vs. (VK) plane. An ellipse is defined by five quantities: the position of its centre (x0,y0), the semi-major axis a, the semi-minor axis b and the angle of inclination α. The shape of an ellipse is described by its eccentricity as .

The centre positions of the loops of the 29 models, which were selected based on their wind properties, are within the range of observed (JK) and (VK) 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 α 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.

thumbnail Fig. 7

Colour-colour loops in the (JK) vs. (VK) plane, with models in grey, the standard model with no phase shift in orange, and observations of seven Mira variables shown in red and blue (see Bladh et al. 2013, for more details). R Oct (blue) seems to be very similar to the synthetic loops.

Open with DEXTER

thumbnail Fig. 8

Ellipse properties of the synthetic (grey) and observed (red and blue) colour-colour loops, in the (VK) vs. (JK) plane against the angle of inclination. Left: eccentricity, a measure of the shape. Middle: semi-major axis, (VK). Right: semi-minor axis, (JK).

Open with DEXTER

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 (JK) plane compared to the variation in the (VK) plane, and common for all the observed loops. The angle of inclination α is between and , 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 (JK) (semi-minor) and (VK) (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.

thumbnail Fig. 9

Line profiles from observations and from synthesis. The number is φv of that panel. Column 1 from the left: time series of an average of 10–20 unblended CO Δv = 3 lines, using FTS spectra of χ Cyg, taken from Lebzelter et al. (2001). Observed heliocentric velocities were converted to systemic velocities RV* using heliocentric centre of mass radial velocity CMRV = − 7.5 km s-1 from Hinkle et al. (1982). Column 2 and 3: synthetic CO Δv = 3 (5–2 P30) line profiles, from a model with Δφtot = 0.2, over a full cycle, shown for resolution 300 000 (Col. 2) and 70 000 (Col. 3). Column 4 and 5: synthetic CO Δv = 3 (5–2 P30) line profiles, from the original model, over a full cycle, shown for resolution 300 000 (Col. 4) and 70 000 (Col. 5).

Open with DEXTER

4.3.1. Investigating the inner atmosphere with CO Δv = 3

The CO Δ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 Δ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 χ Cyg, a S-type Mira variable. The φsplit of χ Cyg (described in Fig. 3) occurs at φv ~ 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, Δφp = 0, Δφ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 (φv ≈ 0.65 − 0.0), while the infalling material results in a red-shifted component (φv ≈ 0.2 − 0.4). The φsplit in such a model occurs at φv ~ 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 φv. This can be seen by comparing panels four and five (the original model, with Δφtot = 0.0) with panels two and three (Δφp = 0.2, Δφs = 0.0, with Δφtot = 0.2). The line profiles of the two models are similar, but different line features occurs at different φ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 φv = 0.0 and the overall shift of the main component of the line, can be reproduced if a model with Δφ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.

thumbnail Fig. 10

Radial velocities (RVs), with circles showing RVs derived from synthetic lines for two models, observation as crosses, and the grey line a running mean for the observations. The observed measurements are a compilation of RVs derived from FTS spectra for a large sample of different Mira stars, with the RVs converted to systematic velocities using the centre of mass radial velocity (CMRV) for each object (data from Lebzelter & Hinkle 2002).

Open with DEXTER

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 φv = 0.0, RV = 0 at around φv = 0.4 and a velocity amplitude ΔRV ≈ 25 km s-1 (Fig. 10, crosses and grey fit, which is the running mean over φ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 φv = 0.0 and φ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 Δv = 3 line synthesis, indicated by orange dots in Fig. 10.

Models with a larger Δφtot approach observations: a line-splitting occurs closer to φv = 0.0. This can be seen in Fig. 10, where the green dots represent a model with a larger total phase shift, Δφtot = 0.2 (Δφp = 0.2, Δφs = 0.0). Increasing Δφ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.

thumbnail Fig. 11

Comparison between synthetic RV curves of models and the running mean of the observations, seen in Fig. 10. The colour indicates the mean distance between the maximum and the minimum of that model’s RV curve and the running mean, described in Eq. (2). The grey squares are models that were previously discarded as they did not reproduce a realistic wind velocity and mass-loss rate combination.

Open with DEXTER

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 (φv(min)), maximum (φv(max)), and the RV zero point (φv(RV = 0)). These points are defined as φv(s,min), φv(s,max), and φv(s,RV = 0) for the synthetic RV curves and φv(o,min), φv(o,max), and φv(o,RV = 0) for the observed mean curves. The differences between the synthetic and observed RV curve for these points are defined as Δmax = | φv(o,max) − φv(s,max) |, Δmin = | φv(o,min) − φv(s,min) | and Δzero = | φv(o,RV = 0) − φv(s,RV = 0) |. We then define the mean difference as (2)The value of Δφ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, φv ≈ 0, the model atmosphere must have a Δφtot ~ 0.2. For models with lower Δφ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 φtot ~ 0.2 need to be used, or a phase shift of Δφv ~ 0.2 (corresponding to need to a difference of ~0.3 if comparing φv with φ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

thumbnail Fig. 12

Comparison between the 1D pulsation models, from Bessell et al. (1996), and the boundary conditions described in Sect. 2.1.1 and Appendix A, defined by Δφs and Δφp. The colours show the relative least-squares error (L2rel), and the points indicates our calculated models and lines show Δφtot = Δφs + Δφp. Left and middle: model Z and model D from Bessell et al. (1996), which are fundamental mode pulsator models of o Ceti. Right: model E from Bessell et al. (1996), a first-overtone pulsator model of o Ceti.

Open with DEXTER

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 ~4000 − 5000 K. 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. (A.7), 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 Δφs ~ − 0.05) with an offset (with Δφp ~ 0.1), resulting in a total phase shift Δφtot ~ 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 Δφs ~ 0.1) and with an offset (with Δφp ~ 0.1), resulting in Δφtot ~ 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; Δφtot ranged between 0 and 0.3, and most of the models had a slightly asymmetric shape with Δφs< 0.0. This sample could be reduced further, by comparing the model RV curves with compilations of CO Δv = 3 RV curves for Mira stars.

The RV comparison indicates that a Δφtot ~ 0.2 is needed for the splitting of this line to occur at correct φv. As no further information about the shape can be derived from this comparison, any model with Δφtot ~ 0.2 should reproduce line splitting at the correct φv.

The (VK) vs. (JK) 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 Δφs between –0.2, and 0.0, with Δφtot ~ 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 Δφtot ~ 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 (Δφs ~ − 0.05) with a small Δφ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, andusing different boundary conditions can result in significantdifferences for the mass-loss rates (about three orders ofmagnitude) 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 × 10-6M 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 φbol and φv for the DARWIN models of ~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 Δφtot ~ 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 Δφtot is smaller than Δφtot ~ 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 φ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.

Acknowledgments

The observational results (RV data as well as FTS spectra of χ 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.

References

Appendix A: Description of the boundary condition

Appendix 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: Rin(t) and Lin(t). The first Rin(t), is the variation in the innermost gas layer. This lower boundary layer is impermeable, meaning no gas flows across it. The second, Lin(t) describes the variation in the luminosity. The inner boundary is places below the stellar photosphere, at around R ≈ 0.9 R. The Rin(t) variation has the form (A.1)where Δup is the velocity amplitude, P is the pulsation period, and R0 the average radial distance of the boundary. This corresponds to a gas velocity of (A.2)Assuming a constant flux at the inner boundary the luminosity becomes . To better match observations of the bolometric flux variation a free parameter fL 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 (A.3)where L0 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 Δφp is introduced between Rin(t) and Lin(t), but the shape is kept unchanged, the luminosity variation becomes (A.4)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 (A.5)where k is the amplitude set to match the amplitude from Eq. (A.3), P is the pulsation period and wφs) is a smoothing factor. If wφs) = 1 we retain the saw-tooth shape; however, when w is increased ΔLin will approach a sinusoidal curve. As seen in Fig. 1, the Δφs is a measure of how asymmetric the curve is, specifying how far from the original Lin the phase of the maximum has been shifted. This measure is analogous to Δφp.

This form also allows for different values of Δφp.

Appendix A.2: Estimation of w

thumbnail Fig. A.1

Exact w = fφs), from the Fourier series Eq. (A.5), calculated with 10 000 terms and the polynomial fit from Eq. (A.6), using the terms from Table A.1. The blue circles represent the values of Δφs used in this work.

Open with DEXTER

Table A.1

Polynomial constants, for Eq. (A.6).

The parameter w, from Eq. (A.5), is a function of Δφs and the relationship between w and Δφs can be seen in Fig. A.1. For a simpler extraction of the w-value a polynomial is fitted, of the form (A.6)A good fit was found for M = 5. Table A.1 shows the constants for this fit.

Appendix A.3: Number of terms used

thumbnail Fig. A.2

Number of Fourier terms to reach an L2 value of less than 0.1%.

Open with DEXTER

When using the Fourier series from Eq. (A.5), the number of terms N needed to avoid significant unwanted overshooting close to the discontinuity (Gibbs phenomenon) depends on the value of Δφs. With a higher Δφ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 Δφs, ΔLin for one period calculated using N Fourier terms is compared to one period with N + 1 terms. When the difference between ΔLin(N) and ΔLin(N + 1) becomes small, it is assumed that ΔLin(N) has converged.

The difference between the two curves is measured with the least-squares error, defined as (A.7)with fi = ΔLin(N + 1)i and gi = ΔLin(N)i. This is a measure of 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 Δφs, is seen in Table A.2 and Fig. A.2.

Table A.2

Number of terms used in this work, for different values of Δφs.

All Tables

Table 1

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

Table A.1

Polynomial constants, for Eq. (A.6).

Table A.2

Number of terms used in this work, for different values of Δφs.

All Figures

thumbnail Fig. 1

Different variations of the inner boundary condition used, showing luminosity variation (black) and the radial expansion and contraction (dotted blue). Upper panel: the original boundary condition, where the shape of the variation is sinusoidal for both luminosity and radial variation, and where the phase is locked. Second panel: shape of the luminosity is unchanged, however an offset between the maximum luminosity and largest expansion is introduced, indicated by Δφp ≠ 0. Third panel: shape of the luminosity is made asymmetric, indicated by Δφs ≠ 0. Lower panel: case with both offset (Δφp ≠ 0) and asymmetry (Δφs ≠ 0).

Open with DEXTER
In the text
thumbnail Fig. 2

Overview of the calculated models, showing the Δφs and Δφp used. The points represent each of the models; grey background means the model finished while the red background means the model did not converge (Sect. 2.1.2). The lines show Δφtot = Δφs + Δφp.

Open with DEXTER
In the text
thumbnail Fig. 3

Schematic of pulsating atmospheric layers, and lines created in this region at different times. The left panel: shows mass-shells, and the colour indicates whether matter is infalling (red) or outflowing (blue). At times when the gas layers are moving inwards (φ ~ φred) the Doppler-shifted line will have a red component (first panel to the right), while outflowing material (φ ~ φblue) leads to a blue-shifted line (last panel to the right). When a shock waves moves through the line-forming region (φ ~ φsplit) there will be both a red component and a blue component (second panel to the right).

Open with DEXTER
In the text
thumbnail Fig. 4

Overview of the calculated models, showing the dependence on Δφs, Δφp and Δφtot. Each square represents one model, and the colours indicate physical properties. The red squares are models that did not converge (see Sect. 2.1.2). Upper left: β, the fraction of momentum available from the luminosity that has gone into driving the wind. Upper right: final radius of the dust grains. Lower left: Log of the mass-loss rate. Lower right: wind velocity of the models.

Open with DEXTER
In the text
thumbnail Fig. 5

Mass shell plots for two models. Pink represents where more than 1% of the available Si has condensed, and indicates the presence of dust. Left: standard model Δφp = 0.0, Δφs = 0.0; dust forms just behind the shock wave, where the density is high, leading to efficient wind driving. Right: model Δφp = 0.4, Δφs = 0.0, with poor wind driving. Dust forms when matter is falling back onto the star, in a low-density region, resulting in significantly less matter condensing into dust, and therefore less material is accelerated in the wind.

Open with DEXTER
In the text
thumbnail Fig. 6

Left: mass-loss rates and wind velocities of models with phase shift at the inner boundary (circles), the original model (orange star) and observations (crosses). The dark grey area is the standard deviation σ of the mass-loss rates of the observations, and the light grey is 2σ. Right: each square represent a model; Δφs, Δφp and Δφtot indicated. The colours are the same as in the left panel; models <1σ in the mass-loss rate and velocity plot are in light green.

Open with DEXTER
In the text
thumbnail Fig. 7

Colour-colour loops in the (JK) vs. (VK) plane, with models in grey, the standard model with no phase shift in orange, and observations of seven Mira variables shown in red and blue (see Bladh et al. 2013, for more details). R Oct (blue) seems to be very similar to the synthetic loops.

Open with DEXTER
In the text
thumbnail Fig. 8

Ellipse properties of the synthetic (grey) and observed (red and blue) colour-colour loops, in the (VK) vs. (JK) plane against the angle of inclination. Left: eccentricity, a measure of the shape. Middle: semi-major axis, (VK). Right: semi-minor axis, (JK).

Open with DEXTER
In the text
thumbnail Fig. 9

Line profiles from observations and from synthesis. The number is φv of that panel. Column 1 from the left: time series of an average of 10–20 unblended CO Δv = 3 lines, using FTS spectra of χ Cyg, taken from Lebzelter et al. (2001). Observed heliocentric velocities were converted to systemic velocities RV* using heliocentric centre of mass radial velocity CMRV = − 7.5 km s-1 from Hinkle et al. (1982). Column 2 and 3: synthetic CO Δv = 3 (5–2 P30) line profiles, from a model with Δφtot = 0.2, over a full cycle, shown for resolution 300 000 (Col. 2) and 70 000 (Col. 3). Column 4 and 5: synthetic CO Δv = 3 (5–2 P30) line profiles, from the original model, over a full cycle, shown for resolution 300 000 (Col. 4) and 70 000 (Col. 5).

Open with DEXTER
In the text
thumbnail Fig. 10

Radial velocities (RVs), with circles showing RVs derived from synthetic lines for two models, observation as crosses, and the grey line a running mean for the observations. The observed measurements are a compilation of RVs derived from FTS spectra for a large sample of different Mira stars, with the RVs converted to systematic velocities using the centre of mass radial velocity (CMRV) for each object (data from Lebzelter & Hinkle 2002).

Open with DEXTER
In the text
thumbnail Fig. 11

Comparison between synthetic RV curves of models and the running mean of the observations, seen in Fig. 10. The colour indicates the mean distance between the maximum and the minimum of that model’s RV curve and the running mean, described in Eq. (2). The grey squares are models that were previously discarded as they did not reproduce a realistic wind velocity and mass-loss rate combination.

Open with DEXTER
In the text
thumbnail Fig. 12

Comparison between the 1D pulsation models, from Bessell et al. (1996), and the boundary conditions described in Sect. 2.1.1 and Appendix A, defined by Δφs and Δφp. The colours show the relative least-squares error (L2rel), and the points indicates our calculated models and lines show Δφtot = Δφs + Δφp. Left and middle: model Z and model D from Bessell et al. (1996), which are fundamental mode pulsator models of o Ceti. Right: model E from Bessell et al. (1996), a first-overtone pulsator model of o Ceti.

Open with DEXTER
In the text
thumbnail Fig. A.1

Exact w = fφs), from the Fourier series Eq. (A.5), calculated with 10 000 terms and the polynomial fit from Eq. (A.6), using the terms from Table A.1. The blue circles represent the values of Δφs used in this work.

Open with DEXTER
In the text
thumbnail Fig. A.2

Number of Fourier terms to reach an L2 value of less than 0.1%.

Open with DEXTER
In the text

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

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

Initial download of the metrics may take a while.