© ESO, 2017

1. Introduction

Winds of asymptotic giant branch (AGB) stars are commonly assumed to be driven by radiative acceleration of dust grains. The pulsation of these stars give rise to shock waves, which propagate through the atmosphere and intermittently lift the upper layers to a few stellar radii. As the compressed gas cools down while moving away from the stellar photosphere, dust formation can start. The dust particles experience radiation pressure from stellar photons, thereby accelerating these grains outwards. Momentum from these accelerated grains is transferred to the surrounding gas through collisions, and thus an outflow of gas and dust from the stellar photosphere is triggered.

For AGB stars with a carbon-rich atmospheric chemistry (C/O > 1) there is a broad consensus that the observed outflows are driven by amorphous carbon grains. For M-type AGB stars (C/O < 1), however, the nature of the wind-driving dust species has been a matter of intense debate. Magnesium-iron silicates (olivine- and pyroxene-type materials) seem to be good candidates, considering the prominence of silicate features in mid-IR spectra and the comparatively high abundances of the constituent elements (Si, Mg, Fe, and O). Detailed models, however, show that silicate grains in the close vicinity of an AGB star, in the layers where the wind originates, can only have a very low content of Fe to avoid destruction by radiative heating (Woitke 2006). This leads to low levels of absorption at visual and near-IR wavelengths, corresponding to a radiative pressure that is insufficient to drive an outflow.

Using detailed atmosphere and wind models, Höfner (2008), Bladh et al. (2013, 2015), and Höfner et al. (2016) have therefore investigated an alternative source of radiation pressure, i.e. scattering of stellar photons on Fe-free silicate grains. This wind-driving mechanism requires grain radii of about 0.1–1 micron to make the dust particles efficient at scattering radiation in the near-IR wavelength region where the stellar flux peaks.

Determining dust particle sizes in the atmospheres and wind formation regions of AGB stars is not a straightforward task. At visual and near-IR wavelengths, where there is a diagnostic potential for the expected grain sizes, the central star is significantly brighter than the dust shell, which – in the best-case scenarios – is separated by only a few tens of milliarcseconds from the stellar surface. However, the brightness of dust shells in the visual and near-IR region is dominated by scattered stellar light, which has a high degree of linear polarization. Therefore, by observing in linearly polarized light the photons scattered by the dust shell can be distinguished from the direct stellar radiation.

Observations by Norris et al. (2012a) provided strong indications that dust particles with sizes where photon scattering can drive an outflow actually exist around AGB stars. These observations used aperture masking polarimetric interferometry (Norris et al. 2012b) on three nearby AGB stars (W Hya, R Dor, and R Leo). A linearly polarized signal was detected approximately at a distance of 2 stellar radii from the centre of the star, which was interpreted as a geometrically thin layer of dust particles scattering light towards the observer. Fitting the observed wavelength-dependent fraction of scattered light to the scattering efficiency of dust grains of various sizes gave a best fit for grains with radii around 0.3 μm.

More recently, resolved direct images of polarized light in the close vicinity of the AGB stars W Hya and R Dor were obtained with the Very Large Telescope (VLT), using the Spectro-Polarimetric High-contrast Exoplanet REsearch (SPHERE) – Zurich imaging polarimeter (ZIMPOL). Ohnaka et al. (2016) have detected three clumpy dust clouds at about 2–3 stellar radii around W Hya, which are comparable to structures resulting from large-scale convection in the 3D star-in-a-box models of Freytag & Höfner (2008). The analysis of the maps of the observed polarized intensity has indicated that they represent an optically thin medium with density enhancements of about a factor of 4 in the clouds and grain radii of 0.4–0.5 microns. A second-epoch measurement, which was taken about eight months later at minimum light (Ohnaka et al. 2017), revealed changes in the morphology of the clumps, and in grain size (showing a dominant grain population with radii of about 0.1 microns), compared to the first observation that was taken at pre-maximum light. These changes are consistent with an interpretation of grains forming around the luminosity minimum and growing during the phases of increasing luminosity; this also fits well with detailed dynamical atmosphere and wind models. Khouri et al. (2016) have obtained high-angular-resolution images of the AGB star R Dor at two epochs, demonstrating that the size and morphology of surface patterns change significantly in less than two months; this is compatible with convection-induced structures found in 3D models. These authors have interpreted these changes as variability in the excitation and/or density of TiO. They also observe light scattered by dust grains close to the star and a strong decrease in the dust density further out, which can be interpreted as rapidly falling densities due to wind acceleration or a change in mass loss rate.

In this paper we explore interpretations of polarized signatures from AGB stars from a theoretical point of view. In contrast to the simple, empirical structures used in the observational papers discussed above, our results are based on structures of the extended atmosphere and dusty circumstellar envelope, which are obtained with detailed, time-dependent dynamical atmosphere and wind models. Polarimetric signatures from the gas and dust shells are modelled and observables quantities that provide diagnostics of dust grain sizes are discussed.

2. Modelling: methods and parameters

2.1. Wind models

The atmospheric structures, for which scattered light is modelled, are snapshots taken from the dynamical models of M-type AGB stars presented in Höfner (2008) and Bladh et al. (2013). In these detailed models of the atmosphere and wind of AGB stars the mass loss is driven by a pulsation-enhanced dust-driven wind, where atmospheric levitation by pulsation-induced shock waves is followed by radiative acceleration of newly formed dust grains.

