© The Authors 2025

1. Introduction

At the final stages of their lives, low- and intermediate-mass stars ascend the asymptotic giant branch (AGB). In this phase, stars consist of an inert core of carbon and oxygen surrounded by nuclear burning shells and a convective envelope (Herwig 2005). The stellar envelope expands dramatically, while pulsations and convection extend the atmosphere, creating conditions suitable for dust condensation. The current mass-loss paradigm consists of radiation pressure acting on dust grains to overcome the gravitational pull of the star and trigger a slow and dense wind (e.g. Höfner & Olofsson 2018).

However, oxygen-rich AGB stars present a particularly challenging case. The nucleation of oxygen-rich dust is thought to begin with materials such as aluminium oxide (Al₂O₃), followed by Fe-free silicate grains (MgSiO₃) (Karovicova et al. 2013; Höfner et al. 2016). A theoretical crisis emerged when it was demonstrated that Fe-free silicates - the dominant species inferred from observations - are insufficiently opaque to drive winds (Woitke 2006). Their work also demonstrated that while Fe-bearing silicates (MgFeSiO₄) could theoretically provide sufficient opacity, these grains would become too warm to exist in the inner envelope where wind acceleration begins. It was instead proposed that Fe-free silicate grains could drive winds if they grow to sizes of at least 0.3 μm, as larger grains effectively scatter stellar photons despite their low absorption cross-sections (Höfner 2008).

High-angular-resolution observations of polarised light have revolutionised our ability to test these models by directly studying light scattered by dust in the dust-formation and wind-acceleration zones (Norris et al. 2012; Khouri et al. 2016, 2020; Ohnaka et al. 2016). These observations confirm the presence of dust grains roughly in the expected size range around oxygen-rich AGB stars, but whether these observed dust populations can actually drive the winds under physically realistic conditions remains to be determined.

R Doradus (R Dor), the closest oxygen-rich AGB star to the Sun at a distance of ∼59 pc (Knapp et al. 2003; Khouri et al. 2024), provides an ideal laboratory for investigating this question. The combination of its proximity and large angular size (diameter ∼62.2 mas; Vlemmings et al. 2018) enables spatially resolved observations of the dust-formation and wind-acceleration regions with unprecedented detail, making it particularly suitable for testing theoretical wind-driving models against observational constraints. R Dor exhibits a semi-regular variable pattern with alternating pulsation modes having periods of 175 and 332 days (Bedding et al. 1998), an effective temperature of ∼3000–3100 K (Dumm & Schild 1998), and a moderate mass-loss rate of ∼10⁻⁷ M_⊙ yr⁻¹ (Van De Sande et al. 2018). High-resolution polarised-light (Khouri et al. 2016) and molecular-emission (Vlemmings et al. 2018; Decin et al. 2018; Nhung et al. 2021; Khouri et al. 2024) observations have revealed complex circumstellar features, including rotation and asymmetric structures, which provide crucial constraints for radiative-transfer modelling.

Earlier studies have focused primarily on fitting polarimetric observations without directly testing whether the observed amount and properties of the dust and gas are consistent with the wind-driving paradigm (e.g. Khouri et al. 2016; Ohnaka et al. 2016). Here, we analyse the dust in the immediate vicinity of R Dor in detail by systematically exploring a six-dimensional parameter space defined by dust grain size, density profile, and stellar radii, to derive constraints on the dust properties in the inner circumstellar environment. For the first time, we combine constraints from the dust model with results for the gas-density distribution (Khouri et al. 2024) to test whether such models can satisfy the physical requirements for wind-driving, combining force-balance calculations with elemental-abundance constraints. Our analysis of R Dor provides a crucial empirical test of the mass-loss mechanism in oxygen-rich AGB stars.

The structure of this paper is as follows. In Section 2, we describe the SPHERE/ZIMPOL observations and data reduction process. In Section 3, we outline the radiative-transfer modelling with RADMC-3D, including the definition of the dust density profile and the parameter space explored, and we present the methodology used to derive the optical properties of dust grains. Section 4 details the model grid and simulation strategy, describing the iterative approach used to optimise computational efficiency and ensure convergence. The results of our modelling, including constraints on dust grain size, composition, and polarimetric fits, are presented in Section 5. In Section 6, we examine elemental and radiative constraints through a detailed analysis of MgSiO₃ models using DHS scattering. Section 7 synthesises the wind-driving potential across all dust compositions and scattering theories. In Section 8, we compare our findings with previous studies and discuss their implications for AGB mass-loss mechanisms. Finally, Section 9 summarises our key results and suggests directions for future work.

2. Observations and data reduction

R Dor was observed using SPHERE/ZIMPOL on the Very Large Telescope on November 28 and 29, 2017 (ESO program 0100.D-0737). The observations were carried out with three filters in polarimetric mode, with cycles consisting of blocks of exposures for each of the four angles of the polarimeter. Observations of a point spread function (PSF) reference star (HD 25170) were obtained before those of R Dor, following the same procedure but with different exposure times. Table 1 summarises the observational setup and integration times for both the target and the PSF reference. The PSF maps are shown in Appendix A.

The data were reduced using the pipeline developed at ETH Zurich (Schmid et al. 2018), with the standard steps of bias subtraction, flat-field correction, and polarimetry modulation-demodulation efficiency correction. We adopted a value of 0.8 for the polarimetric efficiency following Schmid et al. (2018). The pipeline produces images of the Stokes parameters I (total intensity), Q, and U. We calculated the linearly polarised intensity using $I_{\mathrm{p}}=\sqrt{Q^2+U^2}$, the polarisation degree using p = I_p/I, and the angle of the polarised light using θ = arctan(U/Q)/2.

Figure 1 shows the polarisation degree and normalised total intensity maps of R Dor obtained with VLT/SPHERE-ZIMPOL at three wavelengths (0.65, 0.75, and 0.82 μm). These observations reveal the unpolarised stellar core surrounded by a bright arc corresponding to stellar light scattered by dust grains in the circumstellar environment of this AGB star. The scattered-light pattern provides direct evidence for the presence of dust grains in the immediate vicinity of R Dor.

Fig. 1. Polarisation degree maps (top row) and normalised total intensity maps (bottom row) of R Dor observed with VLT/SPHERE-ZIMPOL at three wavelengths: 0.65 μm (left), 0.75 μm (middle), and 0.82 μm (right). The red contours in the bottom right of each panel show the instrument PSF at the 10% and 50% intensity levels (see Appendix A for the complete PSF maps).

3. Radiative transfer modelling

3.1. Dust grains and optical properties

We used the publicly available Optool code (Dominik et al. 2021) to compute the optical properties (e.g. scattering phase functions, absorption efficiencies, and Mueller matrices) of dust grains commonly observed around AGB stars. Optool seamlessly interfaces with RADMC-3D, allowing for straightforward adoption of its outputs in radiative transfer simulations.

Our study considers three main dust compositions: magnesium silicate (MgSiO₃), aluminium oxide (Al₂O₃), and iron-bearing silicate (MgFeSiO₄). Although MgFeSiO₄ is prone to overheating in the inner regions of an oxygen-rich AGB star (Woitke 2006), we included it for comparison to assess the impact of its stronger near-IR absorption on the wind-driving efficiency. For each material, we computed optical properties using two scattering theories: (i) classical Mie theory for homogeneous, perfectly spherical grains and (ii) the DHS method, which approximates irregular grain shapes and can better reproduce observed polarisation signatures (Min et al. 2005).

In addition, we adopted a single dust size to explore the effects of grain size on scattering and absorption. Although real dust populations follow broader size distributions, constraining the observed polarisation and intensity profiles with a single grain size introduces fewer free parameters. This single-size approach thus offers a baseline for characterising the dust in R Dor. The radiative transfer code RADMC-3D (Dullemond et al. 2012) is a versatile tool widely used to model dust and gas in astrophysical environments. The tool uses Monte Carlo and ray-tracing techniques to simulate radiative processes such as scattering, absorption, and emission on a 3D spatial grid, making it ideal for investigating circumstellar environments like that of R Dor. Importantly, RADMC-3D’s ability to compute full Stokes parameters (I, Q, U) makes it particularly well-suited for modelling polarisation and for comparison with SPHERE/VLT observations.

Previous studies (e.g. Maercker et al. 2022; Wiegert et al. 2024) have used RADMC-3D for AGB star environments. Building on this foundation, we applied RADMC-3D to constrain dust grain size, density profiles, and compositions around R Dor. The combination of Optool and RADMC-3D ensured modelling of scattering phase functions and produced polarisation maps, enabling a detailed comparison with observations. The inputs to RADMC-3D were the dust density distribution and the grain optical properties. To consider direction-dependent scattering, the Mueller scattering matrices needed to be provided.

3.2. Wavelength-dependent stellar radii and dust density profile

The stellar radius, R_λ, clearly varies across the three wavelengths of observation (see Figure 1). We attribute this to changes in optical depth and molecular opacity. We consider the stellar radii corresponding to the three wavelengths of interest as free parameters: R_λ1 at λ₁ = 0.65 μm, R_λ2 at λ₂ = 0.75 μm, and R_λ3 at λ₃ = 0.82 μm. Accurately modelling R_λ is essential for capturing the radiative transfer effects that influence dust scattering and polarisation.

The dust density profile in the inner circumstellar envelope of R Dor was modelled using a power-law distribution

$$\begin{array}{c}\hfill \rho (r)={\rho }_0\times {\left(\frac{r}{r_0}\right)}^{p_{\rho }},\end{array}$$(1)

where ρ₀ represents the dust density at the reference radius r₀ = 180 R_⊙. Fixing r₀ allows the normalisation of the dust density to be decoupled from the wavelength-dependent stellar radii R_λ1, R_λ2, and R_λ3. This approach ensures consistency when comparing models across different wavelengths. The power-law exponent, p_ρ, governs the steepness of the density gradient.

3.3. Optimisation of RADMC-3D hyperparameters

Efficient and accurate modelling with RADMC-3D requires careful selection and optimisation of key hyperparameters to balance computational efficiency and the fidelity of the results. A critical aspect is the spatial grid resolution, which must be carefully chosen to adequately capture steep gradients in dust density near the star. Adaptive mesh refinement was not used because preliminary tests showed that a fixed resolution sufficiently captured the dust distribution while minimising computational costs.

For our simulations, we used a computational grid of 112 × 112 × 112 cells spanning a physical size of 23.78 × 23.78 × 23.78 au. The final pixel size in the images was determined by an intermediary step, as we employed a 2D grid for the ray-tracing stage of the radiative transfer modelling in RADMC-3D. In this step, we used a 112 × 112 pixel grid, which corresponds to a pixel size of 3.6 mas, matching the native SPHERE/ZIMPOL pixel scale. This configuration was optimised through a convergence study to ensure that the chosen resolution did not introduce significant numerical artefacts or biases in the results.

Another vital consideration is the number of photon packages used in the Monte Carlo simulations. The number of photon packages was adjusted depending on factors such as optical depth and grain size, ensuring convergence across a wide parameter space. The methodology for determining the appropriate photon count for each model is described in detail in Sect. 4.3.

4. Model grid and simulation strategy

4.1. Overall strategy and parameter space

The simulation workflow for modelling the circumstellar environment of R Dor follows a systematic approach designed to optimise computational efficiency while ensuring physical consistency with observations. Our strategy consists of four main steps: (1) exploring a comprehensive parameter space of dust properties and stellar radii; (2) implementing an iterative radiative transfer procedure with adaptive convergence checks; (3) applying a two-stage filtering process based on total intensity profiles to eliminate unphysical models; and (4) selecting best-fit models through combined polarimetric and intensity constraints. This approach provides a robust framework for constraining the dust grain properties, density profiles, and stellar radii while managing computational resources efficiently.

To systematically explore the parameter space governing dust properties and wind dynamics, we varied key parameters over a well-defined range. Each parameter was chosen to encompass plausible physical values while remaining computationally tractable. The grain size, dust density normalisation, and density profile exponent were sampled logarithmically or linearly, depending on the respective explored range. Similarly, the stellar radii at different wavelengths were varied across a range motivated by observational constraints. The specific parameter ranges and their sampling methods are summarised in Table 2.

Exploring such a large parameter space poses significant computational challenges, as the total number of models scales as nbins⁶. To address this, we carefully balanced grid resolution and computational efficiency, ensuring sufficient parameter coverage while minimising redundant calculations. Based on preliminary tests, we set the number of bins to 15, which was sufficient to accurately describe the different parameter ranges without being excessively large, thereby avoiding an impractical computation time. Additionally, models with extreme parameter combinations - unlikely to fit the observations - were identified and excluded early in the iterative simulation process through our two-stage filtering approach, significantly reducing unnecessary computations.

4.2. Observational constraints and diagnostic quantities

To ensure that our models reproduced both the polarimetric signatures and the overall intensity distribution of R Dor, we employed two complementary sets of diagnostic quantities. These metrics allowed us to quantitatively assess model performance and implement efficient filtering strategies throughout the simulation process.

4.2.1. Total intensity profile constraints

Scattering of radiation on dust grains has been identified as an important source of opacity affecting the sizes of evolved stars at visible wavelengths (Ireland et al. 2004). Additionally, molecular opacity effects, particularly from TiO absorption, contribute to wavelength-dependent variations in stellar radii, although these effects are generally less dominant than dust scattering in the circumstellar regions investigated here (Ireland et al. 2004). The amount of dust required to fit the observed amount of polarised light in R Dor is sufficient to significantly affect the stellar size at the observed wavelengths. Hence, it was necessary to verify whether the models that provided acceptable fits to the observed polarised-light images were consistent with the total-intensity images.

To compare the modelled and observed total intensity profiles, we analysed their normalised profiles by measuring the full widths at various fractions of the maximum intensity. These fractions represent specific intensity levels: 50% of the peak normalised intensity, defined as the full width at half maximum (FWHM), 10%, and 1%. Measuring full widths at different intensity levels provides complementary information: the FWHM traces the compact core of the stellar disc, while the full widths at 10% and 1% probe the extended scattered light distribution. This approach ensured that the focus remained on the overall shape of the intensity profiles.

To quantify the agreement between the modelled and observed intensity profiles, we defined the following diagnostic parameter:

$$\begin{array}{c}\hfill \mathrm{\Delta }{\mathrm{FW}}_j(\%)={\mathrm{FW}}_{\mathrm{mod},j}(\%)-{\mathrm{FW}}_{\mathrm{obs},j}(\%),\end{array}$$(2)

where FW_mod, j(%) and FW_obs, j(%) are the full widths at the specified intensity fraction for the modelled and observed profiles, respectively, at wavelength j. For example, ΔFW₃(50%) represents the difference in FWHM at the longest wavelength (λ₃ = 0.82 μm), while ΔFW₁(1%) corresponds to the difference at 1% of the maximum intensity for the shortest wavelength (λ₁ = 0.65 μm).

This diagnostic parameter served two critical functions in our analysis. First, it provided an early filtering criterion during the iterative simulation procedure: models with ΔFW_j(50%) > 200% (i.e. FWHM more than three times the observed value) at any wavelength were immediately excluded as these represented highly optically thick models incompatible with observations. Second, it served as a final selection criterion, where we required |ΔFW_j(%)| < 30% for all intensity levels and wavelengths. This tolerance threshold accounted for possible gas-emission contributions to the observed profiles, while ensuring that the modelled dust profiles reproduced the overall shape and spatial extent of the observed intensity distribution.

Gas emission and absorption can affect total intensity profiles, particularly in the outer circumstellar regions. Because the radiative transfer simulations focus solely on dust scattering and absorption, they do not account for thermal emission from the gas. This is a shortcoming of the current state-of-the-art modelling of these types of observations, but we expected a relatively small effect in the region where polarised light is observed. The 30% tolerance threshold was chosen to account for these effects while still maintaining meaningful constraints on the dust distribution.

4.2.2. Polarisation degree constraints

The polarisation degree profile as a function of radius served as the primary observable to constrain dust properties. We quantified the comparison between the modelled and observed polarisation profiles using chi-square analysis, with the detailed formulation presented in Section 5. We required the models to achieve χ_red, 3D² < 1.36 (corresponding to the 3σ confidence level) for them to qualify as acceptable fits to the polarimetric data. This constraint ensured that the selected models accurately reproduced the observed polarisation signatures across all three wavelengths.

4.3. Iterative simulation and model-selection procedure

The simulation procedure implemented an adaptive approach that balanced computational efficiency with convergence requirements. The workflow integrated early model filtering, iterative convergence checks, and systematic validation against observational constraints.

4.3.1. Radiative transfer setup and initialisation

The process began with the set of input parameters, including the dust grain size (a₀), density normalisation (ρ₀), density profile exponent (p_ρ), and stellar radii (R_λ1, R_λ2, R_λ3). Using these parameters, we generated the necessary input files for RADMC-3D, including the dust density profile and optical properties computed with Optool. We initialised the simulation with 20 000 photon packages, and the radiative transfer calculations were performed to generate maps of the Stokes parameters I, Q, and U.

4.3.2. Convergence assessment

To ensure a meaningful comparison with observations and assess Monte Carlo convergence, we implemented a systematic convergence check procedure. After each radiative transfer calculation, the modelled 2D intensity and polarisation-degree (PD) maps were first convolved with synthetic symmetrical PSFs. The synthetic PSF was constructed by computing the azimuthal average of the observed PSF at each radius and applying this averaged value uniformly across all azimuthal angles for that radius. This approach ensured that any azimuthal variations in the modelled maps at a fixed radius were purely due to Monte Carlo noise and not instrumental effects, as real PSF azimuthal variations could otherwise be misinterpreted as simulation noise.