The dynamical models, produced with the DARWIN code (Höfner et al. 2016), are spherically symmetric with an inner boundary situated just below the photosphere and an outer boundary at about 20–30 stellar radii for models that develop a wind. The parameters that define the hydrostatic initial models are stellar mass M_⋆, stellar luminosity L_⋆, and the effective temperature of the star T_⋆ along with elemental abundances. The variability of the star is simulated by a sinusoidal variation of velocity and luminosity at the inner boundary of the model, characterized by a pulsation period P and the amplitudes of velocity and luminosity variations (parameters u_p and f_L; see Höfner et al. 2016, for details). The dynamical structure of the atmosphere, including the formation of shock waves and the development of stellar winds, is produced by simultaneously solving the equations of hydrodynamics, frequency-dependent radiative transfer, and grain growth. The abundances of atoms and molecules in the outer cool layers of the atmosphere are given by chemical equilibrium and the models include frequency-dependent opacities of both gas and dust in the radiative transfer to achieve realistic density and temperature structures. Dust formation is modelled with a time-dependent description for growth and evaporation of Fe-free silicate (Mg₂SiO₄) grains. As there is currently no well-developed nucleation theory for oxygen-rich atmospheres, the efficiency of grain formation is modelled by setting the number of seed particles per hydrogen atom (n_gr/n_H) as an input parameter. These seed particles start to grow when the thermodynamic conditions favour grain growth. As a consequence, the models produce a uniform grain size for all dust particles at a given distance from the stellar surface, at a given time.

As a direct result, the DARWIN models produce time series of snapshots of the radial structure of the atmosphere and wind. The models thereby provide information about properties such as velocity, temperature, density, and dust grain size as a function of radial distance and time. Selected snapshots are further processed in an a posteriori radiative transfer calculation to produce spectra and photometry, using the opacities from the COMA code (Aringer 2000; Aringer et al. 2009). Detailed information about the partial pressures of individual atoms and molecules (assuming chemical equilibrium) as a function of radial distance is also produced in this calculation. Three different atmosphere models are used in this paper, where the seed particle abundances (n_gr/n_H) differs and all other input parameters are similar. This generates dust shells with differently sized dust grains; see Table 1.

2.2. Polarimetric radiative transfer

In this section we describe the modelling of scattered light in the dust-enriched extended atmosphere of an AGB star. We calculate the state of polarization of scattered light as a function of wavelength and position within the atmosphere, which we then use to synthesize observations in polarized light with high spatial resolution.

2.2.1. Scattering particles

In order to fully include polarization in the calculations, we use the four Stokes parameters to describe intensity, $\mathit{I}\mathrm{=}\left[\begin{array}{c}\\ \mathit{i}\\ \mathit{q}\\ \mathit{u}\\ \mathit{v}\end{array}\right]\mathit{.}$(1)All scattering processes are described by multiplication of I with a 4 by 4 scattering matrix (F). This matrix has the property that multiplication with incoming intensity and single scattering albedo (a) gives the intensity (including state of polarization) of the scattered light. Scattering matrices are typically functions of the scattering angle θ, defined as the angle between the direction of propagation of the light before and after the scattering event, ${\mathit{I}}_{\mathrm{sca}}\mathrm{=}\mathit{aF}\mathrm{\left(}\mathit{\theta }\mathrm{\right)}{\mathit{I}}_{\mathrm{in}}\mathit{.}$(2)The first step in modelling scattered light is to create scattering matrices for particles in the atmosphere. We divide particles into two categories: small particles such as atoms and molecules and large particles such as dust grains. Small particles (for which we use the prefix m) are significantly smaller than the wavelength of the light. For these we use Rayleigh scattering theory to describe the scattering process. Large particles (prefix d) are comparable in size to the wavelength of the light and Mie scattering theory is used to describe scattered light in that case.

Rayleigh and Mie scattering matrices are calculated as given by Hansen & Travis (1974). For Mie Scattering we use the light scattering code BHMIE, based on Bohren & Huffman (1983). To ensure scattering matrices are correctly computed, we compared properties of our scattering matrices with the same properties generated from the light scattering code MiePlot. Both codes produced identical results. A comparison is shown in Fig. 9.

2.2.2. Atmospheric shells

The extended atmosphere is discretized into 300 homogeneous spherically symmetric shells of equal thickness. Within each shell the properties of gas and dust are known and assumed to be constant. We set the boundary of the innermost shell at 1.2 stellar radii (R_⋆), as there are no dust particles closer to the star than this and the outer boundary at 10 R_⋆, since outside of this region the density is too low to produce a significant amount of scattered light.

We combine the scattering matrices of individual species or sizes of particles into scattering matrices for entire shells, using the methods described by de Haan et al. (1987). For a homogeneous atmosphere, taking only single scattering into account, this can be carried out using $\begin{array}{ccc}& & \mathit{F}\mathrm{\left(}\mathit{\theta ,Z,\lambda }\mathrm{\right)}\mathrm{=}\\ & & \frac{{\sum }_{\mathit{k}}[{\mathit{b}}_{\mathrm{sca}}^{\mathit{m,k}}\mathrm{\left(}\mathit{Z,\lambda }\mathrm{\right)}\mathrm{·}{\mathit{F}}^{\mathit{m,k}}\mathrm{\left(}\mathit{\theta ,\delta }\mathrm{\right(}\mathit{k}\mathrm{\left)}\mathrm{\right)}]\mathrm{+}{\sum }_{\mathit{l}}[{\mathit{b}}_{\mathrm{sca}}^{\mathit{d,l}}\mathrm{\left(}\mathit{Z,\lambda }\mathrm{\right)}\mathrm{·}{\mathit{F}}^{\mathit{d,l}}\mathrm{\left(}\mathit{\theta ,x,y}\mathrm{\right)}]}{{\sum }_{\mathit{k}}{\mathit{b}}_{\mathrm{sca}}^{\mathit{m,k}}\mathrm{\left(}\mathit{Z,\lambda }\mathrm{\right)}\mathrm{+}{\sum }_{\mathit{k}}{\mathit{b}}_{\mathrm{sca}}^{\mathit{d,l}}\mathrm{\left(}\mathit{Z,\lambda }\mathrm{\right)}}\end{array}$(3)as described by Stam et al. (2000). Here Z denotes the index of the atmospheric shell (corresponding to distance from star), b the single scattering optical thickness, and F the scattering matrices (F^m for Rayleigh and F^d for Mie scattering). Summation is carried out over the different molecular species (k) and sizes of dust grains (l) found in the given shell.