The stability of the Monte Carlo simulations was evaluated by monitoring the standard deviation in the intensity and polarisation degree maps over successive iterations. We analysed the 2D total-intensity and PD maps convolved with the synthetic PSF by computing the standard deviation of values over azimuthal angles at each radius. If the standard deviation remained below 5% of the average value at all radii, the simulation was considered to be converged. If this threshold was not met, the number of photon packages was doubled, and the simulation was repeated. This iterative refinement continued for a maximum of 12 iterations, though post-analysis confirmed that convergence was always achieved within fewer steps.

4.3.3. Two-stage filtering process

To optimise computational resources and ensure physical consistency, we implemented a two-stage filtering approach. In the first stage, following the initial convergence assessment, we applied an immediate filter based on the total intensity profile width. If ΔFW_j(50%) > 200% (i.e. the modelled FWHM was more than three times the observed value) at any of the three wavelengths, the model was excluded from further computation. These models were typically highly optically thick, producing stellar radii significantly larger than observed, making it impossible to match the observed intensity profile. This pre-selection prevented the calculation of computationally intensive models that could not reproduce the observations.

In the second stage, for models that passed the early filtering and achieved convergence, a final validation was performed during the best-fit model selection process (detailed in Section 5.2). At this stage, we required |ΔFW_j(%)| < 30% for all intensity levels (50%, 10%, and 1%) and for all wavelengths. Only models satisfying both this intensity constraint and the polarisation chi-squared criterion were retained as viable solutions.

4.3.4. Final processing and output

Once convergence was achieved and the model passed the early filtering stage, the unconvolved intensity and PD maps were convolved with the real, asymmetric instrument PSFs, producing the final simulated maps that were directly comparable to observations. This final convolution step ensured that all instrumental effects were properly accounted for in the model-observation comparison.

The complete iterative simulation procedure provides a systematic and efficient framework for exploring the parameter space while managing computational time¹. By integrating early optical depth filtering, adaptive convergence checks, and multi-parameter validation, the procedure maximises efficiency and ensures that the final models remain physically realistic and consistent with observational constraints. Moreover, it provides a measure of the uncertainty in the models introduced by the fluctuations inherent to Monte Carlo simulations.

5. Results

This section presents a systematic exploration of the parameter space, the constraints derived from polarimetric fitting, and the critical wind-driving analysis.

5.1. Chi-squared analysis

The goodness of fit between the modelled and observed polarisation degree profiles is quantified using the chi-squared (χ²) statistic, which is calculated as

$$\begin{array}{c}\hfill {\chi }^2(a_0,{\rho }_0,p_{\rho },R_1,R_2,R_3)={\displaystyle \sum _{j=1}^3\sum _{i=1}^{50}}\frac{{(⟨P{D}_{\mathrm{obs},i,j}⟩-⟨P{D}_{\mathrm{mod},i,j}⟩)}^2}{{\sigma }_{\mathrm{obs},i,j}^2+{\sigma }_{\mathrm{mod},i,j}^2},\end{array}$$(3)

where ⟨PD_{obs, i, j}⟩ and ⟨PD_{mod, i, j}⟩ are the mean observed and modelled polarisation degrees at the i-th radial bin for the j-th wavelength, respectively. The uncertainties σ_{obs, i, j} and σ_{mod, i, j} represent the standard deviations of the observed and modelled polarisation degrees at each radial bin. The summation was performed over 50 radial bins, corresponding to a maximum radius of 180 mas, given the observational pixel size of 3.6 mas. Beyond this radius, the uncertainties in the observed polarisation degree become prohibitively large, with the standard deviation exceeding the polarisation degree itself.

The mean polarisation degree at a given radius, r_i, and wavelength, λ_j, is computed as

$$\begin{array}{c}\hfill ⟨P{D}_{i,j}⟩=\frac{1}{N_{\theta }}{\displaystyle \sum _{k=1}^{N_{\theta }}}P{D}_{i,j,k},\end{array}$$(4)

where PD_i, j, k represents the polarisation degree at radius r_i, wavelength λ_j, and angular position θ_k, and N_θ is the number of angular bins. This expression applies to both the observed (PD_obs) and modelled (PD_mod) data.

The standard deviation of the polarisation degree at a given radius and wavelength, accounting for azimuthal variations, is given by

$$\begin{array}{c}\hfill {\sigma }_{i,j}=\sqrt{\frac{1}{N_{\theta }-1}{\displaystyle \sum _{k=1}^{N_{\theta }}}{(P{D}_{i,j,k}-⟨P{D}_{i,j}⟩)}^2}.\end{array}$$(5)

This formulation applies separately to both the observed (σ_{obs, i, j}) and modelled (σ_{mod, i, j}) datasets. The reduced chi-square (χ_red²) is then calculated to account for the degrees of freedom in the fit:

$$\begin{array}{c}\hfill {\chi }_{\mathrm{red}}^2(a_0,{\rho }_0,p_{\rho },R_{{\lambda }_1},R_{{\lambda }_2},R_{{\lambda }_3})=\frac{{\chi }^2(a_0,{\rho }_0,p_{\rho },R_{{\lambda }_1},R_{{\lambda }_2},R_{{\lambda }_3})}{N-P},\end{array}$$(6)

where N = 150 is the total number of data points (50 radial bins for each of the three wavelengths) and P = 6 is the number of free parameters (a₀, ρ₀, p_ρ, R_λ1, R_λ2, and R_λ3). To remove the dependence of χ_red² on the stellar radii (R_λ1, R_λ2, R_λ3), we define a new metric, χ_red, 3D², which minimises χ_red² over the three radii and the same values of a₀, ρ₀, and p_ρ:

$$\begin{array}{c}\hfill {\chi }_{\mathrm{red},3\mathrm{D}}^2(a_0,{\rho }_0,p_{\rho })={\displaystyle \underset{R_{{\lambda }_1},R_{{\lambda }_2},R_{{\lambda }_3}}{min}}{\chi }_{\mathrm{red}}^2(a_0,{\rho }_0,p_{\rho },R_{{\lambda }_1},R_{{\lambda }_2},R_{{\lambda }_3}).\end{array}$$(7)

The models are classified based on the value of χ_red, 3D², with thresholds for 1σ, 2σ, and 3σ confidence levels set at 1.05, 1.21, and 1.36, respectively. These thresholds, derived using the chi-squared cumulative distribution function, represent acceptable deviations between the model and observations. In the following analysis, we only consider models that fall within the 3σ confidence level, meaning that only models with χ_red, 3D² < 1.36 are retained.

5.2. Selection of best-fit models

The selection of best-fit models was carried out in two steps: first, by applying thresholds based on the reduced chi-squared statistic (χ_red, 3D²) and second, by filtering models based on the normalised intensity profile differences (ΔFW_j(%)).

We began by selecting all models with χ_red, 3D² < 1.36, ensuring that they provide an acceptable fit to the observed polarisation data.

Among the models that pass this chi-squared threshold, we further refined the selection by comparing the modelled and observed normalised intensity profiles. This comparison was performed using the ΔFW_j(%) parameters, which quantify the differences in full widths at 50%, 10%, and 1% of the maximum intensity for each wavelength (j = 1, 2, 3). Specifically, we retained only models for which

$$\begin{array}{c}\hfill |\mathrm{\Delta }{\mathrm{FW}}_j(50\%)|,|\mathrm{\Delta }{\mathrm{FW}}_j(10\%)|,|\mathrm{\Delta }{\mathrm{FW}}_j(1\%)|<\mathrm{\Delta }{\mathrm{FW}}_{\mathrm{limit}},\end{array}$$(8)

where ΔFW_limit = 30%. We find that using a smaller limit barely affects the mean, minimum, and maximum values of the parameters we aim to constrain (a₀, ρ₀, p_ρ, R_λ1, R_λ2, and R_λ3). However, we adopted a 30% limit to ensure that a sufficient number of models satisfied the criterion for statistical analysis while still maintaining a close match between the modelled and observed intensity profiles.

By combining these two selection criteria, we ensured that the best-fit models provided a comprehensive match to both the polarisation and intensity data. This approach acknowledges that gas emission may contribute to the observed total intensity profiles while our dust-only radiative transfer models focus on reproducing the polarisation signatures. The use of ΔFW_j(%) as an additional filtering parameter strengthens the reliability of the model selection process by addressing discrepancies in intensity profile shapes across multiple wavelengths. In particular, models that provided an acceptable fit to the observed polarisation degree but presented broader total intensity profiles than observed were clearly incompatible with the observations.