The single scattering optical thickness for molecules (${\mathit{b}}_{\mathrm{sca}}^{\mathit{m}}$) is given by ${\mathit{b}}_{\mathrm{sca}}^{\mathit{m}}\mathrm{\left(}\mathit{Z,\lambda }\mathrm{\right)}\mathrm{=}\mathit{\eta }\mathrm{\left(}\mathit{Z}\mathrm{\right)}\frac{\mathrm{24}{\mathit{\pi }}^{\mathrm{3}}{\mathit{N}}_{\mathrm{av}}\mathrm{\left(}{\mathit{n}}^{\mathrm{2}}\mathrm{\right(}\mathit{\lambda }\mathrm{)}\mathrm{-}\mathrm{1}{\mathrm{)}}^{\mathrm{2}}\mathrm{(}\mathrm{6}\mathrm{+}\mathrm{3}\mathit{\delta }\mathrm{)}}{{\mathit{\lambda }}^{\mathrm{4}}{\mathit{N}}_{\mathrm{L}}^{\mathrm{2}}\mathit{R}\mathrm{\left(}{\mathit{n}}^{\mathrm{2}}\mathrm{\right(}\mathit{\lambda }\mathrm{)}\mathrm{+}\mathrm{2}{\mathrm{)}}^{\mathrm{2}}\mathrm{(}\mathrm{6}\mathrm{-}\mathrm{7}\mathit{\delta }\mathrm{)}}\underset{{\mathit{Z}}_{\mathit{i}}}{\overset{{\mathit{Z}}_{\mathit{j}}}{\mathrm{\int }}}\frac{\mathit{p}\mathrm{\left(}\mathit{z}\mathrm{\right)}}{\mathit{T}\mathrm{\left(}\mathit{z}\mathrm{\right)}}\mathrm{d}\mathit{z,}$(4)where N_av is the Avagado constant, N_L is the Loschmidt constant, R is the gas constant, n is the refraction index, δ is the depolarization factor, p(z) is the pressure as a function of distance from the star, and T(z) is temperature as a function of distance from the star. To incorporate the number density of different atoms and molecules, the expression is multiplied with the volume mixing ratio η(Z) of the atomic or molecular species (assumed to be constant within the shell). The integration over height is carried out assuming exponential decrease in P/T through the shell, where P(z_i), T(z_i), P(z_j), and T(z_j) correspond to tabulated values at the top and bottom of the shell.

The single scattering optical thickness for dust grains (${\mathit{b}}_{\mathrm{sca}}^{\mathit{d}}$) are given by ${\mathit{b}}_{\mathrm{sca}}^{\mathit{d}}\mathrm{\left(}\mathit{Z,\lambda }\mathrm{\right)}\mathrm{=}\mathit{N}\mathrm{\left(}\mathit{r,Z}\mathrm{\right)}\mathrm{·}\mathrm{(}{\mathit{z}}_{\mathit{i}}\mathrm{-}{\mathit{z}}_{\mathit{j}}\mathrm{)}\mathrm{·}{\mathit{\sigma }}_{\mathrm{sca}}\mathrm{\left(}\mathit{r,\lambda }\mathrm{\right)}\mathit{.}$(5)Here N(r,Z) is the number density of particles with radius r in shell Z, which is assumed to be constant within the shell; σ_sca is the scattering cross section of particles with radius r, calculated using the BHMIE code.

With both scattering and absorption optical thickness of molecules and dust grains, see Eqs. (6) and (7), we can compute the combined optical thickness of an atmospheric shell (Eq. (8)), ${\mathit{b}}_{\mathrm{abs}}^{\mathit{m}}\mathrm{\left(}\mathit{Z,\lambda }\mathrm{\right)}\mathrm{=}\mathit{\eta }\mathrm{\left(}\mathit{Z}\mathrm{\right)}\frac{{\mathit{N}}_{\mathrm{av}}{\mathit{\sigma }}_{\mathrm{abs}}^{\mathit{m}}\mathrm{\left(}\mathit{Z,\lambda }\mathrm{\right)}}{\mathit{R}}\underset{{\mathit{Z}}_{\mathit{i}}}{\overset{{\mathit{Z}}_{\mathit{j}}}{\mathrm{\int }}}\frac{\mathit{p}\mathrm{\left(}\mathit{z}\mathrm{\right)}}{\mathit{T}\mathrm{\left(}\mathit{z}\mathrm{\right)}}\mathrm{d}\mathit{z}\mathit{.}$(6)The integral is evaluated assuming exponential decrease through the shell. The molecular absorption cross section σ_abs is assumed to be constant within the shell. ${\mathit{b}}_{\mathrm{abs}}^{\mathit{d}}\mathrm{\left(}\mathit{Z,\lambda }\mathrm{\right)}\mathrm{=}\mathit{N}\mathrm{\left(}\mathit{r,Z}\mathrm{\right)}\mathrm{·}\mathrm{(}{\mathit{z}}_{\mathit{i}}\mathrm{-}{\mathit{z}}_{\mathit{j}}\mathrm{)}\mathrm{·}{\mathit{\sigma }}_{\mathrm{abs}}\mathrm{\left(}\mathit{r,\lambda }\mathrm{\right)}$(7)$\begin{array}{ccc}& & \mathit{b}\mathrm{\left(}\mathit{Z,\lambda }\mathrm{\right)}\mathrm{=}\\ & & \sum _{\mathit{k}}{\mathit{b}}_{\mathrm{sca}}^{\mathit{m,k}}\mathrm{\left(}\mathit{Z,\lambda }\mathrm{\right)}\mathrm{+}\sum _{\mathit{k}}{\mathit{b}}_{\mathrm{abs}}^{\mathit{m,k}}\mathrm{\left(}\mathit{Z,\lambda }\mathrm{\right)}\mathrm{+}\sum _{\mathit{l}}{\mathit{b}}_{\mathrm{sca}}^{\mathit{d,l}}\mathrm{\left(}\mathit{Z,\lambda }\mathrm{\right)}\mathrm{+}\sum _{\mathit{l}}{\mathit{b}}_{\mathrm{abs}}^{\mathit{d,l}}\mathrm{\left(}\mathit{Z,\lambda }\mathrm{\right)}\mathit{.}\end{array}$(8)

2.2.3. Radiative transfer

With scattering matrices for all atmospheric shells, we can now perform the radiative transfer. In the first step we consider each shell individually. For each shell, we create reflection and transmission matrices (de Haan et al. 1987). These matrices have the property that when multiplied with the incoming light arriving at the bottom of the atmosphere, the product is the light emerging from the upper (reflection matrix) or lower (transmission matrix) boundary of the layer. Thus, given the incident light on a shell and its angle to the local zenith, we can calculate the intensity (including state of polarization) of the outgoing light at a given angle to local zenith. The separation into reflection and transmission matrix is needed owing to the different paths light has to take to reach either upper or lower boundary of the shell; see Fig. 1.

Fig. 1 Possible paths for light moving through an atmospheric shell. The transmission matrix describes light emerging at the upper boundary of this shell and the reflection matrix describes light emerging from the lower boundary of this shell.

For transmission and reflection matrices, the direction of propagation of the incoming and outgoing light is needed. This is not uniquely defined by scattering angle, and we therefore introduce μ and ϕ to define direction of propagation; μ = | cosΘ | is the angle between the direction of propagation of the light and the local zenith, where μ_o is for the outgoing light and μ_i is for incoming light. The azimuth angle ϕ_o (outgoing) and ϕ_i (incoming) are measured from an arbitrary plane containing the local zenith and clockwise when looking downwards; see Fig. 2. Because of rotational symmetry only the difference in azimuth angles (ϕ = ϕ_i−ϕ_o) is relevant.

Fig. 2 Direction of propagation of incoming and outgoing light is defined by μ and ϕ.

When multiple scattering is neglected and no internal light sources are present inside the shell, the intensity of light at optical depth τ in an atmospheric shell can be calculated using $\begin{array}{ccc}& & {\mathit{I}}_{\mathrm{down}}\mathrm{\left(}\mathit{\tau ,}{\mathit{\mu }}_{\mathrm{o}}\mathit{,\varphi ,\lambda }\mathrm{\right)}\mathrm{=}(\frac{\mathrm{1}}{\mathit{\pi }}\underset{\mathrm{0}}{\overset{\mathrm{1}}{\mathrm{\int }}}{\mathit{\mu }}_{\mathit{i}}\mathrm{d}{\mathit{\mu }}_{\mathit{i}}\underset{\mathrm{0}}{\overset{\mathrm{2}\mathit{\pi }}{\mathrm{\int }}}\mathrm{d}\mathit{\varphi }\mathrm{·}[\frac{\mathrm{1}}{\mathrm{4}{\mathit{\mu }}_{\mathrm{o}}{\mathit{\mu }}_{\mathit{i}}}\underset{\mathrm{0}}{\overset{\mathit{\tau }}{\mathrm{\int }}}\mathrm{d}{\mathit{\tau }}^{\mathrm{\prime }}\\ & & {\mathrm{e}}^{\frac{\mathrm{-}{\mathit{\tau }}^{\mathrm{\prime }}}{{\mathit{\mu }}_{\mathit{i}}}}{\mathrm{e}}^{\frac{\mathrm{-}\mathrm{(}\mathit{\tau }\mathrm{-}{\mathit{\tau }}^{\mathrm{\prime }}\mathrm{)}}{{\mathit{\mu }}_{\mathrm{o}}}}\mathit{a}\mathrm{\left(}{\mathit{\tau }}^{\mathrm{\prime }}\mathrm{\right)}\mathrm{\times }\mathit{F}\mathrm{\left(}\mathit{\theta }\mathrm{\right(}{\mathit{\mu }}_{\mathit{i}}\mathit{,}{\mathit{\mu }}_{\mathrm{o}}\mathit{,\varphi }\mathrm{\left)}\mathit{,\lambda }\mathrm{\right)}]\mathrm{+}{\mathrm{e}}^{\frac{\mathrm{-}\mathit{\tau }}{{\mathit{\mu }}_{\mathrm{o}}}})\mathrm{·}{\mathit{I}}_{\mathit{i}}\mathrm{\left(}\mathrm{0}\mathit{,}{\mathit{\mu }}_{\mathit{i}}\mathit{,\varphi ,\lambda }\mathrm{\right)}\end{array}$(9)$\begin{array}{ccc}& & {\mathit{I}}_{\mathrm{up}}\mathrm{(}\mathit{\tau ,}\mathrm{-}{\mathit{\mu }}_{\mathrm{o}}\mathit{,\varphi ,\lambda }\mathrm{)}\mathrm{=}(\frac{\mathrm{1}}{\mathit{\pi }}\underset{\mathrm{0}}{\overset{\mathrm{1}}{\mathrm{\int }}}{\mathit{\mu }}_{\mathit{i}}\mathrm{d}{\mathit{\mu }}_{\mathit{i}}\underset{\mathrm{0}}{\overset{\mathrm{2}\mathit{\pi }}{\mathrm{\int }}}\mathrm{d}\mathit{\varphi }\mathrm{·}[\frac{\mathrm{1}}{\mathrm{4}{\mathit{\mu }}_{\mathrm{o}}{\mathit{\mu }}_{\mathit{i}}}\underset{\mathit{\tau }}{\overset{\mathit{b}}{\mathrm{\int }}}\mathrm{d}{\mathit{\tau }}^{\mathrm{\prime }}\\ & & {\mathrm{e}}^{\frac{\mathrm{-}{\mathit{\tau }}^{\mathrm{\prime }}}{{\mathit{\mu }}_{\mathit{i}}}}{\mathrm{e}}^{\frac{\mathrm{-}\mathrm{(}{\mathit{\tau }}^{\mathrm{\prime }}\mathrm{-}\mathit{\tau }\mathrm{)}}{{\mathit{\mu }}_{\mathrm{o}}}}\mathit{a}\mathrm{\left(}{\mathit{\tau }}^{\mathrm{\prime }}\mathrm{\right)}\mathrm{\times }\mathit{F}\mathrm{\left(}\mathit{\theta }\mathrm{\right(}{\mathit{\mu }}_{\mathit{i}}\mathit{,}\mathrm{-}{\mathit{\mu }}_{\mathrm{o}}\mathit{,\varphi }\mathrm{\left)}\mathit{,\lambda }\mathrm{\right)}])\mathrm{·}{\mathit{I}}_{\mathit{i}}\mathrm{\left(}\mathit{\tau ,}{\mathit{\mu }}_{\mathit{i}}\mathit{,\varphi ,\lambda }\mathrm{\right)}\mathit{.}\end{array}$(10)Here I_down is the intensity of light moving downwards at optical depth τ and at an angle given by μ_o and ϕ, and likewise I_up is the intensity of light moving upwards at a direction given by μ_o and ϕ. The integration over τ′ is to account for reflections along the depth of the shell adding to the intensity, which for light moving downwards is from 0 (top of the layer) to τ, and for light moving upwards from b (bottom of layer) to τ. The exponentials give the attenuation when the light moves through the atmosphere at the angle given by μ_i (before scattering) or μ_o (after scattering).