Table 3 presents the constrained parameter ranges derived from our comprehensive modelling analysis for the three dust compositions using both DHS and Mie scattering theories. These results represent the ensemble of models that successfully reproduce the observed polarimetric signatures while maintaining physical consistency with the intensity profiles. Detailed analysis of individual best-fit models, including their specific parameter combinations and goodness-of-fit statistics, are provided in Appendix B. There is a remarkable consistency between Mie and DHS results for dust grain size, dust density, and stellar radii parameters, which strengthens our conclusion that the constraints we derive on these physical parameters are not significantly dependent on assumptions about grain morphology. As an illustrative example, Figure 2 demonstrates the excellent agreement between our best MgSiO₃ DHS model and observations across the three wavelengths. Similar comparisons for the other dust materials and scattering theories are presented in Appendix E.

Fig. 2. Radial profiles of the polarisation degree (top row) and normalised total intensity (bottom row) for R Dor at three wavelengths: 0.65 μm (left), 0.75 μm (middle), and 0.82 μm (right). The profiles are averaged over 360 degrees. The green curves with error bars represent the observational data. The red curves show one of the acceptable models obtained in Fig. B.1, while the blue shading indicates the density of overlapping acceptable models.

Having successfully identified dust configurations that reproduce the observed polarimetric signatures of R Dor, we now address the fundamental challenge that motivates this study by evaluating whether these observationally constrained dust populations can drive the stellar wind under physically realistic conditions. Our assessment incorporates two essential criteria that any viable wind-driving model must satisfy: (1) consistency with elemental gas depletion constraints derived from solar abundances and gas-phase observations, ensuring that the required dust masses do not exceed the available elemental budget; and (2) sufficient radiative pressure to overcome gravitational forces, quantified through the force balance parameter Γ.

6. Elemental and radiative constraints: MgSiO₃ analysis with DHS

Using solar photospheric elemental abundances, we derived (see Appendix C) absolute lower limits to the gas-to-dust mass ratio, considering full depletion of the limiting element (silicon for silicates and aluminium for alumina). Assuming complete depletion of these elements, we find minimum gas-to-dust mass ratios of 420 for MgSiO₃, 240 for MgFeSiO₄, and 9400 for Al₂O₃. Observations of SiO line emission from R Dor indicate that only 8–15% of the silicon is actually depleted (Van De Sande et al. 2018), raising the realistic limits for the gas-to-dust mass ratios to 2800 for MgSiO₃ and 1600 for MgFeSiO₄. For Al₂O₃, no observational constraints on aluminium depletion are available from Van De Sande et al. (2018), so we can only apply the theoretical complete depletion limit. We now examine whether our dust models that successfully fit the polarimetric observations can satisfy these physical constraints and drive stellar winds.

Figure 3 compares our observationally constrained dust density profiles considering DHS MgSiO₃ grains with gas density profiles from CO observations (Khouri et al. 2024). Following Khouri et al. (2024), we defined two regions based on the gas density profile transition: the inner region (R < 1.6 × R₈₈₇) and the outer region (R > 1.6 × R₈₈₇), where R₈₈₇ is the stellar radius at 887 μm as defined in Khouri et al. (2024). The analysis reveals stark differences between these regions. In the inner region, models with larger grain sizes (0.037, 0.056, and 0.085 μm) satisfy observational constraints, requiring silicon depletion levels consistent with or below the 8–15% (Van De Sande et al. 2018). However, the outer region presents a more restrictive scenario: none of our models align with the observed depletion constraints, as all would require silicon depletion levels exceeding the observationally determined values.

Fig. 3. Dust and gas properties around R Dor considering MgSiO3 dust using DHS scattering theory. Left panel: Density profiles showing the gas density (dashed black line with grey uncertainty region) from Khouri et al. (2024) and dust densities for different grain sizes (coloured lines) corresponding to models that fit the observational constraints. The gas density exhibits a transition at 1.6 × R887 (Khouri et al. 2024). Right panel: Gas-to-dust (GTD) ratios derived from the density profiles. The hatched area (GTD ≤ 420) represents unrealistic values where more silicon would be locked in dust than expected in the stellar atmosphere. The red horizontal band shows constraints from SiO-depletion measurements (Van De Sande et al. 2018). The background colour gradient indicates the percentage of available silicon locked in dust grains, from 1% (light blue) to 100% (dark blue).

Having established which dust models satisfy elemental abundance constraints, we now examine their wind-driving capability. The ratio Γ of radiative to gravitational acceleration (see Appendix D) serves as the critical test for dust-driven winds, with Γ > 1 indicating sufficient radiative pressure to overcome gravity.

Figure 4 illustrates Γ as a function of radius for MgSiO₃ dust grains calculated using the DHS approximation. The results are unambiguous: regardless of grain size or region, no models achieve Γ > 1, demonstrating that radiative pressure on dust cannot overcome gravitational forces. Specifically, models that satisfy the elemental abundance constraints in the inner region do not generate sufficient radiative pressure for wind acceleration.

Fig. 4. Radiative-to-gravitational force ratio (Γ) as a function of radius for MgSiO3 dust grains calculated using DHS scattering theory. Different coloured lines represent dust grain sizes. The green shaded area (Γ > 1) indicates regions where radiation pressure exceeds gravity. The hatched area (Γ < 1) represents regions where radiation pressure alone is insufficient to overcome gravity.

This analysis for MgSiO₃ using the DHS theory reveals a fundamental disconnect: while some dust models are consistent with elemental abundance constraints, none accounts for the expected radiation pressure to drive the observed stellar wind. Similar comprehensive analyses for all other dust compositions and scattering theories are synthesised in the following section to provide a complete assessment of wind-driving potential across all dust species.

7. Synthesis of wind-driving potential across dust compositions and scattering theories

We now synthesise results from all dust compositions and scattering theories by applying the two essential criteria established in the previous section. Figure 5 shows Γ as a function of the gas-to-dust mass ratio for all models in both inner (squares) and outer (circles) envelope regions across the three dust materials and two scattering theories.

Fig. 5. Radiative-to-gravitational force ratio (Γ) vs gas-to-dust mass ratio for dust models fitting the polarisation observations. Squares and circles denote inner and outer regions (R ≶ 1.6R887 from Khouri et al. 2024); colours indicate grain sizes. Five zones are shown, based on physical viability: (A) Γ > 1 but > 100% elemental depletion is impossible; (B) Γ < 1 and > 100% and depletion is impossible; (C) < 100% depletion but Γ < 1, representing insufficient radiation pressure; (D) Γ > 1 with high but possible depletion inconsistent with observations; and (E) Γ > 1 with observed depletion levels, the ideal wind-driving zone. The horizontal line marks Γ = 1.

We divided the parameter space into five distinct zones: (A) Γ > 1 but gas-to-dust mass ratios requiring excessive depletion (> 100% elemental depletion); (B) Γ < 1 and gas-to-dust mass ratios requiring excessive depletion (> 100% elemental depletion); (C) Gas-to-dust ratios requiring < 100% elemental depletion but Γ < 1; (D) Γ > 1 and gas-to-dust mass ratios requiring elemental depletion between observed levels and 100% (depletion potentially viable but inconsistent with observations); (E) Γ > 1 and gas-to-dust mass ratios consistent with observed depletion levels (ideal wind-driving zone).

Across all dust compositions and scattering theories, we find that no models fall within zone E, indicating that none can drive the outflow of R Dor under observationally realistic conditions.

Both MgSiO₃ and Al₂O₃ exhibit similar patterns regardless of the scattering theory. All models fall within zones B or C, with no instances reaching the wind-driving threshold of Γ > 1. In the inner region, both materials show acceptable gas-to-dust mass ratios (zone C) but insufficient radiation pressure, with Γ values typically ∼0.01. In the outer region, Al₂O₃ models fall into zone B due to excessive depletion requirements, while MgSiO₃ models remain in zone C. These results suggest a two-layer dust structure: Al₂O₃ can exist in the inner region but cannot explain the observed polarised light in the outer region, while MgSiO₃ can be present throughout both regions and accounts for the outer region polarisation.

In contrast, MgFeSiO₄ presents a different scenario. Some models with grain sizes of ∼0.056 μm appear to reach zone D in the inner region, suggesting potential wind-driving capability. However, these models face two critical limitations: (1) dust temperatures would exceed sublimation thresholds due to strong near-IR absorption (Woitke 2006), making their existence impossible in the acceleration zone; and (2) they require gas-to-dust mass ratios at the theoretical minimum for complete silicon depletion, contradicting the ∼10% depletion levels observed (Van De Sande et al. 2018).

Scattering theory comparison reveals minimal differences between Mie and DHS approaches. While Mie theory permits slightly larger maximum grain sizes (0.12 vs 0.09 μm), this difference is insufficient for any model to reach zone E. This demonstrates that grain model uncertainties do not affect our fundamental conclusion about wind-driving insufficiency.

The comprehensive analysis reveals a systematic failure across all dust compositions and scattering theories: while these models successfully reproduce observed polarimetric signatures, no model fits observational constraints while providing sufficient radiation pressure to drive stellar winds in R Dor.