The integration within the square brackets can be evaluated analytically for a homogeneous atmospheric shell; a is considered to be constant within the shell. We call these matrices reflection matrix R when light entering the boundary of the layer is scattering back towards the same boundary (0° ≤ Θ_o< 90°) and transmission matrix T when light is emerging at the opposite boundary (90° ≤ Θ_o ≤ 180°). Integration over τ is carried out from 0 to b in both cases, i.e. we calculate the intensity of light at the top or bottom of the shell for light moving upwards and downwards, respectively, as follows:$\mathit{R}\mathrm{\left(}{\mathit{\mu }}_{\mathit{i}}\mathit{,}{\mathit{\mu }}_{\mathrm{o}}\mathit{,\varphi ,\lambda }\mathrm{\right)}\mathrm{=}\frac{\mathit{a}}{\mathrm{4}\mathrm{(}{\mathit{\mu }}_{\mathrm{o}}\mathrm{+}{\mathit{\mu }}_{\mathit{i}}\mathrm{)}}(\mathrm{1}\mathrm{-}{\mathrm{e}}^{\frac{\mathrm{-}\mathit{b}}{{\mathit{\mu }}_{\mathrm{o}}}\mathrm{+}\frac{\mathrm{-}\mathit{b}}{{\mathit{\mu }}_{\mathit{i}}}})\mathrm{\times }\mathit{F}\mathrm{\left(}\mathit{\theta }\mathrm{\right(}{\mathit{\mu }}_{\mathit{i}}\mathit{,}{\mathit{\mu }}_{\mathrm{o}}\mathit{,\varphi }\mathrm{\left)}\mathit{,\lambda }\mathrm{\right)}\mathrm{;}$(11)if μ_o ≠ μ_i, $\mathit{T}\mathrm{\left(}{\mathit{\mu }}_{\mathit{i}}\mathit{,}{\mathit{\mu }}_{\mathrm{o}}\mathit{,\varphi ,\lambda }\mathrm{\right)}\mathrm{=}\frac{\mathit{a}}{\mathrm{4}\mathrm{(}{\mathit{\mu }}_{\mathrm{o}}\mathrm{-}{\mathit{\mu }}_{\mathit{i}}\mathrm{)}}({\mathrm{e}}^{\frac{\mathrm{-}\mathit{b}}{{\mathit{\mu }}_{\mathrm{o}}}}\mathrm{-}{\mathrm{e}}^{\frac{\mathrm{-}\mathit{b}}{{\mathit{\mu }}_{\mathit{i}}}})\mathrm{\times }\mathit{F}\mathrm{\left(}\mathit{\theta }\mathrm{\right(}{\mathit{\mu }}_{\mathrm{0}}\mathit{,\mu ,\varphi }\mathrm{\left)}\mathit{,\lambda }\mathrm{\right)}\mathrm{;}$(12)if μ_o = μ_i, $\mathit{T}\mathrm{\left(}{\mathit{\mu }}_{\mathit{i}}\mathit{,}{\mathit{\mu }}_{\mathrm{o}}\mathit{,\varphi ,\lambda }\mathrm{\right)}\mathrm{=}\frac{\mathit{a}}{\mathrm{4}{\mathit{\mu }}_{\mathrm{o}}^{\mathrm{2}}}{\mathrm{e}}^{\frac{\mathrm{-}\mathit{b}}{{\mathit{\mu }}_{\mathrm{o}}}}\mathrm{\times }\mathit{F}\mathrm{\left(}\mathit{\theta }\mathrm{\right(}{\mathit{\mu }}_{\mathit{i}}\mathit{,}{\mathit{\mu }}_{\mathrm{o}}\mathit{,\varphi }\mathrm{\left)}\mathit{,\lambda }\mathrm{\right)}\mathit{.}$(13)With transmission and reflection matrices for all shells and all combinations of angles for incoming and outgoing light, we can now perform radiative transfer calculations for the entire extended atmosphere. We consider only a plane of the atmosphere (the 3D structure is reintroduced in Sect. 2.2.4), and calculate the amount of scattered light from each point on this plane for all shells. When calculating the amount of scattered light reaching the observer from some point in the atmosphere, there are two distinct cases, A and B in Fig. 3. For forwards scattered light (case A), we first trace stellar light through the inner shells (using their transmission matrices where incoming and outgoing angles are the same) to the scattering point. We then compute the intensity of light scattered in a direction towards the observer using the transmission matrix of this shell and the relevant incoming and outgoing angles. Light is then traced through the outer layers as it travels towards the outer boundary of our model using the transmission matrices of the shell, where incoming and outgoing angles are the same. For backscattered light (case B), we trace stellar light through the inner shells to the scattering point and compute the intensity of light scattered in a direction towards the observer using the reflection matrix of this shell. Light is then traced through first the inner and then the outer shells as it travels towards the observer. Light that is absorbed or scattered in directions not towards the observer is assumed to be lost, i.e. no multiple scattering is included, which implies that the results are only valid for optically thin atmospheres. Using this approach should be acceptable when modelling typical M-type ABG stars. Their spectra shows very little absorption from dust in the visual and near-infrared (near-IR) wavelength regions and very little circumstellar reddening, resulting in an optically thin atmosphere.

Fig. 3 Sketch of the scheme for radiative transfer calculations. Stellar light moves outwards from the star to a point at which some of the light is scattered towards a distant observer. From this point light travels through the outer (case A) or first through the inner and then the outer (case B) before it reaches the boundary of our model.

2.2.4. Synthesizing observations

These calculations generate the scattered light, including state of polarization, from each point within a plane of the shell. From this data a 1D image is constructed by combining scattered light along the same line of sight, but emerging from different shells. Since the extended atmosphere is assumed to be spherically symmetric, we can use this to construct a 2D image by postulating that the intensity is a function of distance from the centre of the star. We remove the intensity emerging from closer than 1.2 R_⋆ (projected distance as seen by a distant observer); i.e. we remove the direct stellar light.

This radiative transfer scheme was tested on a simple extended atmosphere with no gas and only a single thin layer of dust grains situated at 2.5 R_⋆; these dust grains all have a grain radius of 0.3 μm, which is similar to the set-up used by Norris et al. (2012a). The results from this test is shown in Fig. 4. As expected, when the circumstellar dust shell is viewed in linearly polarized light one can observe bright arcs at the distance where dust particles are present. The intensity outside these arcs is zero and the intensity inside the arcs is decreasing towards the star; the intensity from the inner parts comes from backwards and forwards scattered light.

Fig. 4 Synthesized image in linearly polarized light (λ = 2.3 μm) of a simplified dust shell model (a single thin layer containing dust particles with a grain radius of 0.3 μm and no gas). Blue represents the horizontally polarized light. Red represents the vertically polarized light. Just as expected, we observe bright arcs of polarized light where dust grains are situated.

3. Results

3.1. Synthetic polarized images of detailed models

The radiative transfer method described in Sect. 2.2 is applied to snapshots of detailed wind models taken at different pulsation phases. We use model B2 for in-depth studies (see Table 1). Synthetic spectra for this model at four different phases are shown in Fig. 5. The corresponding density and grain-size profiles, as a function of radial distance from the star, as well as the flow velocities, are shown in Fig. 6. Images in linearly polarized light were synthesized for the same four phases and in three different wavelength regions (1.2 μm, 1.5 μm, and 2.1 μm). Results are shown in Fig. 7.

Fig. 5 Spectral energy distribution for the model B2 at four different pulsation phases.

Fig. 6 Density structure of the dust component as a function of radial distance for the model B2, taken at four different pulsation phases (same as in Fig. 5). Colour-coded solid lines show density and dust grain size. Dashed grey lines show the flow velocity.

Fig. 7 Linearly polarized light from the spherically symmetric model B2. Colour coding the same as in Fig. 4. Scattered light computed at the maximum phase (0.0), minimum phase (0.5), and intermediate phases (0.25 and 0.75) is shown.

The overall structure of polarized light from model B2 is similar to the simplified model, shown in Fig. 4, with bright arcs where the density of the dust is high, and decreasing intensity closer to the star. The most notable difference between Figs. 7 and 4 is the non-zero intensity outside the bright arcs in Fig. 7. This is due to the more extended distribution of dust in the detailed model compared to the geometrically thin dust shell in the simple test case. Clear variations in the polarized light with pulsation phase can also be seen in Fig. 7, both regarding the overall intensity and the radial position of the arcs.

This variation in polarized light is caused by temporal changes in grain size and density structure during the pulsation cycle, as illustrated in Fig. 6. A new dust layer starts to condense at the inner edge of the existing dust shell around the luminosity minimum (panel 3, φ = 0.5) when the overall temperature of the atmosphere is lowest. In panel three of Fig. 6, a small peak in dust density appears with grain sizes around 0.15 μm. During the consecutive phases the newly formed dust layer moves outwards and the grains grow, reaching a size of 0.45 μm near the luminosity maximum (see panel 1, φ = 0.00, which follows panel 4 because of the cyclic behaviour of the model). The change in the radial position of the density peak is reflected in the location of the arcs in Fig. 7.

The top left panel of Fig. 7 (phase 0 at 1.2 μm) shows another interesting phenomenon, namely a dependence of polarization state on the distance from the star. This is caused by the fact that light seen at different distances has been scattered at different angles and that the polarization state depends on the scattering angle. The angle, and thereby the distance, at which this change in polarization state occurs depends on the size of the dust grains relative to the wavelength of the light. This polarization flip can be seen in Fig. 7 for λ = 1.2 μm as a change in colour from red to blue (or vice versa) with radial distance. Light scattered at all angles has the same state of polarization at longer wavelengths (λ = 1.5 μm and λ = 2.1 μm). In principle, the detection of arcs with a flip in their polarization state could be used as a way to determine grain sizes, although it requires very good spatial resolution.

3.2. Simulating non-sphericity

Up to this point we only considered spherically symmetric models, assuming the same properties of the dust and gas shell in all directions. However, imaging of circumstellar dust close to AGB stars shows significant deviations from spherical symmetry. Recently, Ohnaka et al. (2016, 2017) and Khouri et al. (2016) have obtained direct images of polarized light around two nearby AGB stars (W Hya and R Dor) with SPHERE/ZIMPOL, demonstrating the variable, non-spherical nature of the circumstellar material. The images of both stars show clumpy structures, which vary on timescales of several months. In particular, Ohnaka et al. (2016, 2017) have interpreted the observed change in dominant grain size between pre-maximum and minimum phases of W Hya as an effect of grain growth during the phases of increasing luminosity, following the formation of new grains around the luminosity minimum. Both, the grain sizes, of about 0.1 and 0.5 micron at minimum and pre-maximum phases, respectively, and the timing of grain growth derived from these observations, fit very well with our model B2, as shown in Fig. 6.