8. Discussion

8.1. Two-layer dust envelope structure

Our zone analysis is consistent with a two-layer dust envelope: an inner gravitationally bound shell containing Al₂O₃ and an outer unbound envelope dominated by MgSiO₃. This stratification aligns with previous observations of W Hydrae (Khouri et al. 2015) and theoretical models (Höfner et al. 2016). Mid-infrared interferometry (Karovicova et al. 2013) supports this structure in oxygen-rich Mira variables, where alumina grains condense near the stellar surface (∼2 R_⋆) and serve as seed particles for silicates forming at larger distances ((4–5) R_⋆). Our detection of grains in the inner envelope and their predicted survival only at close distances (≤2 R_⋆) are consistent with theoretical models of dust-formation zones around AGB stars, where both Al₂O₃ and MgSiO₃ grains could reproduce the observed polarisation signatures. For the polarised-light observations to be reproduced, we find that silicate grains must exist relatively close to the star ≲2 R_⋆, using as a reference the stellar radius of 27.2 mas in the near-IR reported by Norris et al. (2012).

8.2. Dust grain size and density profile constraints

Understanding the optimal dust grain size for AGB wind-driving has evolved from theoretical predictions of ∼1 μm Fe-free silicate grains driving winds through scattering (Höfner 2008) to observational evidence of smaller grains ∼0.3 μm contributing to wind acceleration (Norris et al. 2012). Our R Dor results reveal even smaller grains of ∼0.1 μm, consistent with the findings of Khouri et al. (2016). This grain size variation appears common across oxygen-rich AGB stars: W Hydrae exhibits larger grains (0.4–0.5 μm) (Ohnaka et al. 2016), while Khouri et al. (2020) report a range from ≥0.1 μm (W Hya) to ≤0.1 μm (R Crt and SW Vir), suggesting grain growth processes depend on individual stellar characteristics.

Our zone analysis provides a crucial insight: despite confirming ∼0.1 μm grains in R Dor, these grains cannot independently drive winds under realistic elemental depletion constraints. The analysis considering MgFeSiO₄ demonstrates that elemental depletion requirements prevent even iron-bearing silicates from being potential wind drivers. These results for R Dor suggest that additional mechanisms must complement dust scattering to produce the observed mass loss.

While our modelling assumes single grain sizes, realistic dust populations follow size distributions. Our polarimetric observations are most sensitive to grains in the ∼0.01–0.3 μm range, where scattering efficiency peaks at visible wavelengths. Undetected portions of the size distribution are unlikely to significantly alter wind-driving conclusions: larger grains would produce detectable signatures if present in wind-driving quantities, while smaller grains contribute negligibly to radiation pressure. Most critically, elemental abundance constraints limit the total dust mass regardless of how it is distributed across grain sizes.

We also find a steep dust density profile (r^−3.4 to r^−4.1) that aligns with findings for R Dor (Khouri et al. 2016; Ohnaka et al. 2016). This steep profile is likely caused by the observed grains being located at the interface between gravitationally bound material and the outflow. The outward acceleration of the grains can also contribute to the observed steepness of the density profile, but this is probably a minor effect, given that the decrease is approximately 40 times steeper than for a constant-velocity outflow.

9. Summary and conclusion

Our multi-criteria analysis of R Dor indicates that MgSiO₃ and Al₂O₃ grains around this star cannot independently drive its stellar wind under the conditions observed during our study. By integrating polarimetric observations with physical constraints of elemental abundances and force balance, we find that even under the most generous assumptions of complete elemental depletion, the radiative force remains insufficient for wind acceleration.

Although Fe-bearing silicates (MgFeSiO₄) could theoretically provide adequate radiative pressure in the outer region, two significant challenges prevent them from driving winds: such grains would sublimate in the critical acceleration regions due to overheating (Woitke 2006), and they would require unrealistic gas-to-dust mass ratios that lock 80–90% of available silicon into dust grains – far exceeding the 10% depletion observed (Van De Sande et al. 2018).

While these results represent the first comprehensive test of dust-driven wind theory using observationally constrained models, we emphasise that they apply specifically to R Dor at the time of our observations. The wind-driving properties of dust could potentially vary with the stellar pulsation cycle, as shock waves propagating through the circumstellar environment may alter dust grain properties (size, composition, and spatial distribution) and local physical conditions (temperature, density, and velocity structure) that affect radiation pressure efficiency. Future multi-epoch observations spanning different pulsation phases would be valuable to assess whether dust-driven wind viability changes throughout the pulsation cycle. Our analysis reveals that dust contributes weakly to the momentum budget in this particular case and cannot independently sustain the observed outflow.

Our results suggest that for R Dor, dust likely plays a supporting rather than dominant role in mass loss, with pulsation-induced processes and other physical mechanisms contributing to the observed outflows. This finding motivates expanded studies of other oxygen-rich AGB stars with varying properties to determine whether similar constraints apply more broadly across this stellar population. Such comprehensive surveys will be essential to assess whether modifications to current mass-loss theories are needed and to develop improved models that capture the full complexity of AGB wind-driving physics.

Acknowledgments

We thank the referee for their very careful reading of our paper, which has significantly improved its quality. Financial support from the Knut and Alice Wallenberg foundation is gratefully acknowledged through grant no. KAW 2020.0081.

References

Appendix A: PSFs

Figure A.1 shows the point spread functions obtained from the reference star HD 25170 for the three wavelengths used in our analysis, displayed on a logarithmic scale.

Fig. A.1. Point Spread Function (PSF) maps for VLT/SPHERE-ZIMPOL observations of R Dor at three wavelengths: 0.65 μm (left), 0.75 μm (middle), and 0.82 μm (right), displayed on a logarithmic scale.

Appendix B: Best-fit models

B.1. Results for MgSiO₃ using DHS

Understanding how the six key parameters are constrained is crucial for interpreting the dust properties around R Dor. Table 3 summarises the constrained values for the dust grain size (a₀), dust density normalisation (ρ₀), density profile exponent (p_ρ), and the three wavelength-dependent stellar radii (R_λ1, R_λ2, R_λ3). The table includes results for three dust materials (MgSiO₃, Al₂O₃, and MgFeSiO₄), each modelled using both Mie theory and the DHS approximation.

Given the extensive parameter space explored across different materials and scattering theories, we first focus on MgSiO₃ using DHS to illustrate how we obtained these constraints. By analyzing the parameter space for this specific case, we show the process that led to the values presented in Table 3 before expanding the discussion to other dust compositions and scattering theories.

Figure B.1 illustrates the parameter space explored in this study for MgSiO₃ using DHS, highlighting the distributions and correlations between the six constrained parameters: log₁₀(a₀), log₁₀(ρ₀), p_ρ, and R_λ1, R_λ2, and R_λ3. The full parameter space consists of 15 bins per parameter, leading to a total of 15⁶ (1.14 × 10⁷) computed models. However, only 21 models satisfy the selection criteria described in the previous section, and it is these 21 models that are represented in Figure B.1.

The diagonal panels in Figure B.1 show the marginal distributions for each parameter, highlighting the minimum and maximum values for models that satisfy our constraints. It is crucial to emphasise that all 21 models represented in these distributions are equally valid solutions that successfully reproduce the observational constraints. For this reason, we focus exclusively on the minimum and maximum parameter values to define the uncertainty ranges, rather than reporting means or standard deviations which would be misleading in this context. Each model within these boundaries represents an equally good fit to the data, making the full range of values physically meaningful. The off-diagonal panels display scatterplots of pairwise parameter relationships, complemented by Spearman correlation coefficients to quantify the strength and direction of these correlations.

Fig. B.1. Parameter correlations and distributions for MgSiO3 for the selected models that satisfy the criteria of reduced chi-squared (χred2 < 1.36) at the 3σ confidence level with a normalised intensity tolerance of ΔFWlimit = ±30%. The diagonal panels show the marginal distributions of the six parameters: log10(a0) (dust grain size in microns), log10(ρ0) (dust density at r0 = 180 R⊙ in g/cm3), pρ (density power-law exponent), and stellar radii (Rλ1, Rλ2, Rλ3) in solar radii for the wavelengths 0.65, 0.75, and 0.82 microns, respectively. Vertical and horizontal dashed lines in the histograms represent the median parameter values, and shaded regions indicate parameter ranges. The off-diagonal panels display pairwise scatterplots of the parameters, with denser regions indicating higher model density.

B.1.1 Correlations between parameters

The scatter plots in the off-diagonal panels reveal several key correlations and interdependencies between parameters. A notable correlation is observed between ρ₀ and a₀, where an increase in a₀ corresponds to a decrease in ρ₀. This correlation reflects a degeneracy between these two parameters: a decrease in ρ₀ reduces the total dust cross section, which is compensated by an increase in the grain size, a₀. Physically, this suggests that larger grains require lower densities to match the observed scattering and absorption properties, as their interaction with radiation differs from that of smaller grains. Despite this degeneracy, a₀ never exceeds an upper limit of approximately 0.09 microns for DHS and 0.12 microns for Mie, regardless of the value of ρ₀, providing a strong constraint on the maximum dust grain size. The lower limit on a₀, however, depends on the upper limit of ρ₀, which we set at 10⁻¹⁴ g/cm³ based on constraints from the gas density (Khouri et al. 2024) and assuming a very optimistic gas-to-dust ratio of 100. This ratio, typical of the diffuse interstellar medium, is considerably higher than what would be expected in AGB star environments. Therefore, 10⁻¹⁴ g cm⁻³ represents a genuine physical upper bound for the dust density. While the lower limit is influenced by the upper limit of ρ₀, the key result here is the constraint on the maximum value of a₀.

A strong correlation is also observed between p_ρ and R_λ1, R_λ2, and R_λ3. As p_ρ increases (indicating a less steep density profile), the stellar radii decrease. This relationship likely stems from the interplay between the dust density distribution and the optical depth, which affects the observed stellar radii at different wavelengths. A less steep density profile leads to a more extended dust distribution, reducing the contribution of dust scattering near the star and effectively shifting the apparent stellar radii inward. Finally, there is a strong correlation between R_λ1, R_λ2, and R_λ3: an increase in one implies an increase in the others. This is expected because the stellar radii at different wavelengths are intrinsically linked through the temperature and opacity structure of the stellar atmosphere, and larger radii at shorter wavelengths indicate a globally extended stellar envelope.

B.1.2 Independent parameter distributions

While correlations reveal interdependencies between parameters, examining their independent distributions provides insights into their physically meaningful ranges. For MgSiO₃ with DHS scattering, we find the dust size, a₀, spans from 0.02 to 0.09 microns. The dust density normalisation, ρ₀, spans approximately three orders of magnitude, from 4 × 10⁻¹⁷ to 1 × 10⁻¹⁴ g/cm³. The power-law exponent governing the dust density profile, p_ρ, ranges from −4.4 to −2.8, indicating a steeper decline than the classical r⁻² profile expected for steady-state winds with constant outflow velocity.

For the stellar radii, we find consistent patterns across wavelengths. The radius at 0.65 μm (R_λ1) ranges from 288 to 414 R_⊙, while R_λ2 (0.75 μm) spans from 270 to 414 R_⊙, and R_λ3 (0.82 μm) varies between 216 and 378 R_⊙. The pattern R_λ1 > R_λ2 > R_λ3 reflects the wavelength-dependent scattering opacity in the stellar atmosphere, where shorter wavelengths correspond to higher optical depths and larger effective radii.

Additionally, we confirm that the absolute minima for all parameters fall within the boundaries of our exploration grid, ensuring that the parameter space was well-sampled and that the constraints derived from the models are trustable. This holds true for nearly all parameters, with one notable exception: the dust density parameter (ρ₀) consistently reaches its maximum value of 10⁻¹⁴ g cm⁻³, which corresponds to the upper boundary of our exploration grid. However, as discussed in our analysis of parameter correlations, this limitation primarily affects the minimum value of the grain size parameter (a₀) and does not significantly impact our main conclusions.

The polarisation degree panels display a characteristic pattern where polarisation initially increases with radius until reaching a maximum. This behavior occurs because close to the star, most emission originates from the unpolarised stellar surface. As the radius increases, the contribution of polarised scattered light from dust grains becomes more significant until reaching a maximum. Beyond this peak, polarisation degree decreases with increasing radius, corresponding to the decline in dust density, which is primarily responsible for the polarisation effect. Importantly, the observations consistently fall within the range defined by the overlapping models, indicating a reliable fit across all wavelengths. Similarly, the normalised total intensity profiles show excellent agreement between observations and models, exhibiting a smooth decline with radius.

B.2. Mie vs. DHS

For MgSiO₃, the results obtained using Mie and DHS scattering theories show notable similarities (Table 3). Both methods constrain the dust grain size a₀ to the same lower limit of 0.02 μm, with only a marginal difference in upper limits (0.09 μm for DHS versus 0.12 μm for Mie) - a distinction that has negligible impact on our conclusions. This consistency suggests that, for MgSiO₃, both scattering theories can reproduce the observed features with similar grain size populations despite their different treatments of grain shapes.

For the dust density normalisation, ρ₀, we observe some differences - 4 × 10⁻¹⁷ g cm⁻³ for the lower bound with DHS versus 10⁻¹⁷ g cm⁻³ for MIE - while both share the same upper bound of 10⁻¹⁴ g cm⁻³. The power-law exponent p_ρ spans nearly identical ranges for both scattering theories (-4.4 to -2.8), indicating that both approaches can reproduce the observed polarisation patterns with similar dust density distributions. The stellar radii parameters (R_λ1, R_λ2, and R_λ3) also show consistent ranges between the two methods, with only minor differences in the lower bound of R_λ1 (288 R_⊙ for DHS versus 306 R_⊙ for Mie). These similarities indicate that both scattering theories can effectively model the observed polarised light using physically comparable dust configurations.

For MgFeSiO₄ and Al₂O₃, we observe similar patterns of consistency between Mie and DHS results. The substantial overlap in parameter ranges across all dust species and both scattering theories suggests that our conclusions about dust properties in R Dor are reliable regardless of which scattering theory is applied.

The inclusion of both Mie and DHS calculations in our analysis serves an important methodological purpose. While Mie theory models perfectly spherical dust grains, the DHS approach approximates irregular grain shapes through a distribution of hollow spheres, offering a more realistic representation of actual cosmic dust particles, which are likely non-spherical in nature. By employing both methods, we test whether our conclusions about dust properties depend on assumptions about grain geometry.

B.3. Detailed parameter constraints for MgFeSiO₄ and Al₂O₃

Beyond MgSiO₃, we extended our analysis to MgFeSiO₄ and Al₂O₃ (Table 3). The results reveal notable distinctions in the behavior of these dust species.

The comparison between the modelled and observed polarisation degree and total intensity profiles is presented in Fig. E.1 for MgSiO₃ using Mie, in Fig. E.2 for Al₂O₃ using DHS and MgFeSiO₄ using Mie, and in Fig. E.3 for MgFeSiO₄ using DHS.

For MgFeSiO₄, the grain size parameter a₀ is constrained to 0.06 to 0.09 μm for DHS and 0.06 to 0.12 μm for Mie theory, as shown in Table 3. Notably, the minimum grain size for MgFeSiO₄ is larger than for MgSiO₃ (0.02 μm) across both scattering theories. The dust density normalisation ρ₀ for MgFeSiO₄ shows identical ranges for both scattering theories: 2 × 10⁻¹⁷ to 10⁻¹⁴ g/cm³. The power-law exponent p_ρ ranges from -4.8 to -3.0 for DHS and -4.4 to -3.0 for Mie, indicating a slightly steeper density gradient is possible when modelling with DHS for this material.

For Al₂O₃, the grain size parameter a₀ is constrained to 0.04 to 0.09 μm for DHS and 0.04 to 0.12 μm for Mie. The dust density normalisation ρ₀ for Al₂O₃ reaches the lowest values among the three materials: 6.3 × 10⁻¹⁸ to 10⁻¹⁴ g/cm³ for DHS and 10⁻¹⁷ to 10⁻¹⁴ g/cm³ for Mie. The power-law exponent p_ρ for Al₂O₃ spans -4.8 to -3.0 for both scattering theories. The stellar radii for Al₂O₃ show ranges of 306-414 R_⊙ for R_λ1 with DHS (306-432 R_⊙ with Mie), 216-396 R_⊙ for R_λ2 with DHS (198-414 R_⊙ with Mie), and 198-360 R_⊙ for R_λ3 with both theories.

For both Al₂O₃ and MgFeSiO₄, the stellar radii parameters (R_λ1, R_λ2, and R_λ3) show consistent wavelength-dependent patterns across all materials. MgFeSiO₄ models exhibit a somewhat broader range of stellar radii values, particularly at the shortest wavelength (0.65 μm), where DHS modelling yields radii as low as 234 R_⊙ compared to the 306 R_⊙ lower bound found with Mie. This difference may reflect how the scattering and absorption properties of iron-bearing silicates interact with stellar radiation differently than iron-free materials. As shown in Table 3, the stellar radii for MgFeSiO₄ range from 234-432 R_⊙ for R_λ1 with DHS (306-432 R_⊙ with Mie), 198-432 R_⊙ for R_λ2 with both theories, and 198-360 R_⊙ for R_λ3 with both theories, consistent with the values presented in our parameter exploration.

The overall consistency in parameter ranges between Mie and DHS results, despite some slight variations in upper and lower bounds, suggests that our conclusions about dust properties in R Dor are valid across different assumptions about grain geometry. A key finding across all three dust materials is the convergence on similar grain size constraints, with upper limits consistently around 0.1 μm regardless of the scattering theory employed or the dust composition. Equally significant is the power-law exponent p_ρ, which ranges from -4.8 to -2.8 across all materials - markedly steeper than the r⁻² profile typically expected in steady-state wind models.