To mimic a non-spherical structure, as in the observations, we combine a range of pulsation phases into a single synthetic image, using snapshots from model B2 as a radial structure in a certain direction, rather than a global, spherically symmetric structure. Two maximum phases are placed at opposite sides of the star, but slightly of-axis (at 0° and 200°). Intermediate phases are spread evenly between the maximum phases in clockwise direction with interpolation between phases to create a smooth transition. See left panel of Fig. 8 for a sketch of where phases are located. The middle and right panels of this figure show the polarized light for the same underlying structure but at different wavelengths. A difference compared to a spherically symmetric model is that the images now produce a net polarization.

Fig. 8 Non-spherically symmetric version of model B2 (composite of different phases) in polarized light. The upper panel represents the horizontally and vertically polarized light; the colour coding is the same as in Fig. 4. The lower panel represents the combined polarized intensity.

3.3. Dust size diagnostics

Grain size is a critical parameter, which determines the contribution of pure scattering to the radiative acceleration experienced by the dust, and therefore has important consequences for the wind mechanism. The state of polarization of scattered light depends on both the size parameter (i.e. the relative size of scattering particles to wavelength of light) and the scattering angle. Scattering matrix elements S₂₁ (shown in Fig. 9) gives the intensity and state of polarization of scattered light, where the negative values correspond to polarization perpendicular to direction of incoming light and positive values correspond to polarization parallel to incoming light, as function of scattering angle. As the size parameter changes, the state of polarization for a given scattering angle can change. For model B2, with dust grains with radius around 0.45 μm (solid black line in Fig. 9), S₂₁ changes from positive to negative at scattering angles around 90° when the wavelength increases.

To see how the state of polarization of scattered light changes with grain size and wavelength, we calculate the polarized light for a set of simplified models with a thin spherically symmetric dust shell situated at 2.5 R_⋆, consisting of dust particles of uniform size (see Fig. 10). Synthetic polarized light was computed for 10 different dust sizes and in wavelengths ranging from 1.0 to 2.5 μm. In the wavelength region close to where the state of polarization switches from vertical to horizontal, at scattering angles around 90° (i.e. from positive to negative in Fig. 9), one can observe multiple bright arcs of polarized light with alternating states of polarization. At shorter wavelengths than this threshold, the polarization is dominated by light polarized parallel to local stellar radiation, while at longer wavelengths, polarization is dominated by light polarized perpendicular to local stellar radiation. In all cases, the switch from one dominant state of polarization to the other happens at wavelengths around three times the radius of dust grains. Thus, if such change in state of polarization could be observed, the size of the scattering particles could easily be estimated.

To observe changes in state of the polarization, high spatial resolution, sufficient to resolve individual arcs of polarized light, would be needed. Alternately, this could also be observed in unresolved images if the star is asymmetrical in such way that a net polarized signal from the shell could be detected. One such example, based on model B2, was discussed in Sect. 3.2 and shown in Fig. 8. In a similar way, we compute images for non-spherical structures constructed from all three models in Table 1, in wavelengths from 1.0 to 2.5 μm; see Fig. 11. Along with a spatially resolved image of each model, the net state of polarization of the combined dust shell is indicated with an arrow, assuming no net polarization from the star itself. In all cases net polarization changes orientation at wavelength within ± 1 μm as predicted by simplified models; see Fig. 10.

Fig. 9 Matrix element S21 of Mie-scattering matrices for Mg2SiO4 dust grain with optical data from Jäger et al. (2003). The top panel indicates λ = 1.2 μm, middle panel indicates λ = 1.5 μm, and bottom panel indicates λ = 2.1 μm. The dashed lines show S21 for a distribution of dust grain sizes from 0.1 to 1.0 μm, and the solid lines show S21 for sizes corresponding to the average size in models. Grey dotted lines show the same quantity calculated for 0.45 μm sized forsterite grains using the light scattering code MiePlot.

Fig. 10 Linearly polarized light from a grid of simplified models. Colour coding is the same as in Fig. 4.

Fig. 11 Linearly polarized light from non-spherically symmetric detailed models with varying grain size and at different wavelengths. The dominating state of polarization is indicated by arrows at the top left corner of each model. Colour coding is the same as in Fig. 4. Row B3 shows ⟨ agr ⟩ = 0.33 μm, row B2 shows ⟨ agr ⟩ = 0.45 μm , and row B1 shows ⟨ agr ⟩ = 0.66 μm.

4. Discussion and conclusions

In this paper we have presented a tool for computing synthetic images of scattered light from detailed dynamical models of gas and dust shells around AGB stars, and we have discussed the diagnostic potential of polarization regarding dust grain sizes. The study presented here was originally inspired by the work of Norris et al. (2012a), who used aperture masking polarimetric interferometry at near-IR wavelengths to determine grain sizes in the close vicinity of W Hya, R Dor, and R Leo. In particular, our intention was to check how assumptions about the radial structure of the circumstellar material, for example regarding density and grain size, affect the resulting scattered light. Norris et al. (2012a) used the assumption of a geometrically thin, spherically symmetric shell of dust particles with a single grain size for interpreting their observed data. In contrast, our results are based on extended, time-dependent dynamic atmosphere and wind models, computed with the DARWIN code, which feature a detailed treatment of grain growth (Höfner 2008; Bladh et al. 2013; Höfner et al. 2016)

From the corresponding synthetic images the importance of accounting for temporal changes with pulsation phase becomes apparent. The intensity of the polarized light, the position of the inner edge of the dust shell (where the polarized intensity is highest), and the size of the dust grains near the inner edge are all changing with the luminosity phase, which significantly affects the resulting polarized signal. In addition to the temporal variations caused by radial pulsations, non-spherical structures in the dust shells, for example due to large-scale convective flows (Freytag & Höfner 2008), can also have an impact on the polarized light. In the present paper, we simulated this effect by combining a range of pulsation phases of a 1D model into a single 3D structure, before computing synthetic images. We therefore conclude that the interpretation of data using time-independent, spherically symmetric models should be made with caution.

While the complex, variable morphology of the circumstellar dust shell can make the interpretation of resolved data for nearby stars more difficult, it might turn out to be a blessing in disguise for unresolved objects. An asymmetry of the circumstellar envelope creates a net polarization of the total scattered light from the dust shell, which could give diagnostics of grain sizes. The ratio between the size of the scattering particles and the observed wavelength is the dominant parameter determining at what wavelengths net polarization switches from parallel to perpendicular (see Fig. 9). By determining at what wavelength the dominant state of polarization changes, average particle size could be constrained. It should, however, be noted that a non-detection of such a change in net polarization cannot be used to put constraints on grain sizes. Apart from grains not being within the corresponding size regime for the given wavelength range, there could be numerous explanations, for example, no significant deviations from spherical symmetry, superposition of grains with a large spread in sizes, polarized light dominated by gas from the inner parts of the shell, etc.

Finally, it should be mentioned that the chemical composition of the dust grains cannot be reliably determined from the scattered light observations at visible and near-IR wavelengths alone, but only in combination with measurements in the mid-IR regime where the solid particles show characteristic lattice modes. Considerations of thermal stability and radiative heating in the close vicinity of the star suggest Fe-free silicates (Mg₂SiO₄, MgSiO₃) and/or Al₂O₃ as likely candidates for the light-scattering dust grains. Spectro-interferometric observations of AGB stars indicate that such grains can exist around 2 stellar radii, or even below, with Al₂O₃ forming closer to the star than silicates (e.g. Wittkowski et al. 2007; Zhao-Geisler et al. 2011, 2012; Karovicova et al. 2013; Sacuto et al. 2013). This is consistent with an interpretation that the smaller grains found by Ohnaka et al. (2017) around W Hya near minimum light may consist of Al₂O₃ and that they grow silicate mantles at later phases. The scenario of winds driven by grains consisting of Al₂O₃ cores and Mg₂SiO₄ mantles was recently studied with DARWIN models by Höfner et al. (2016), resulting in both, realistic wind properties and variations of visual and near-IR colours in good agreement with observations. In this context, the Al₂O₃ particles prior to silicate condensation experience a radiative pressure that is insufficient for driving an outflow, and they form a gravitationally bound dust shell close to the star in accordance with the results of Ireland et al. (2005) and Khouri et al. (2015, 2016).

Acknowledgments

This work has been supported by the Swedish Research Council (Vetenskapsrådet) and by the ERC Consolidator Grant funding scheme (project project STARKEY, G.A. No. 615604). The computations of spectra and photometry were performed on resources provided by the Swedish National Infrastructure for Computing (SNIC) at UPPMAX.

References

All Tables

All Figures

	Fig. 1 Possible paths for light moving through an atmospheric shell. The transmission matrix describes light emerging at the upper boundary of this shell and the reflection matrix describes light emerging from the lower boundary of this shell.
In the text

	Fig. 2 Direction of propagation of incoming and outgoing light is defined by μ and ϕ.
In the text

	Fig. 3 Sketch of the scheme for radiative transfer calculations. Stellar light moves outwards from the star to a point at which some of the light is scattered towards a distant observer. From this point light travels through the outer (case A) or first through the inner and then the outer (case B) before it reaches the boundary of our model.
In the text

	Fig. 4 Synthesized image in linearly polarized light (λ = 2.3 μm) of a simplified dust shell model (a single thin layer containing dust particles with a grain radius of 0.3 μm and no gas). Blue represents the horizontally polarized light. Red represents the vertically polarized light. Just as expected, we observe bright arcs of polarized light where dust grains are situated.
In the text

	Fig. 5 Spectral energy distribution for the model B2 at four different pulsation phases.
In the text

	Fig. 6 Density structure of the dust component as a function of radial distance for the model B2, taken at four different pulsation phases (same as in Fig. 5). Colour-coded solid lines show density and dust grain size. Dashed grey lines show the flow velocity.
In the text

	Fig. 7 Linearly polarized light from the spherically symmetric model B2. Colour coding the same as in Fig. 4. Scattered light computed at the maximum phase (0.0), minimum phase (0.5), and intermediate phases (0.25 and 0.75) is shown.
In the text

	Fig. 8 Non-spherically symmetric version of model B2 (composite of different phases) in polarized light. The upper panel represents the horizontally and vertically polarized light; the colour coding is the same as in Fig. 4. The lower panel represents the combined polarized intensity.
In the text

Fig. 9 Matrix element S21 of Mie-scattering matrices for Mg2SiO4 dust grain with optical data from Jäger et al. (2003). The top panel indicates λ = 1.2 μm, middle panel indicates λ = 1.5 μm, and bottom panel indicates λ = 2.1 μm. The dashed lines show S21 for a distribution of dust grain sizes from 0.1 to 1.0 μm, and the solid lines show S21 for sizes corresponding to the average size in models. Grey dotted lines show the same quantity calculated for 0.45 μm sized forsterite grains using the light scattering code MiePlot.

In the text

	Fig. 10 Linearly polarized light from a grid of simplified models. Colour coding is the same as in Fig. 4.
In the text

	Fig. 11 Linearly polarized light from non-spherically symmetric detailed models with varying grain size and at different wavelengths. The dominating state of polarization is indicated by arrows at the top left corner of each model. Colour coding is the same as in Fig. 4. Row B3 shows ⟨ agr ⟩ = 0.33 μm, row B2 shows ⟨ agr ⟩ = 0.45 μm , and row B1 shows ⟨ agr ⟩ = 0.66 μm.
In the text