Appendix C: Elemental-abundance constraints

When modelling dust-driven winds in AGB stars, it is essential to consider the elemental abundances available for dust formation. The total abundance of elements (both in gas phase and locked in dust grains) around R Dor is expected to be similar to solar abundances (Asplund et al. 2009). This places fundamental constraints on the maximum amount of dust that can form from the available gas, and consequently on the minimum gas-to-dust mass ratio that is physically realistic.

For the three dust species considered in this study, the limiting elements are silicon (Si) for silicate grains (MgSiO₃ and MgFeSiO₄) and aluminium (Al) for alumina (Al₂O₃). These elements have the lowest abundances relative to their stoichiometric requirements in the respective dust compounds, making them the bottleneck for dust formation.

The gas-to-dust mass mass ratio limit can be derived from the depletion fraction of the limiting element. For a given dust species, this relationship is expressed as:

$$\begin{array}{c}\hfill {\left(\frac{{\rho }_{\mathrm{g}}}{{\rho }_{\mathrm{d}}}\right)}_{\min }=\frac{1}{X_{\mathrm{dust}}}=\frac{1}{f_{\mathrm{X}}\times A_{\mathrm{X}}\times \frac{m_{\mathrm{dust}}}{m_{\mathrm{X}}}},\end{array}$$(C.1)

where X_dust is the dust mass fraction, f_X is the fraction of the limiting element X (Si or Al) depleted from the gas phase into dust grains, A_X is the abundance of element X relative to hydrogen by mass (based on solar abundances), and $\frac{m_{\mathrm{dust}}}{m_{\mathrm{X}}}$ is the ratio of the total dust mass to the mass of element X in the dust compound.

Using the solar photospheric abundances from Asplund et al. (2009), where log(N_Si/N_H)+12 = 7.51 and log(N_Al/N_H)+12 = 6.45, we can calculate the corresponding mass fractions relative to the total gas mass. This conversion involves several steps:

First, we convert the logarithmic abundance notation to linear number ratios:

$$\begin{array}{cc}\hfill \frac{N_{\mathrm{Si}}}{N_{\mathrm{H}}}& ={10}^{(7.51-12)}\approx 3.24\times {10}^{-5}\hfill \end{array}$$

$$\begin{array}{cc}\hfill \frac{N_{\mathrm{Al}}}{N_{\mathrm{H}}}& ={10}^{(6.45-12)}\approx 2.82\times {10}^{-6}\hfill \end{array}$$

Next, we convert these number ratios to mass ratios relative to hydrogen by multiplying by the ratio of atomic masses:

$$\begin{array}{cc}\hfill \frac{m_{\mathrm{Si}}}{m_{\mathrm{H}}}& =\frac{N_{\mathrm{Si}}}{N_{\mathrm{H}}}\times \frac{M_{\mathrm{Si}}}{M_{\mathrm{H}}}\approx 3.24\times {10}^{-5}\times \frac{28.09}{1.008}\approx 9.03\times {10}^{-4}\hfill \end{array}$$

$$\begin{array}{cc}\hfill \frac{m_{\mathrm{Al}}}{m_{\mathrm{H}}}& =\frac{N_{\mathrm{Al}}}{N_{\mathrm{H}}}\times \frac{M_{\mathrm{Al}}}{M_{\mathrm{H}}}\approx 2.82\times {10}^{-6}\times \frac{26.98}{1.008}\approx 7.55\times {10}^{-5}\hfill \end{array}$$

where M_Si = 28.09 amu, M_Al = 26.98 amu, and M_H = 1.008 amu are the atomic weights.

Finally, we account for the fact that hydrogen constitutes exactly 73.81% of the total gas mass in the solar composition, with the remainder consisting of helium (24.85%) and heavier elements (1.34%) (Asplund et al. 2009). The mass fractions relative to the total gas mass are therefore:

$$\begin{array}{cc}\hfill A_{\mathrm{Si}}& =\frac{m_{\mathrm{Si}}}{m_{\mathrm{total}}}=\frac{m_{\mathrm{Si}}}{m_{\mathrm{H}}}\times \frac{m_{\mathrm{H}}}{m_{\mathrm{total}}}\approx 9.03\times {10}^{-4}\times 0.7381\approx 6.7\times {10}^{-4}\hfill \end{array}$$

$$\begin{array}{cc}\hfill A_{\mathrm{Al}}& =\frac{m_{\mathrm{Al}}}{m_{\mathrm{total}}}=\frac{m_{\mathrm{Al}}}{m_{\mathrm{H}}}\times \frac{m_{\mathrm{H}}}{m_{\mathrm{total}}}\approx 7.55\times {10}^{-5}\times 0.7381\approx 5.6\times {10}^{-5}\hfill \end{array}$$

For MgSiO₃ (pyroxene), the mass fraction of silicon can be calculated from the atomic weights:

$$\begin{array}{cc}\hfill \frac{m_{\mathrm{Si}}}{m_{{\mathrm{MgSiO}}_3}}& =\frac{28.09}{24.31+28.09+3\times 16.00}\hfill \end{array}$$

$$\begin{array}{cc}& =\frac{28.09}{100.40}\approx 0.28\text{or}28\%\hfill \end{array}$$

For MgFeSiO₄ (olivine), the mass fraction of silicon is:

$$\begin{array}{cc}\hfill \frac{m_{\mathrm{Si}}}{m_{{\mathrm{MgFeSiO}}_4}}& =\frac{28.09}{24.31+55.85+28.09+4\times 16.00}\hfill \end{array}$$

$$\begin{array}{cc}& =\frac{28.09}{172.25}\approx 0.163\text{or}16.3\%\hfill \end{array}$$

For Al₂O₃ (alumina), the mass fraction of aluminum is:

$$\begin{array}{cc}\hfill \frac{m_{\mathrm{Al}}}{m_{{\mathrm{Al}}_2{\mathrm{O}}_3}}& =\frac{2\times 26.98}{2\times 26.98+3\times 16.00}\hfill \end{array}$$

$$\begin{array}{cc}& =\frac{53.96}{101.96}\approx 0.53\text{or}53\%\hfill \end{array}$$

Assuming complete depletion of the limiting element (f_X = 1), the minimum gas-to-dust mass ratios would be approximately 420 for MgSiO₃, 240 for MgFeSiO₄, and 9400 for Al₂O₃.

These constraints on the gas-to-dust mass ratio are incorporated into our assessment of wind-driving potential. Models that predict gas-to-dust mass ratios below these minimum values would require unrealistic levels of element depletion and are therefore excluded from consideration, regardless of their ability to fit the observed polarisation data.

Substituting the calculated values into our equation, we can express the minimum gas-to-dust mass ratio as a function of the depletion fraction for each dust species:

$$\begin{array}{cc}\hfill {\left(\frac{{\rho }_{\mathrm{g}}}{{\rho }_{\mathrm{d}}}\right)}_{\min ,{\mathrm{MgSiO}}_3}& \approx \frac{420}{f_{\mathrm{Si}}}\hfill \end{array}$$(C.2)

$$\begin{array}{cc}\hfill {\left(\frac{{\rho }_{\mathrm{g}}}{{\rho }_{\mathrm{d}}}\right)}_{\min ,{\mathrm{MgFeSiO}}_4}& \approx \frac{240}{f_{\mathrm{Si}}}\hfill \end{array}$$(C.3)

$$\begin{array}{cc}\hfill {\left(\frac{{\rho }_{\mathrm{g}}}{{\rho }_{\mathrm{d}}}\right)}_{\min ,{\mathrm{Al}}_2{\mathrm{O}}_3}& \approx \frac{9400}{f_{\mathrm{Al}}}\hfill \end{array}$$(C.4)

These expressions provide a direct way to estimate the minimum physically plausible gas-to-dust mass ratio for a given depletion fraction of the limiting element in each dust species.

Observations of gas-phase abundances in AGB stars suggest that depletion is typically incomplete. For R Dor specifically, we can estimate the silicon depletion fraction using the SiO abundance measurements from Van De Sande et al. (2018). Their observations indicate that the SiO abundance relative to H₂ in the circumstellar envelope of R Dor is between 5.5 × 10⁻⁵ and 6.0 × 10⁻⁵. Comparing this to the expected total silicon abundance (relative to H₂) of approximately 6.5 × 10⁻⁵ based on solar values (Asplund et al. 2009), we can calculate the fraction of silicon depleted into dust grains as follows:

$$\begin{array}{cc}\hfill f_{\mathrm{Si}}& =\frac{{[\mathrm{Si}]}_{\mathrm{total}}-{[\mathrm{Si}]}_{\mathrm{gas}}}{{[\mathrm{Si}]}_{\mathrm{total}}}\hfill \end{array}$$(C.5)

$$\begin{array}{cc}& =\frac{{[\mathrm{Si}]}_{\mathrm{total}}-{[\mathrm{SiO}]}_{\mathrm{gas}}}{{[\mathrm{Si}]}_{\mathrm{total}}}\hfill \end{array}$$(C.6)

$$\begin{array}{cc}& =\frac{6.5\times {10}^{-5}-(5.5\text{to}6.0)\times {10}^{-5}}{6.5\times {10}^{-5}}\hfill \end{array}$$(C.7)

$$\begin{array}{cc}& =\frac{(0.5\text{to}1.0)\times {10}^{-5}}{6.5\times {10}^{-5}}\hfill \end{array}$$(C.8)

$$\begin{array}{cc}& \approx 0.08\text{to}0.15\hfill \end{array}$$(C.9)

This calculation shows that only about 8-15% of the available silicon is depleted into dust grains. Substituting these depletion fractions (f_Si = 0.08 and 0.15) into Equations (C.2) and (C.3) yields minimum gas-to-dust mass ratios of approximately 2800-5300 for MgSiO₃ and 1600-3000 for MgFeSiO₄.

To assess whether our dust models are physically plausible, we compare the gas-to-dust mass ratios derived from our polarimetric observations with the minimum values calculated above. Our primary criterion is that the gas-to-dust mass ratio should exceed the minimum theoretical value derived from elemental abundances. This minimum threshold corresponds to 100% of the limiting element (silicon or aluminum) being locked into dust grains - a scenario that represents an absolute physical limit rather than a likely reality. We deliberately adopt this flexible criterion to give our dust models the maximum opportunity to satisfy physical constraints. Using this threshold allows us to identify which dust models are physically impossible (those requiring more than 100% of available elements) versus those that are at least theoretically possible, if perhaps implausible. When comparing our models to the more realistic depletion fractions derived from gas-phase observations (approximately 10% for silicon according to Van De Sande et al. (2018)), we can further distinguish between models that are merely possible and those that align with actual circumstellar conditions.

Appendix D: Force balance and Γ calculation

Dust-driven winds rely on the balance between stellar radiation pressure and gravitational forces acting on dust grains. To assess the wind-driving potential of the dust with the grain sizes, densities, and compositions derived in this study, we evaluate the force balance parameter Γ - the ratio of radiative pressure to gravitational force. Following Höfner & Olofsson (2018), we use an expression for Γ that accounts for wind dynamics by incorporating gas and dust velocities:

$$\begin{array}{c}\hfill \mathrm{\Gamma }=\frac{⟨{\kappa }_{\mathrm{pr}}⟩L_{\star }}{4\pi cG{M}_{\star }}\frac{v_{\mathrm{g}}}{r_{\mathrm{gd}}{v}_{\mathrm{d}}},\end{array}$$(D.1)

where L_⋆ is the stellar luminosity, c is the speed of light, G is the gravitational constant, M_⋆ is the stellar mass, v_g is the gas velocity, v_d is the dust velocity, r_gd represents the gas-to-dust mass ratio (defined as ρ_g/ρ_d), and ⟨κ_pr⟩ is the flux-mean radiative pressure opacity per unit mass of dust. This flux-mean opacity is defined as:

$$\begin{array}{c}\hfill ⟨{\kappa }_{\mathrm{pr}}⟩=\frac{\int {\kappa }_{\mathrm{pr}}(\lambda ,a)B_{\lambda }(T_{\star })d\lambda }{\int B_{\lambda }(T_{\star })d\lambda }.\end{array}$$(D.2)

Here, κ_pr(λ, a) is the radiative pressure cross-section at wavelength λ, and B_λ(T_⋆) is the Planck function at stellar temperature T_⋆. The radiative pressure cross-section accounts for both absorption and scattering contributions:

$$\begin{array}{c}\hfill {\kappa }_{\mathrm{pr}}(\lambda ,a)={\kappa }_{\mathrm{abs}}(\lambda ,a)+(1-g){\kappa }_{\mathrm{sca}}(\lambda ,a),\end{array}$$(D.3)

where κ_abs and κ_sca are the absorption and scattering cross-sections, respectively, and g is the scattering asymmetry parameter.

Appendix E: Additional figures

Fig. E.1. Same as Fig. 2 but for MgSiO3 and using the Mie theory (top) and for Al2O3 and using the Mie theory (bottom).

Fig. E.2. Same as Fig. 2 but for Al2O3 and using the DHS theory (top) and for MgFeSiO4 and using the Mie theory (bottom).

Fig. E.3. Same as Fig. 2 but for MgFeSiO4 and using the DHS theory.

All Tables

All Figures

	Fig. 1. Polarisation degree maps (top row) and normalised total intensity maps (bottom row) of R Dor observed with VLT/SPHERE-ZIMPOL at three wavelengths: 0.65 μm (left), 0.75 μm (middle), and 0.82 μm (right). The red contours in the bottom right of each panel show the instrument PSF at the 10% and 50% intensity levels (see Appendix A for the complete PSF maps).
In the text

Fig. 2. Radial profiles of the polarisation degree (top row) and normalised total intensity (bottom row) for R Dor at three wavelengths: 0.65 μm (left), 0.75 μm (middle), and 0.82 μm (right). The profiles are averaged over 360 degrees. The green curves with error bars represent the observational data. The red curves show one of the acceptable models obtained in Fig. B.1, while the blue shading indicates the density of overlapping acceptable models.

In the text

Fig. 3. Dust and gas properties around R Dor considering MgSiO3 dust using DHS scattering theory. Left panel: Density profiles showing the gas density (dashed black line with grey uncertainty region) from Khouri et al. (2024) and dust densities for different grain sizes (coloured lines) corresponding to models that fit the observational constraints. The gas density exhibits a transition at 1.6 × R887 (Khouri et al. 2024). Right panel: Gas-to-dust (GTD) ratios derived from the density profiles. The hatched area (GTD ≤ 420) represents unrealistic values where more silicon would be locked in dust than expected in the stellar atmosphere. The red horizontal band shows constraints from SiO-depletion measurements (Van De Sande et al. 2018). The background colour gradient indicates the percentage of available silicon locked in dust grains, from 1% (light blue) to 100% (dark blue).

In the text

	Fig. 4. Radiative-to-gravitational force ratio (Γ) as a function of radius for MgSiO3 dust grains calculated using DHS scattering theory. Different coloured lines represent dust grain sizes. The green shaded area (Γ > 1) indicates regions where radiation pressure exceeds gravity. The hatched area (Γ < 1) represents regions where radiation pressure alone is insufficient to overcome gravity.
In the text

Fig. 5. Radiative-to-gravitational force ratio (Γ) vs gas-to-dust mass ratio for dust models fitting the polarisation observations. Squares and circles denote inner and outer regions (R ≶ 1.6R887 from Khouri et al. 2024); colours indicate grain sizes. Five zones are shown, based on physical viability: (A) Γ > 1 but > 100% elemental depletion is impossible; (B) Γ < 1 and > 100% and depletion is impossible; (C) < 100% depletion but Γ < 1, representing insufficient radiation pressure; (D) Γ > 1 with high but possible depletion inconsistent with observations; and (E) Γ > 1 with observed depletion levels, the ideal wind-driving zone. The horizontal line marks Γ = 1.

In the text

	Fig. A.1. Point Spread Function (PSF) maps for VLT/SPHERE-ZIMPOL observations of R Dor at three wavelengths: 0.65 μm (left), 0.75 μm (middle), and 0.82 μm (right), displayed on a logarithmic scale.
In the text

Fig. B.1. Parameter correlations and distributions for MgSiO3 for the selected models that satisfy the criteria of reduced chi-squared (χred2 < 1.36) at the 3σ confidence level with a normalised intensity tolerance of ΔFWlimit = ±30%. The diagonal panels show the marginal distributions of the six parameters: log10(a0) (dust grain size in microns), log10(ρ0) (dust density at r0 = 180 R⊙ in g/cm3), pρ (density power-law exponent), and stellar radii (Rλ1, Rλ2, Rλ3) in solar radii for the wavelengths 0.65, 0.75, and 0.82 microns, respectively. Vertical and horizontal dashed lines in the histograms represent the median parameter values, and shaded regions indicate parameter ranges. The off-diagonal panels display pairwise scatterplots of the parameters, with denser regions indicating higher model density.

In the text

	Fig. E.1. Same as Fig. 2 but for MgSiO3 and using the Mie theory (top) and for Al2O3 and using the Mie theory (bottom).
In the text

	Fig. E.2. Same as Fig. 2 but for Al2O3 and using the DHS theory (top) and for MgFeSiO4 and using the Mie theory (bottom).
In the text

	Fig. E.3. Same as Fig. 2 but for MgFeSiO4 and using the DHS theory.
In the text