© The Authors 2023

1 Introduction

Asymptotic giant branch (AGB) stars are the late evolutionary stage of low- and intermediate-mass stars (0.8 M_⊙ ≲ M_⋆ ≲ 8 M_⊙). These stars exhibit mass-loss rates that range from 10⁻⁸ up to 10⁻⁵ M_⊙ yr⁻¹, with terminal wind speeds of 5–30 km s⁻¹ (Habing & Olofsson 2004; Ramstedt et al. 2008). To achieve such outflows, a mechanism is needed to overcome the stellar gravitational attraction. For AGB stars, the wind is believed to be a pulsation-enhanced dust-driven wind. A complex interplay between strong convection in the AGB atmosphere, and large-amplitude long-period pulsations, forms shock waves in the atmosphere (Freytag & Höfner 2008, 2023; Freytag et al. 2017). These shocks levitate the gas into sufficiently cool regions where it is able to condensate into dust. Dust particles can efficiently absorb stellar radiation, such that they are pushed outwards by the radiation force. When moving outwards, the dust collides with the surrounding gas and drags it along. This creates an efficient mechanism for mass loss around AGB stars (Lamers & Cassinelli 1999; Höfner & Olofsson 2018).

High-resolution observations of AGB stars have revealed complex structures in their outflows (Ramstedt et al. 2014; Kervella et al. 2016; Decin et al. 2020). One of the leading hypotheses to explain these morphologies is the presence of a binary companion, gravitationally shaping the outflow into complex morphologies. To investigate this hypothesis, 3D hydro-dynamic studies have been performed using both grid-based and smoothed particle hydrodynamics (SPH) codes (e.g. Theuns & Jorissen 1993; Mastrodemos & Morris 1999; Kim & Taam 2012; Saladino et al. 2018, 2019; Maes et al. 2021; Malfait et al. 2021; Aydi & Mohamed 2022; Lee et al. 2022). Most of these studies use the so-called free-wind approximation where the gravity of the mass-losing star is ignored, and none of the complexities of the wind-launching mechanism are included. To improve the modelling of the outflow simulations, recent attempts have been made to include more of the underlying physical mechanisms. Using a grid-based code, Chen et al. (2017, Chen et al. 2020) considered pulsations and an approximate form of dust opacity and radiative transfer, showing that these implementations alter the resulting outflow morphologies and could lead to the formation of circumbinary disks. In SPH codes, Aydi & Mohamed (2022) simulated pulsations to launch the wind following the 1D modelling of Bowen (1988), while Siess et al. (2022) implemented dust nucleation.

Previous studies (e.g. Theuns & Jorissen 1993; Mastrode-mos & Morris 1999; Kim & Taam 2012; Maes et al. 2021; Malfait et al. 2021) showed that the structures, emerging from the interaction of a companion with the AGB wind, depend on the relative magnitude of the wind and orbital velocities. Understanding structure formation with a mass-losing AGB star is therefore tightly linked with understanding the velocity of the wind. To model this wind velocity realistically, an accurate acceleration prescription is required, including a realistic description of the radiation force.

Radiative transfer is thus a key ingredient for stellar wind models, as it provides the driving force in dust-driven winds. However, except for ‘star-in-the-box’ type simulations (Freytag & Höfner 2023), performing full radiative transfer calculations on-the-fly in spatially extended simulations is not yet feasible due to the large computational cost. Therefore, one needs to resort to approximate radiative transfer descriptions to investigate the effect of radiation on the shaping of the wind on large scales, potentially perturbed by a companion.

In this paper, we present the implementation of a ray tracer in the SPH code PHANTOM, and use it to define four different radiative transfer approximations, each of which yielding a different prescription for the radiation force. Each approximation increases the complexity, from no radiative transfer (free-wind approximation), to geometrical dilution (geometrical approximation), to accounting for radiation not exclusively coming from the AGB star (Lucy approximation), and attenuation of the stellar radiation (attenuation approximation). We analyse the effects of the approximations on the velocity profiles and morphological structures, considering two mass-loss rate regimes, both in the case of a single AGB star, as well as in an AGB binary system. Finally, we compare the Lucy and attenuation approximations to the results of the 3D radiative transfer code MAGRITTE.

The outline of this paper is as follows. In Sect. 2, we present the four different radiative transfer prescriptions investigated in this study, as well as their numerical implementation. In Sect. 3, the effects of these prescriptions on the velocity profile of a single AGB star are investigated. In Sect. 4, the changes on the morphological structures of a binary system with a primary mass-losing AGB star are analysed. A comparison between the approximations and the full 3D radiative transfer code MAGRITTE is given in Sect. 5, and the main results are summarized in Sect. 6.

2 Model and setup

2.1 Smoothed particle hydrodynamics

The smoothed particle hydrodynamics (SPH) code PHANTOM (Price & Federrath 2010; Lodato & Price 2010; Price et al. 2018) is used for our 3D hydrodynamic simulations. In the framework of SPH (Gingold & Monaghan 1977; Lucy 1977), the density distribution of a particle, the equations of motion, and energy conservation read (Price et al. 2018): $${\rho }_i={\displaystyle \sum _j{m}_{j}W\left(\left|r_i-r_j\right|,h_i\right),}$$(1) $$\frac{\text{d}{\upsilon }_i}{\text{d}t}=-{\displaystyle \sum _j{m}_j\left[\frac{P_i+q_i}{{\rho }_i^2\text{Ω}}{\nabla }_i{W}_{ij}\left(h_i\right)+\frac{P_j+q_j}{{\rho }_j^2{\text{Ω}}_j}{\nabla }_i{W}_{ij}\left(h_j\right)\right]+a_{\text{ext},i},}$$(2) $$\frac{\text{d}{u}_i}{\text{d}t}=\frac{P_i}{{\rho }_a^2{\text{Ω}}_i}{\displaystyle \sum _j{m}_j}\left({\upsilon }_i-{\upsilon }_j\right)\cdot {\nabla }_a{W}_{ab}\left(h_a\right)+{\text{Λ}}_{\text{shock}}-\frac{{\text{Λ}}_{\text{cool}}}{\rho },$$(3)

where i denotes the SPH particle, and j its neighbouring particles. Particles are defined by their (constant) mass m, position r, velocity υ, and specific internal energy u. The calculation of the local density ρ and pressure P requires knowledge of the neighbours, which are found within the support of the smoothing kernel W (particles within R_kernh, where R_kern is the (dimen-sionless) kernel radius and h the smoothing length). We use the M₄ cubic spline kernel, which vanishes outside a radius r = 2h. This leads on average to a number of neighbours equal to N_neigh = 57.9 (Price et al. 2018). In Eq. (2), α_ext,i represents the external accelerations applied to the particle i. In this setup, the forces are the gravitational attraction of the star(s) and of the additional radiation force on the particle. In a binary system, it is given by $$a_{\text{ext,}i}=-\frac{G{M}_{\text{AGB}}}{r_{i,1}^2}\left(1-{\text{Γ}}_i\right){\widehat{r}}_{i,1}-\frac{G{M}_{\text{comp}}}{r_{i,2}^2}{\widehat{r}}_{i,2},$$(4)

where r_i,1 and r_i,2 are the distances from the position of the i’th particle to the AGB star and the companion, respectively. Γ_i is the Eddington factor, which is given by $${\text{Γ}}_i=\frac{\left({\kappa }_d+{\kappa }_g\right)F_i{r}_{i,1}^2}{G{M}_{\text{AGB}}c},$$(5)

in which k_d and k_ɡ are the opacity of the dust and gas, respectively, and F_i is the flux coming from the AGB star, reaching the i’th particle. The widely adopted free-wind approximation implies setting Γ_i equal to one, not taking into account the potentially very complex behaviour of Γ (see Sect. 2.2.1 for a more detailed description).

The calculation of the dust opacity is complex. Several prescriptions are available in PHANTOM, including the complex nucleation theory (Siess et al. 2022), which is based on the theory of moments developed by Gail & Sedlmayr (2013). However, for the purpose of this paper, to reduce the computational cost and the complexity of the interpretation of the results, we use the simplified analytic dust opacity given by Bowen (1988) $${\kappa }_d\left(T_{\text{eq}}\right)=\frac{{\kappa }_{\max }}{1+\exp \left[\left(T_{\text{eq}}-T_{\text{cond}}\right)/\delta \right]},$$(6)

where T_eq is the dust temperature, T_cond = 1200 K the dust condensation temperature, K_max = 6 cm² g⁻¹ the maximal dust opacity, and δ = 60 K the temperature range over which dust condensation occurs. These parameters can be changed to model different types of dust, where this setup describes carbon-rich dust (Bowen 1988). For the gas opacity, a constant value of K_ɡ = 2 × 10⁻⁴ cm² g⁻¹ is chosen (Bowen 1988).

In Eq. (3), Λ_shock represents the energy dissipation rate required to give the correct entropy increase in shocks. It consists of viscous shock heating, artificial thermal conductivity, and artificial resistivity if magnetic fields are included (for details, see Price et al. 2018). The cooling rate is represented by Λ_cool. Two processes have been included in our simulations. They correspond to the Bowen (1988) prescription which is expressed as $${\text{Λ}}_{\text{Bowen}}=\frac{3R}{2\mu }\frac{T_g-T_{\text{eq}}}{C\prime },$$(7)

where R is the gas constant, µ the mean molecular weight, and C′ the parametric cooling rate, taken to be 3 × 10⁻⁵ g s cm⁻³ (Bowen 1988). This expression, which is proportional to the temperature difference between the gas and the dust, was constructed to mimic diffusion between the two species. Furthermore, cooling by neutral hydrogen, Λ_H, is also considered, following the formula given by Spitzer (1978). The total cooling rate is thus given by Λ_cool = Λ_Bowen + Λ_H.

2.2 Radiative transfer approximations

The transport of energy via radiation can be considered along specific directions, often referred to as rays. The frequency-dependent (indicated with the subscript ν) radiative transfer equation along a ray reads $$\frac{\text{d}{I}_v}{\text{d}s}={\eta }_v-{\alpha }_v{I}_v.$$(8)

Here, I_ν represents the intensity (i.e. the observable quantity), ds is a path element along the ray, η_ν is the emission coefficient, which is a measure for the radiative energy that is gained along to the ray, and α_ν is the absorption coefficient, which quantifies the radiative energy lost along the ray (α_ν = K_νρ). By defining the frequency-dependent optical depth, τ_ν, as $$\text{d}{\tau }_v={\alpha }_v\text{d}s,$$(9)

the transfer equation can be re-written as $$\frac{\text{d}{I}_v}{\text{d}{\tau }_v}=S_v-I_v,$$(10)

where S_ν = η_ν/α_ν is referred to as the source function. The transfer equation can formally be solved, yielding $$I_v\left(s\right)=I_v\left(0\right)e^{-{\tau }_v\left(s\right)}+{\displaystyle {\int }_0^{{\tau }_v\left(s\right)}{S}_v{e}^{-\left({\tau }_v\left(s\right)-\tau \right)}\text{d}\tau }.$$(11)

For a homogeneous medium, the source function does not change throughout the medium, and can be moved outside of the integral, yielding the solution $$I_v=I_v\left(0\right)e^{-{\tau }_v}+S_v\left(1-e^{-{\tau }_v}\right).$$(12)

From the intensity, we define the mean intensity J_ν and flux F_ν $$J_v=\frac{1}{4\pi }{\displaystyle \oint I_v\left(\theta ,\varphi \right)\text{dΩ,}}$$(13) $$F_v={\displaystyle \oint I_v\left(\theta ,\varphi \right)\text{cos\hspace{0.17em}}\theta \text{\hspace{0.17em}dΩ,}}$$(14)

where dΩ = sin θdθdф is the differential solid angle, and θ the angle between the direction of the intensity and a normal vector of the surface through which we are considering the flux. Under the condition of local thermodynamic equilibrium (LTE), the source function equals the Planck function, S_ν = B_ν(T), and $${\displaystyle {\int }_0^{\infty }{\kappa }_v{B}_v{\text{d}}_v={\displaystyle {\int }_0^{\infty }{\kappa }_v{J}_v{\text{d}}_v.}}$$(15)

In general, knowledge of the mean intensity, J_ν, requires to solve the radiative transfer equations throughout the entire medium. However, some approximations can allow us to avoid this costly computation. In our study, we restrict ourselves to the grey case (i.e. ignoring any frequency dependence). Assuming radiative equilibrium, the frequency-integrated mean intensity and flux are given by $$J={\displaystyle {\int }_0^{\infty }{J}_v\text{d}v={\displaystyle {\int }_0^{\infty }{B}_v\text{d}v=\frac{{\sigma }_{\text{sb}}}{\pi }}{T}^4,}$$(16) $$F={\displaystyle {\int }_0^{\infty }{F}_v\text{d}v,}$$(17)

where σ_sb is the Stefan-Boltzmann constant. The computation of the variables F and T_eq, using a full radiative transfer description, is too computationally demanding. In order to alleviate this problem, we investigate four different approximate descriptions, and evaluate their applicability for 3D stellar wind models.

2.2.1 Free-wind approximation

The free-wind approximation is the most drastic approximation, since no explicit treatment of the radiation transport is included, and Γ (Eq. (5)) is simply set equal to one, implying the SPH particles do not feel the gravitational pull of the mass-losing star, as it is artificially balanced by the radiation force. Despite the crudeness of this approximation, it provides a simple way to launch a wind without requiring a full treatment of pulsations and dust formation. This approximation was introduced by Theuns & Jorissen (1993) and has been widely adopted since, because of its simplicity (e.g. Mastrodemos & Morris 1999; Kim & Taam 2012; Liu et al. 2017; Saladino et al. 2019; Maes et al. 2021; Malfait et al. 2021; Lee et al. 2022).

2.2.2 Geometrical approximation

Assuming a spherical star of radius R_⋆, an SPH particle located at a distance r from the source ‘sees’ the star over an opening angle θ_M, given by sinθ_M = R_⋆/r. Furthermore, assuming that the star is isotropically emitting a constant intensity I_ν, such that I_ν = J_ν, the flux received at the particle’s location is $$\begin{array}{l}{F}_v\left(r\right)=I_v{\displaystyle {\int }_0^{2\pi }\text{d}\varphi {\displaystyle {\int }_0^{{\theta }_{\text{M}}}\cos \theta \sin \theta \text{d}\theta }}\hfill \\ =\pi I_v{\sin }^2{\theta }_{\text{M}}=F_v\left(R_{★}\right)\frac{R_{★}^2}{r^2}.\hfill \end{array}$$(18)

In the grey approximation, neglecting all frequency dependencies, the expression for Γ (Eq. (5)) takes the standard form $$\text{Γ=}\frac{\left({\kappa }_d+{\kappa }_g\right)L_{\text{AGB}}}{4\pi cG{M}_{\text{AGB}}},$$(19)

where $L_{\text{AGB}}=4\pi R_{*}^{2}F\left(R_{*}\right)$ is the AGB luminosity. In this expression, the Eddington factor depends on the dust opacity, which requires the dust temperature. From Eq. (13), the mean intensity reads $$\begin{array}{l}{J}_v\left(r\right)=\frac{I_v}{4\pi }{\displaystyle {\int }_0^{2\pi }\text{d}\varphi {\displaystyle {\int }_0^{{\theta }_{\text{M}}}\sin \theta \text{d}\theta \text{=}\frac{1}{2}{I}_v\left(1-\cos {\theta }_{\text{M}}\right)}}\hfill \\ =\frac{1}{2}\left[1-\sqrt{1-{\left(\frac{R_{★}}{r}\right)}^2}\right]I_v=W\left(r\right)I_v,\hfill \end{array}$$(20)

where W(r) is often referred to as the geometrical dilution factor. Assuming radiative equilibrium, the local mean intensity can be related to the radiative equilibrium temperature by $$T_{\text{eq}}^4\left(r\right)=\frac{\pi }{{\sigma }_{\text{sb}}}J\left(r\right)=\frac{1}{2}\left(1-\sqrt{1-{\left(\frac{R_{★}}{r}\right)}^2}\right)T_{{}_{★}}^4.$$(21)

This temperature is associated to the dust, because it absorbs most of the radiation due to its high intrinsic opacity. In this simple approximation, Γ and T_eq can both be evaluated locally, making it easy to implement in a hydrodynamics code. However, we should emphasize that these expressions heavily rely on the assumptions of spherical symmetry, radiative equilibrium, and that the local mean intensity is dominated by the intensity of the AGB star. This approach was used, for instance, by Aydi & Mohamed (2022).

2.2.3 Lucy approximation

Assuming spherical symmetry, local thermodynamic equilibrium (LTE; which implicitly assumes that collisions dominate over radiative processes, see e.g. Gail & Sedlmayr 2013, Sects. 8.1 and 8.2), and an optically thin extended envelope, an improved prescription for J can be obtained, which lifts the assumption that J is dominated by the intensity of the AGB star. Following Lucy (1971, 1976; also described in Appendix A.1.2 of Gail & Sedlmayr 2013), the frequency-integrated mean intensity J, as a function of distance r from the central star, reads (Eq. (12) from Lucy 1971) $$J\left(r\right)=\left[\frac{1}{2}\left(1-\sqrt{1-{\left(\frac{R_{★}}{r}\right)}^2}\right)+\frac{3}{4}{\tau }_{\text{L}}\right]J\left(R_{★}\right),$$(22)

where the Lucy optical depth τ_L is given by $${\tau }_{\text{L}}={\displaystyle {\int }_r^{\infty }\left({\kappa }_d+{\kappa }_g\right)\rho {\left(\frac{R_{★}}{r\prime }\right)}^2\text{d}r\prime }.$$(23)

The Lucy optical depth can physically be interpreted as quantifying the radiation that is absorbed by the surrounding envelope, outside a given distance r from the AGB star. This radiation is re-emitted isotropically (because of radiative equilibrium) and provides a positive feedback contribution to the mean intensity at the considered location r. From Eq. (16), we immediately get the equilibrium temperature (Eq. (3) of Lucy 1976): $$T_{\text{eq}}^4\left(r\right)=\frac{\pi }{{\sigma }_{\text{sb}}}J\left(r\right)=\left[\frac{1}{2}\left(1-\sqrt{1-{\left(\frac{R_{★}}{r}\right)}^2}\right)+\frac{3}{4}{\tau }_{\text{L}}\right]T_{{}_{★}}^4.$$(24)

If the spherical symmetry is broken, it is still possible to estimate the dust temperature at a given position by replacing the angle-independent radial coordinate (r) by (r, θ, φ) in the calculation of τ_L. In this generalized form, Eq. (24) writes $$T_{\text{eq}}^4\left(r,\theta ,\varphi \right)=\left[\frac{1}{2}\left(1-\sqrt{1-{\left(\frac{R_{★}}{r}\right)}^2}\right)+\frac{3}{4}{\tau }_{\text{L}}\left(r,\theta ,\varphi \right)\right]T_{{}_{★}}^4,$$(25)

where θ and φ are the azimuthal and polar angles, respectively, indicating the direction of the ray, originating from the AGB star. This approach can better account for local dusty regions in the simulation. The Lucy approximation is often used in wind simulations to estimate the dust temperature (e.g. Bowen 1988; Saladino et al. 2018, 2019; Lee et al. 2022).

2.2.4 Attenuation approximation

If the medium between the AGB star and an SPH particle is opaque and only absorbs the stellar radiation without radiating itself (i.e. S_ν = 0), then Eq. (12) simplifies to $$I_v=I_v\left(0\right)e^{-{\tau }_v},$$(26)

which does not depend on the symmetry of the problem. In the grey approximation, the flux at the particle’s location then becomes $$F\left(r,\theta ,\varphi \right)=F\left(R_{★}\right)\frac{R_{★}^2}{r^2}{e}^{-\tau \left(r,\theta ,\varphi \right)}=\frac{L_{\text{AGB}}}{4\pi r^2}{e}^{-\tau \left(r,\theta ,\varphi \right)},$$(27)

where the optical depth τ is given by $$\tau \left(r,\theta ,\varphi \right)={\displaystyle {\int }_{R_{★}}^r\left({\kappa }_d+{\kappa }_g\right)\rho \text{d}r\prime }.$$(28)

Likewise, the expression for the dust temperature now reads $$T_{\text{eq}}^4\left(r,\theta ,\varphi \right)=\frac{1}{2}\left(1-\sqrt{1-{\left(\frac{R_{★}}{r}\right)}^2}\right)T_{{}_{★}}^4{e}^{-\tau \left(r,\theta ,\varphi \right)}.$$(29)

Here, the Lucy optical depth (Sect. 2.2.3) has disappeared because there is no emission, and, hence, no external source of radiation can heat up the medium. This prescription was used, for instance, by Chen et al. (2017, 2020) in their study of the morphology of AGB outflows, perturbed by a stellar companion.

2.3 Ray-tracing implementation in Phantom

The Lucy and attenuation approximations, described above, require the calculation of the optical depth τ or τ_L. These are non-local quantities to be evaluated along the line of sight, connecting an SPH particle and the AGB star. For the computation of these quantities, we implemented the ray tracer described in this section¹. More details, for instance, about the trade-offs made during development, can be found in Appendix A.

2.3.1 Algorithm

Since SPH particles are not located on a mesh, we opt for the meshless ray tracer that is implemented in MAGRITTE (De Ceuster et al. 2020a,b). MAGRITTE is a 3D radiative transfer code, developed to handle both mesh-based and meshless data, and can therefore work with SPH.

The algorithm is visually explained in Fig. 1. Starting from a position P₁, a ray is traced in the R-direction. The algorithm searches for the nearest neighbours, using the information provided| by the kd-tree of the SPH particles available in PHANTOM. From these neighbours, it selects the point closest to the ray, provided that the point lies in the direction of the ray and is within a distance R_kernh, where h is the smoothing length and R_kern = 2. This search is repeated until a given location, or the edge of the numerical domain, is reached.

To calculate optical depths, a segment will be defined as the distance between two subsequent points projected on the ray, for instance, between K₁ and K₂ (the projections of P₁ and P₂, respectively) in Fig. 1. The optical depth can then be approximated by a sum over all segments created. using the ray-tracer. As such, the optical depth can be discretized as $$\tau ={\displaystyle {\int }_{R_{★}}^r\text{d}\tau \approx {\displaystyle \sum _i\text{Δ}{\tau }_i={\displaystyle \sum _i\left(\frac{{\kappa }_i{\rho }_i+{\kappa }_i{}_{+1}{\rho }_i{}_{+1}}{2}\right)\text{Δ}{s}_i},}}$$(30)

where the index, i, ranges over all points encountered along the ray, form the stellar surface at R_⋆ to the radial distance from the AGB star, r, Δs_i = dist(K_i+1, K_i) is the length of the segment, K_i, and ρ_i are the opacity and density, respectively, at point K_i. Δs_i can easily be determined from geometrical considerations, while K_i ρ_i requires a little more attention. It is evaluated using the SPH smoothing kernel, considering only the particles within two smoothing lengths, yielding $${\kappa }_i{\rho }_i={\displaystyle \sum _{j:\left|r_i-r_j\right|<2h_i}{m}_j}{\kappa }_{j}W\left(r_i-r_j,h_j\right).$$(31)

The opacity K_i ρ_i is computed at each point, K_i, along the ray and can then be used in Eq. (30) to give the integrated optical depth.

Fig. 1 Visual representation of the ray-tracing algorithm. Starting from the upper right corner at point P1 and following the direction of the ray, the subsequent point P2 is determined within the sphere of influence of radius Rkernh. The procedure is then repeated until the boundary of the domain is reached (see text for explanations).

2.3.2 3D optical depth interpolation

To calculate the optical depth throughout the numerical domain, the most physically correct approach is to trace a ray from each SPH particle to the AGB star. This approach is similar to that of Kessel-Deynet & Burkert (2000). However, it is very computationally intensive and inefficient, because the optical depths obtained for the particles close to the AGB star are re-calculated whenever considering particles further out. one way around this, is to trace a predefined number of rays originating from the star, and interpolate the optical depth from the values of τ on these rays. The idea of tracing only a predefined number of rays, is common to most 3D ray-tracing radiative transfer algorithms (see e.g. Altay & Theuns 2013; De Ceuster et al. 2020a). Starting from the surface of the star (at radius R_⋆), rays are traced outwards, and the optical depth increments, Δτ_i, are calculated, accumulated, and stored for each segment and for each ray. Now, to obtain the optical depth at each SPH particle, the following interpolation scheme is used. First, for each SPH particle, located at a distance r from the star, the four rays passing closest to the particle are identified. Then, the optical depth along each of these four rays is linearly interpolated at the distance r. Finally, the optical depth values of the four rays (all evaluated at r) are used to estimate τ at the particle’s location. The accuracy of this approach highly depends on the number of rays, as well as on their spatial distribution (see Appendix A.3 for more details).

To obtain a uniform distribution of rays in 3D, we use the HEALPIX package (Górski et al. 2005), specifically designed for this problem. HEALPIX divides the 2-sphere into isolat-erally distributed pixels of equal area. We use the unit vectors pointing to the centre of each pixel as the directions for the rays. In the nested scheme, HEALPIX can refine the spatial division by increasing the number of pixels, such that each pixel is split into four. Each time this happens, the order of the scheme increases by one. So, at order zero, the 2-sphere is divided into 12 pixels, and for a general order o, the number of pixels is n_rays = 12 × 4°. HEALPIX also provides functions to obtain the pixel associated with every point on the 2-sphere. This means that when the pixel unit vectors are used as directions along which to trace our rays, every point in the simulation can automatically be associated with a ray. In our simulations we use HEALPIX-order 5, which corresponds to n_rays = 12288.

3 Single-star models

To investigate the impact of each of the four radiative transfer prescriptions on the outflow velocity profile, we perform eight single star simulations adopting a low (10⁻⁸ M_⊙ yr⁻¹) and high (3 × 10⁻⁶ M_⊙ yr⁻¹) mass-loss rate. These mass-loss rates roughly cover the range observed in AGB stars (Habing & olofsson 2004; Ramstedt et al. 2008). Higher mass-loss rates where not considered, because in the Lucy and attenuation approximations a wind could not be launched with these parameters. The properties of the models are listed in Table 1, where we adopt the stellar parameters from Chen et al. (2020). In the geometrical, Lucy, and attenuation approximation, the wind injection velocity is set to v_inj = 33 km s⁻¹. With this initial velocity, the particles can reach the dust condensation radius, which is a necessary condition to launch a wind, and have reasonable terminal velocities in agreement with observations. For the free-wind approximation, which does not take into account dust condensation, an injection velocity of v_inj = 25.2 km s⁻¹ is chosen to match the terminal wind speed of the other cases.

Fig. 2 Velocity profiles from the SPH simulation (points) and 1D semi-analytical wind solution (lines) in a single-star configuration for the low (ṀAGB = 10−8 M⊙ yr−1, upper panel) and the high (ṀAGB = 3 × 10−6 M⊙ yr−1, bottom panel) mass-loss rate models. The free-wind approximation (red SPH points and 1D solid line profile), the geometrical approximation (light-green and dashed line), the Lucy approximation (orange and dash-dotted line), and the attenuation approximation (blue and dashed-double dotted line) are shown. The geometrical, Lucy, and attenuation approximation overlap in the upper panel. The vertical lines in the lower panel indicate the corresponding dust condensation radius for the geometrical, Lucy, and attenuation approximation. The dotted line shows the escape velocity.

3.1 Low mass-loss rate

The velocity profiles obtained in the low mass-loss rate regime (Ṁ_AGB = 10⁻⁸ M_⊙ yr⁻¹) are displayed in the upper panel of Fig. 2. In the free-wind approximation (red), the material is immediately accelerated outwards, and propagates at almost constant velocity. As Γ is set equal to one, no external force is applied, and only the pressure gradient close to the star influences the wind velocity. For the geometrical approximation (light-green), close to the wind injection zone, the material is too hot to condensate into dust (see Fig. 3, upper panel) and without the radiative acceleration, the gas pressure gradient is insufficient to drive the wind. Therefore, the velocity of the material initially decreases until it reaches the dust condensation radius (where T_eq = T_cond at R_Ge,dust = 3.6 au, see Fig. 3 upper panel). At this point, dust forms and the radiation force on the particles accelerates the material outwards. This results in a different velocity profile than the free-wind approximation. The velocity profiles in the Lucy and attenuation approximations (orange and blue) look identical to that of the geometrical approximation, because in the low mass-loss regime, densities in the wind are low. As a consequence, the optical depths τ_L (Eq. (23)) and τ (Eq. (28)) are very small, such that their effects are negligible.

Fig. 3 Radiative equilibrium temperature of the SPH particles (points) and 1D semi-analytical wind solutions (lines) in a single-star configuration for the low mass-loss (ṀAGB = 10−8 M⊙ yr−1, upper panel) and the high (ṀAGB = 3 × 10−6 M⊙ yr−1, bottom panel) mass loss-rates models. The free-wind approximation (red SPH points and 1D solid line profile), the geometrical approximation (light-green and dashed line), the Lucy approximation (orange and dash-dotted line), and the attenuation approximation (blue and dashed-double dotted line) are shown. The geometrical, Lucy, and attenuation approximation overlap in the upper plot. The vertical lines in the lower panel indicate the corresponding dust condensation radius for the geometrical, Lucy, and attenuation approximation.

3.2 High mass-loss rate

The velocity profiles obtained in the high mass-loss rate regime (Ṁ_AGB = 3 × 10⁻⁶ M_⊙ yr⁻¹) can be seen in the bottom panel of Fig. 2. In the free-wind and geometrical approximation (red and light-green), the equations of motion (Eqs. (1)–(3) are independent of the mass-loss rate, even with the Bowen cooling prescription (Eq. (7)) activated. However, this conclusion does not hold anymore if HI cooling is considered, because of the non-linear dependence of this rate on the density. In these singlestarmodels, no shock waves are present, the temperature remains below 3000 K, and cooling due to HI is inefficient. Therefore, the velocity profiles in the free-wind and geometrical approximation appear to be identical to the low mass-loss rate case (upper panel of Fig. 2). However, this is not the case for the Lucy and the attenuation approximation (orange and blue), because the density is sufficiently high that τ_L (Eq. (23)) and τ (Eq. (28)) become non-negligible, or even larger than the dilution factor W(r). This directly impacts the dust temperature profiles (see Fig. 3, bottom panel) and modifies the dust condensation radii (vertical lines in Figs. 2 and 3). An illustration of the optical depths profiles is presented in Figs. B.1 and B.2.

In the Lucy approximation, the dust temperature is increased compared to the geometrical case (Eq. (25)). Although τ_L is small, it is multiplied by, which can make its contribution significant. The condensation temperature is thus reached further out from the AGB star, at R_Lu,dust = 4.85 au (Fig. 3). For the attenuation approximation, the temperature decreases rapidly, owing to the e^−τ factor (Eq. (29)), and as a consequence, dust forms closer to the AGB star (at R_Ge,dust = 3.50 au). Thus, the radiation force is activated at higher velocities compared to the Lucy simulation. Once the wind material has passed the dust condensation radius, the radiation force becomes independent of the approximation, used to determine T_eq, since the dust opacity is now almost constant (K_d ≈ K_max). Beyond the condensation radius, the shape of the velocity profile in the Lucy approximation is similar to that of the geometrical approximation, but with a lower asymptotic value. This is different in the attenuation approximation. Here, not only does the dust temperature decrease exponentially with τ, but so does the radiation force (Eq. (27)). Thus, the acceleration of the material drops more rapidly, and the terminal wind velocity reaches lower velocities, yielding a flatter velocity profile in this case.

4 Binary-star models

An interesting application of the present implementation is studying the impact of a companion star or planet on the dynamics and morphology of the AGB wind. It is important that we investigate the applicability and differences of the four approximations in three-dimensional, non-spherically symmetric models, in which a binary companion is perturbing the outflow. We used the same setup as before, but add a companion with mass M_comp = 0.51 M_⊙ and an accretion radius R_acc = 0.1 au at an orbital separation of 6 au in a circular orbit (adopted from Chen et al. 2020). The simulations are evolved for a total of six orbital periods, reaching a quasi-steady state. As explained in Sect. 3, the Lucy and attenuation approximation give almost identical results as the geometrical approximation in the low mass-loss rate regime. Further, for the geometrical and free-wind approximation, the effect of changing the mass-loss rate is weak, hence we limit this discussion to the high mass-loss rate regime.

4.1 Wind-companion interaction strength

Maes et al. (2021) found that the ratio of the energy density of the companion to the kinetic energy of the wind gives a good indication of the complexity of the outflow and magnitude of the wind-companion interaction. This parameter is defined as $$ℰ=\frac{e_{\text{grav}}}{e_{\text{kin}}}=\frac{{\left(24G^3{M}_{\text{comp}}^2{M}_{\text{AGB}}\right)}^{1/3}}{{\upsilon }_{\text{w}}^{2}a\left(1-e\right)},$$(32)

where v_w is an estimate of the wind velocity at the location of the companion, which is calculated as $${\upsilon }_{\text{w}}=\sqrt{{\upsilon }_{\text{single}}^2\left(r=a\right)+{\upsilon }_{\text{AGB}}^2},$$(33)

with υ_single(r) the wind velocity at radius r in the corresponding single-star model, υ_AGB the orbital velocity of the AGB star, and e the eccentricity (here zero). The higher ε, the stronger the interaction of the companion with the wind, and the more complex the resulting outflow is expected to be. The only effect of our different approximations on this ε value, is the change in υ_single(r = 6 au), for which the four different values can be read from the bottom panel of Fig. 2. This velocity is ~26,20,18, and 12 km s⁻¹ in case of the free-wind, geometrical, attenuation, and Lucy approximation, respectively. This results in ε values of 0.4,0.7,0.8, and 1.7, respectively. This indicates that the interaction of the companion with the wind is expected to be the weakest for the free-wind approximation, and strongest for the Lucy approximation.

4.2 Wind structures

To illustrate the morphology and wind structures that result from the wind-companion interaction, Figs. 4 and 5 display the density and velocity maps, respectively, in slices through the orbital plane for the four approximations. The density profile in the meridional plane can be found in Fig. B.3. These figures reveal that the wind structures and global morphologies depend sensitively on the treatment of the radiative transfer.

In the free-wind approximation, the relatively weak wind-companion interaction strength (low ε < 1), creates a two-edged spiral structure attached to the companion, that shapes the wind into an approximate Archimedean spiral, as can be seen in the density profile in the orbital plane (Fig. 4, upper left). The creation of such Archimedean spirals is well studied and explained in detail by, for example, Malfait et al. (2021) and Maes et al. (2021). Due to the high wind velocity with respect to the orbital velocity, this thin high-density spiral propagates rapidly outwards with a radial velocity of about 27 km s⁻¹, which results in relatively wide inter-arm low-density gaps, with a width approximately equal to the distance travelled by the spiral structure in one orbital period (see Fig. 5, upper left). The meridional plane density distribution (Fig. B.3, upper left) shows that in the edge-on view, these spirals appear as thin, concentric arcs.

In the geometrical and attenuation approximation, the wind structure close to the stars is similar to the free-wind approximation, with a two-edged spiral structure attached to the companion (upper and lower right panels in Fig. 4). Although the wind velocity around the companion, and thereby the ε value, is similar for these approximations, the morphology of these systems is different. As the velocity in the geometrical approximation accelerates up to about ~27 km s⁻¹, reaching the same terminal velocity as the free-wind model (see Figs. 2 and 5), the 2-edged spiral remains again relatively thin, and the inter-arm separation large.

This is not the case for the attenuation approximation, where the wind material in the high-density spirals only reaches a velocity of ~20 km s⁻¹, and the low-density material in between the spirals has a velocity of ~10-15 km s⁻¹ (Fig. 5). This makes that the spiral structure more compressed, with smaller low-density inter-arm gaps. Because there is a velocity dispersion within the spiral, the outer frontward spiral edge has a higher radial velocity than the inner backward spiral edge, such that after ~1.75 orbital periods the outer spiral edge catches up and interacts with the previous inner spiral edge, that originated one orbital period earlier (around x = 0 au, y = 80 au). After this interaction, one approximate Archimedean spiral structure remains, and the inter-arm low-density gaps disappear (Malfait et al. 2021; Maes et al. 2021). In the meridional plane, the spiral appears again as arcs (Fig. B.3, lower right). This plot also shows that the high-density structure is more compressed, with smaller low-density inter-arc regions, and that the outer edge of the widening arcs catches up and interacts with the previous inner edge (overlap first visible at x = 100 au, y = 0 au).

Due to the larger dust condensation radius in the Lucy approximation, the wind velocity around the location of the companion is significantly lower than in the previously discussed models (bottom panel in Fig. 2, at r = 6 au, and r < 10 au region in Fig. 5). This allows the companion to compress more wind material around it, such that instead of a 2-edged spiral structure, there is one spiral originating behind the companion, and a second 2-edged bow shock spiral originating in front of the companion (Malfait et al. 2021; Maes et al. 2021). This can be seen in more detail in the zoomed-in density plot in the upper left plane of Fig. 6. Moreover, due to the strong compression of gas around the companion, as well as sufficient cooling to reduce the thermal pressure, an accretion disk has formed. In the Lucy approximation, there is no radiative force active close to the companion (see Sect. 4.3), so the forming accretion disk is not blown away, facilitating its formation. This accretion disk is shown in more detail in Fig. B.4, where the density distribution is overplotted with velocity vectors. Figure B.4 displays how material spirals in towards the companion sink particle through a high-density disk. For a more elaborate description of accretion disks, see Lee et al. (2022) and Malfait et al. (in prep). In the meridional plane, the bow shock spiral translates into an expansion of the edge-on arcs (Fig. B.3, lower right).

Chen et al. (2020) used the attenuation approximation in their simulations, and report the formation of both an accretion and a circumbinary disk. The absence of these features in our computations stems from the fact that the terminal wind velocity and cooling prescriptions are different between these two works. In our simulations, the terminal velocity is ~20km s⁻¹, higher than the value of ~15km s⁻¹ found by Chen et al. (2020). With a faster wind, the particles interact less with the companion (Sect. 4.1) and pass over the companion preventing the formation of the disk. The disk may, however, become visible when decreasing the accretion radius of the companion (Lee et al. 2022). The same is true for the circumbinary disk, as a higher terminal wind velocity prevents the formation of such structures. Chen et al. (2020) also include molecular cooling, associated with H₂, H₂O, and Co, which contribute to reduce the heat (pressure) and to provide more favourable conditions for gas condensation. These processes mostly influence the formation of circumbinary disks, since in a circumstellar accretion disk, the temperature is higher, and atomic cooling is expected to dominate (Mastrodemos & Morris 1999; see also Woitke et al. 1996 for information on the dominant cooling processes).

Fig. 4 Density distributions in a slice through the orbital plane for the four simulations, each of which using different radiative transfer prescription: free-wind (top left), geometrical (top right), Lucy (bottom left), and attenuation (bottom right) approximation, for the high mass-loss rate case with a binary companion. Both stars are on the x-axis, where the primary AGB star is on the left, and the companion on the right.

Fig. 5 Same as Fig. 4, but for the velocity distribution.

4.3 Impact of the approximations on the radiation force

In the free-wind approximation, the radiation force is not explicitly calculated, but is assumed to be equal to the gravitational force of the AGB star. For the geometrical approximation, its expression does not depend on the presence of a companion (Eqs. (19) and (21)). This is, however, not the case for the Lucy and attenuation approximation, because of the directional dependence of τ_L (Eq. (23)) and τ (Eq. (28)), that enter the evaluation of T_eq and J. The density distribution ρ, optical depths (τ and τ_L), dust temperature T_eq, wind velocity υ, opacity κ = κ_d + κ_g, and Eddington factor Γ in the orbital plane are shown in Figs. 6 and 7 for the Lucy and attenuation approximations, respectively.

In the Lucy simulation, close to the companion, inside the 2-edged bow shock, the density and τ_L are high (upper middle panel). The mean intensity is thus increased in that region (Eq. (22)) and under the condition of LTE, the equilibrium temperature is locally higher (Eq. (16)). This brakes the symmetry and the dust condensation surface, inside which no dust forms and which is defined as the region where T(r) = T_cond, is not spherical anymore. This 3D surface is shown as the black 2D contour in the various panels. The opacity increases rapidly across this surface, which reflects directly on the Eddington factor Γ (bottom right panel).

These features are different in the simulation with the attenuation approximation (Fig. 7). Because the dust condensation radius is significantly smaller than the orbital separation, the dust condensation surface remains approximately spherically symmetric, as shown by in the various panels. The effect of this attenuation is clearly visible behind the companion, where the optical depth starts to deviate from spherical symmetry. This effect can be seen as a shadow behind the companion (x > 4 au side along y = 0 au). This shadow is cast as the material close to the companion forms a higher-density spiral structure, that absorbs the radiation from the AGB star. This effect also reduces the dust temperature in the region behind the companion (upper right plot). As the dust temperature is already sufficiently low for dust to form, this does not influence the opacity. But since the equation for the radiation force (described by Γ, see Eq. (27)) contains a factor e^−τ, the radiation force completely vanishes behind the companion.

Fig. 6 Relevant properties of the high mass-loss rate binary simulation using the Lucy approximation. The density is plotted in the upper left panel, the Lucy optical depth in the upper middle panel, the dust temperature in the top right panel, the velocity in the lower left panel, the opacity in the lower middle panel, and the Eddington factor in the lower right panel. The thin solid black contour indicates the location of the dust condensation surface.

Fig. 7 Same as Fig. 6 but for the attenuation approximation. τL is replaced by τ in the upper middle panel.

5 Discussion

5.1 Accuracy of the ray-tracing approximations

To gauge the accuracy and quantify the benefit of our raytracer implementation, we compare the dust temperatures T_eq and the radiation forces, obtained with our prescriptions, to the results obtained with the 3D ray-tracing radiative transfer solver MAGRITTE. We focus only on the Lucy and attenuation approximation, as those are the only two prescriptions that leverage the newly implemented ray-tracer. To make this comparison, we take the final snapshots of their respective simulations (discussed in Sect. 4.3), and use these as an input for MAGRITTE.

Although MAGRITTE is usually advertised as a line radiative transfer code (De Ceuster et al. 2020a,b), we only use its core ray-tracer and solver functions here. More specifically, we supply it with the grey opacity (Eq. (6)) and emissivity η = kρJ (with J defined as in Eq. (16)), and let it solve the radiative transfer equation (Eq. (8)) along 2 700 uniformly distributed rays, originating from each SPH particle. From the resulting intensities along those rays, we then derive the mean intensity J and flux F. However, the emissivity depends on the dust temperature, which itself depends on the radiation field. Therefore, to obtain a self-consistent solution between the dust temperatures and the radiation field, we need to compute them in an iterative way. Starting from the analytic radiative equilibrium temperature of the geometrical approximation (Eq. (21)), the dust opacities are computed (Eq. (6)), and then MAGRITTE is used to compute the mean intensity J. From this, the dust temperatures can then be recomputed using (Eq. (16)). This process is repeated until the change in the dust temperatures becomes negligible (the resulting mean relative temperature differences in the final iteration is 0.02%, and the maximal 2%). To reduce the computational cost, only the SPH particles within a radius of 30 au are included. This should not alter the results, as densities (and thus interactions) are diluted significantly beyond this radius.

5.1.1 Lucy approximation

Figure 8 shows the dust temperature T_eq, calculated using MAGRITTE, as a function of distance from the AGB star. Here, the dust temperature in the geometrical approximation (used as initial temperature in MAGRITTE) is shown in magenta, results obtained with the Lucy approximation (as calculated in PHANTOM) are shown in orange, and results from MAGRITTE are shown in green. We see that inside the orbit (r ≲ 6 au), the 1D Lucy dust temperature nicely follows the MAGRITTE prediction. Just before the location of the companion, a sudden drop in the MAGRITTE dust temperature appears, as well as a slight increase at the top of the green curve (at r = 5. 5 and 6 au). The increase coincides with the edge of the accretion disk and spiral arm, which are heated more efficiently. In the Lucy simulation, this effect is much more pronounced, due to the underlying assumption of spherical symmetry, such that, whenever a direction with high optical depth is encountered, an entire sphere at this optical depth is assumed. The particles from the MAGRITTE post-processed model with a low dust temperature, are the particles located just behind the heated regions. The second drop in the dust temperature around r ≈ 6–6.3 au occurs behind the companion, and is due to the bow shock acting as a shadow (see also Fig. 9).

Figure 9 shows the difference between the dust temperature computed in the Lucy approximation (Fig. 6, upper right panel) and the one obtained with MAGRITTE, for a slice through the orbital plane. While there are clear deviations from the Lucy temperature in the direction of the companion, in other directions, where the density profile is less perturbed, differences with respect to the Lucy temperature remain small. Hence, the Lucy approximation performs best in regions where the underlying assumption of spherical symmetry remains approximately valid. In the direction of the companion, and especially in the accretion disk around the companion (x = 4 au, y = 0 au), the dust temperature is artificially heated, due to the underlying assumption of spherical symmetry, as explained above. Behind the companion, the shadow is not captured in the Lucy approximation, resulting in the deviations at x > 5 au and y < 0 au (red region).

The analysis of the radiation force (ƒ_rad) is more complex, as this is a vector quantity. In all four prescriptions, the radiation force is assumed to be radial, thus pointing from the AGB star to the position of the SPH particle under consideration. To verify whether this is also true for the radiation flux computed with MAGRITTE, the non-radial fraction of the radiation force (cos (θ)) can be computed $$\cos \left(\theta \right)=\frac{f_{rad}\cdot r}{\Vert f_{\text{rad}}\Vert \Vert r\Vert }.$$(34)

This quantity is shown as a function of distance from the AGB star in Fig. 10. A clear spike is visible at the location of the companion (r = 6 au), due to the dense accretion disk emitting a significant amount of radiation, which dominates the radiation force near the accretion disk. Although these forces are highly non-radial, they are relatively small in magnitude when compared to the local gravitational attraction of the AGB star. This quantity, $$\frac{\cos \left(\theta \right)f_{\text{rad}}}{f_{\text{grav}}}=\frac{\kappa \left(T_{\text{eq}}\right)\cos \left(\theta \right)\Vert F\Vert /c}{G{M}_{\text{AGB}}/r^2},$$(35)

is displayed in Fig. 11. This shows that the assumption of a radial radiation force is reasonable. In the remainder of this analysis, we assume a radial radiation force and only consider its magnitude. The small spike around the inner boundary at r = R_AGB = 1.24 au is a numerical artefact of MAGRITTE that reveals the discretization of the (spherical) stellar surface. This feature is also present in Fig. 16, but it is inconsequential.

The Eddington factor Γ (Eq. (5)) is shown for a slice through the orbital plane in Fig. 12. First, we compare Γ resulting from MAGRITTE (Fig. 12) to Γ resulting from the Lucy approximation (Fig. 6, lower right panel). In the MAGRITTE case, the dust condensation surface is a perfect sphere, and since it is smaller than the orbital separation, it is not perturbed by the companion, and remains spherically symmetric. This is in contrast to the dust condensation surface in Fig. 6, which is extended beyond the companion and even engulfs it. The radiation force inside the accretion disk around the companion is negligible (Γ ~ 0), since the high local densities make the radiation field isotropic, such that the contributions to the radiation force from different directions cancel each other out. Behind the companion ( x > 5 au and y < 0 au), there is a low-Γ region, and the resulting shadow looks similar to the shadow cast in that attenuation simulation (Fig. 7, lower right panel). The radiation force increases again when moving further out, due to the fact that material above and below the shadow in the orbital plane shine and accelerate the material again.

Fig. 8 Radiative equilibrium temperature as a function of distance from the primary AGB star for the simulation using the Lucy approximation. Magenta represents the temperature calculated using the geometrical approximation, orange using the Lucy approximation, and green using MAGRITTE. The lower envelope of the Lucy simulation (orange points) which is not visible as it falls behind the green points, closely follows the Teq,Lucy dot-dashed line.

Fig. 9 Difference in the orbital plane between the radiative equilibrium temperature calculated with the Lucy approximation and MAGRITTE. The temperature is similar in the two cases except in the directions of the companion where the Lucy approximation yields higher temperatures.

Fig. 10 2D histogram of the relative non-radial component of the radiation force, as a function of distance from the AGB star, given by MAGRITTE for the snapshot of the model using the Lucy prescription. The colourbar is in log-scale.

Fig. 11 2D histogram of the non-radial component of the radiation force, relative to the gravitational attraction of the AGB star (see Eq. (35)) as a function of distance from the AGB star, given by MAGRITTE for the snapshot of the model using the Lucy prescription. The colourbar is in log-scale.

Fig. 12 Eddington factor Γ (Eq. (5)), calculated using the magnitude of the flux obtained with MAGRITTE, in a slice through the orbital plane for the snapshot of the model, using the Lucy prescription.

5.1.2 Attenuation approximation

Figure 13 shows the dust temperature T_eq, calculated using MAGRITTE as a function of distance from the AGB star. The MAGRITTE results for the snapshot of the simulation, using the attenuation prescription, show the same pattern as in the Lucy approximation. Within the orbit, the MAGRITTE dust temperature follows closely the 1D Lucy approximation (and not the attenuation approximation), and at the location of the companion a drop appears due to the interaction of the high density spiral arm. In contrast to the Lucy approximation, there is only one drop and one peak, as there is only one spiral arm in this simulation, and no accretion disk or bow shock are formed. In the attenuation approximation, the dust temperature decreases faster than in the geometrical and Lucy approximation, similar to the 1D profile (Fig. 3). This is caused by the fact that the attenuation approximation only accounts for absorption and not for re-emission, resulting in an underestimation of the intensities and dust temperatures. Looking at the region behind the companion, both MAGRITTE and the attenuation approximation show a drop in the dust temperature. This drop starts at slightly lower radii in the attenuation approximation, compared to the MAGRITTE calculation. In the attenuation approximation, the shadow immediately forms when high densities are encountered, as the optical depth increases. In a full radiative transfer treatment, when a high density region is encountered, the first layers absorb a lot of photons, locally trapping some of the radiation and producing a local heating that is seen in the peak in T_eq,Magritte at 6 au. After this peak, the temperature drops at ≈7 au, as radiation escapes isotropically. The temperature decrease in the attenuation approximation is also too strong, and this is caused by the radiation being only blocked, while in MAGRITTE, re-emission is accounted for.

Figure 14 shows the difference between the dust temperature, resulting from the attenuation approximation (see Sect. 2.2.4; Fig. 7, upper right panel), and the one obtained with MAGRITTE, for a slice through the orbital plane. Close to the AGB star, the difference in dust temperature remains small, but the differences increase rapidly further out. In the shadow region behind the companion (x > 5 au and y < 0 au), the temperature in the attenuation approximation is too cold, but this scheme is able to reproduce the shadow, but the effect is exacerbated.

The non-radial fraction of the radiation force, shown in Fig. 15, is globally lower than in the case of the Lucy approximation, mainly because no accretion disk is forming in the attenuation simulation, and thus there is no region where the photons are trapped. When normalizing the radiation force to its maximum value (see Fig. 16), we see a small non-radial contribution, similar to the Lucy case This renders the conclusion that the assumption of a radial radiation force is reasonably good.

The Eddington factor Γ (Eq. (5)) in the orbital plane is shown in Fig. 17. Due to the higher temperature estimate in MAGRITTE, the dust condensation radius is shifted farther out than in the attenuation approximation (Fig. 7), resulting in a lower radiation force close to the AGB star. Further out and away from the shadow, the Eddington factor remains approximately constant according to MAGRITTE, while it is monotonically decreasing in the attenuation approximation. Thus, the global radiation force is not correctly described by the attenuation approximation, for the same reason as the dust temperature. The shadow behind the companion (x > 5 au and y < 0 au) is reproduced by the attenuation approximation. However, due to the weak radiation force in this region, it is unclear whether this feature is modelled accurately or overcompensated.

In summary, for regions that (locally) resemble a spherically symmetric wind (i.e. away from the companion), the Lucy approximation provides a correct estimate for the dust temperature and radiation force, taking into account the radiation that is absorbed and re-emitted isotropically. This effect heats up the material compared to the geometrical approximation. Close to the companion, dust is artificially heated due to the spherical symmetric assumption underlying the Lucy approximation. The Lucy approximation fails to account for the shadow cast behind the companion. On the other hand, in the attenuation approximation, the regions away from the companion are not modelled correctly. The material is artificially cooled, causing a reduction of the radiation force. Close to the companion, a shadow forms, which is captured by the attenuation approximation both in the dust temperature as well as in the radiation force. However, the drastic artificial decrease of the radiation force, and the fact that most of the morphological simulation is spherically symmetric, favours the Lucy approximation. This suggests that a combination of the Lucy and attenuation approximation might yield even better results, combining the strengths of both approximations.

Fig. 13 Radiative equilibrium temperature as a function of distance from the primary AGB star, for the simulation using the attenuation prescription. Magenta represents the temperature calculated, using the geometrical prescription, orange using the attenuation prescription, and green using MAGRITTE.

Fig. 14 Difference in the orbital plane between the radiative equilibrium temperature, calculated with the attenuation prescription MAGRITTE. Significant differences persist over the entire domain.

Fig. 15 2D histogram of the relative non-radial component of the radiation force, as a function of the distance from the AGB star, given by MAGRITTE, for the snapshot of the simulation using the attenuation prescription. The colourbar is in log-scale.

Fig. 16 2D histogram of the non-radial component of the radiation force, relative to the gravitational attraction of the AGB star (see Eq. (35)), as a function of distance from the AGB star, given by by MAGRITTE, for the snapshot of the simulation using the attenuation prescription. The colourbar is in log-scale.

Fig. 17 Eddington factor Γ (Eq. (5)), calculated using the magnitude of the flux vector, in a slice through the orbital plane, using MAGRITTE for the snapshot of the simulation, using the attenuation prescription.

5.2 Future work

One of the ingredients still missing in this study is a treatment of pulsations. In our approach, pulsations are neglected and material is launched at a preset velocity from the stellar surface. In future studies, pulsations will be modelled by simulating a radially oscillating inner boundary acting as a piston. This method was originally used in the 1D study by Bowen (1988), but has recently been implemented in the SPH context by Aydi & Mohamed (2022). After including pulsations and adopting the Lucy approximation from this study, the full Bowen (1988) study can be replicated in 3D, and include a companion.

Furthermore, the treatment of the dust opacity can still be improved. Siess et al. (2022) already implemented carbon dust formation in PHANTOM, which can be used in combination with this study. Using this formalism, the dust opacities can be made more consistent with the physical and chemical environments of stellar outflows, which will result in a more accurate description of the radiation force.

An accurate treatment of both chemical processes and cooling are of crucial importance as well (e.g. Boulangier et al. 2019; Van de Sande & Millar 2019). The chemistry can alter the poly-tropic index, as well as the mean molecular weight, which now are assumed to be constant. To account for the complex chemistry taking place in AGB outflows, without trading too much of the required computation time, machine learning techniques will be used to emulate the chemical network (e.g. Holdship et al. 2021; Grassi et al. 2022), reducing the computation time to allow for an on-the-fly simulation of the chemistry (Maes et al. in prep). Furthermore, we could improve on the cooling prescriptions. Cooling has a significant impact on the energy equation, since insufficient cooling can prevent the formation of accretion disks around the companion (Theuns & Jorissen 1993; Malfait et al., in prep.). Moreover, it can alter the transfer rates of both mass and angular momentum.

A final improvement for the simulations itself is the gas-dust interaction. In this study position coupling is assumed, although in a dust-driven wind the gas is dragged by the faster moving dust (Mattsson & Sandin 2021), which might result in significant drift velocities. In some cases this drag can lead to a decrease in both the mass-loss rate as well as the wind velocity (Sandin & Mattsson 2020), which in turn influences the morphology as well. In PHANTOM, two approaches are already implemented to account for this drag. In the first approach, the two fluids are treated separately (Laibe & Price 2012), while in the second, both fluids can be combined using the so-called one-fluid approximation (Laibe & Price 2014a,b; Price & Laibe 2015).

Once truly realistic models are acquired, these forward models can be used to compare to observations. This can be done using radiative transfer codes like MCFOST (Pinte et al. 2006; Tessore et al. 2021) or MAGRITTE (De Ceuster et al. 2020a,b) to produce synthetic observations. These synthetic observations are constructed by tracing chemical species in the simulation, and creating (synthetic) spectral line maps. These synthetic observations can then be compared to real observations, unraveling the origin of the complexities of these winds.

6 Summary

In this paper, we present the implementation of a ray tracer for radiative transfer in the SPH code PHANTOM. Using this ray tracer, a 3D map of the optical depth can be obtained, a key quantity that is needed, for example, in AGB wind simulations. This technique enables us to investigate four different radiative transfer approximations: the free-wind approximation (no radiative transfer), the geometrical approximation (including geometrical dilution, as well as dust formation), the Lucy approximation (including radiation not exclusively coming from the AGB star), and the attenuation approximation (including the absorption of the stellar flux without re-emission).

We performed simulations of a single AGB star for each of the four prescriptions, with a low and high mass-loss rate, resulting in eight different models. The velocity profile, resulting from the free-wind approximation, remains relatively constant throughout the entire spatial domain. In the other approximations, where the radiation force and dust formation are accounted for, there is an initial decrease in velocity close to the AGB star, since this region is too hot for dust to condensate. Beyond the dust condensation radius, the material is accelerated outwards, allowing the velocity to escape the escape value. In the low mass-loss rate regime, the effects of the Lucy and attenuation approximations are small and the velocity profiles are almost identical. In the high mass-loss regime, the effects and their differences become more pronounced. For the Lucy approximation, the dust condensation radius moves outwards, resulting in a lower terminal velocity compared to the geometrical approximation. For the attenuation approximation, the dust condensation radius is shifted inwards, activating the radiation force earlier in the wind, but as radiation is attenuated, the radiation force also decreases, resulting in a flat velocity profile.

We also considered binary systems and investigated the effect of the different radiative transfer approximations on the wind morphology. In the free-wind approximation, a thin, two-edged spiral structure forms, while in the geometrical approximation, a thicker spiral is present, with increased interaction of the wind close to the companion. Again, in the low mass-loss rate regime the Lucy and attenuation approximation follow the geometrical approximation, while for the high mass-loss rate regime differences appear. In the Lucy approximation, an accretion disk is able to form around the companion because of the lower wind velocity. This accretion disk creates an additional bow shock, increasing the morphological complexity. In the attenuation approximation, the two-edged spiral is compressed significantly, such that the spiral arms interact, forming a single Archimedean spiral structure.

In order to gauge the accuracy of these radiative transfer approximations, and hence the applicability, we compared the dust temperature and the radiation force from the simulations, using the Lucy and attenuation approximation, to results obtained with the full 3D radiative transfer code MAGRITTE. We showed that the non-radial component of the radiation force is small, which implies that the often made assumption of a radial radiation force is adequate. The Lucy approximation can reproduce the dust temperature and the radiation field accurately in the parts of the simulation that resemble a spherically symmetric outflow. However, in the direction of the companion, the density is higher and the radiation field is overestimated due to the underlying assumption of spherical symmetry in the Lucy approximation. The Lucy approximation can model the radiation force correctly in regions that resemble a spherically symmetric outflow, but this approximation does not reproduce the shadow cast behind the companion because of excessive re-emission. In the attenuation approximation, the dust temperature and the radiation force are underestimated, because this approximation only models the extinction of radiation, but not the re-emission. Although, as expected, the attenuation approximation does account for the shadow cast behind the companion. In conclusion, the Lucy approximation turns out to be the most adequate radiative transfer prescription for AGB binary simulations, since most of the domain in the considered simulations resembles a spherically symmetric outflow, and the acceleration close to the AGB star is important for the companion interaction.

Acknowledgements

M.E., J.M., S.M. and L.D. acknowledge support from the Research Foundation Flanders (FWO) grant G099720N. L.S. is a senior FRS-F.N.R.S research associate. F.D.C is a Postdoctoral Research Fellow of the Research Foundation – Flanders (FWO), grant 1253223N. T.K. and L.D. acknowledge support from the KU Leuven IDN grant IDN/19/028. T.C. is a PhD Fellow of the Research Foundation – Flanders (FWO), grant 1166722N.

Appendix A Ray tracer

In order to investigate and validate the trade-offs made during the development of the ray-tracer (see Sect. 2.3), we compare its performance against a second ray tracer algorithm that is optimized for accuracy, but not for computational speed. This ray tracer would be too slow to be used on-the-fly in PHANTOM simulations, but is useful to gauge the performance of our implementation. In this appendix, we gradually build up the reasoning that led to the ray tracer that was eventually implemented in PHANTOM, starting from this ideal (but slow) ray tracer (IRS; ideal ray-tracing scheme).

The ideal ray tracer works as follows. Starting from a specific location in the simulation, it calculates the nearest neighbours, out of which it will select the next particle to be considered on the ray. There, it recalculates the nearest neighbours, to again decide which particle should be considered next, and, as such, moves along the ray. To minimize the step size, and obtain the optimal discretization for integrals along the ray, the smallest set of neighbours is used. This set is obtained using a Delau-nay tetrahedralization, which is used in MAGRITTE as well (De Ceuster et al. 2020b). While walking through the selected particles along a ray, these particles are projected onto the ray, and, at these projected points, the integrand is evaluated (using Eq. 31). To calculate the optical depth, τ_IRS, in this ideal scenario, a ray is traced starting from the particle in question, in the direction of the AGB star. This ensures that each of the rays properly represents the actual integral that is calculated. The resulting algorithm is similar to Kessel-Deynet & Burkert (2000).

To speed up the calculation, while trying to remain as accurate as possible, different improvements are investigated in the rest of this Section². To quantify the performance, the optical depth, τ, calculated with our improved scheme, is compared to its IRS value, τIRS, by calculating a relative error $${\delta }_{\text{rel}}=\frac{1}{N_{\text{part}}}{\displaystyle \sum _{i=1}^{N_{\text{part}}}\left|\frac{{\tau }_{\text{IRS,i}}-{\tau }_i}{\left({\tau }_{\text{IRS,i}}+{\tau }_i\right)/2}\right|.}$$(A.1)

To test the algorithm, we used an SPH model of an AGB binary, more specifically, dump 600 of model v05e50 (υ_inj = 5km s⁻¹, e = 0.5 using the free-wind approximation) of Malfait et al. (2021). This is a late snapshot of the most complex morphology in that study, using ≈ 1.2 × 10⁶ SPH particles.

A.1 Nearest neighbours

A first bottleneck in the IRS algorithm, is the computation of the Delaunay nearest neighbours. Since PHANTOM already works with a kd-tree to store its particles (a tree-like data structure that significantly speeds up nearest neighbour determinations), calculating nearest neighbours leveraging this tree is much faster. The set of neighbours includes ~ 60 particles, such that the sampling along the ray can take larger steps, and hence the integration is done more crudely. Using this nearest neighbour set and comparing the results to the IRS, using Eq. A.1 results in a relative error of only 1%, while giving a speed-up of a factor 2. This speedup does not yet include the calculation of the Delaunay nearest neighbours (taking a significant amount of time), such that the actual speed-up is even higher.

Fig. A.1 Relative error (Eq. A.1) and computation time, required for calculating τ as a function of the number of rays traced in the reference model. The solid red line represents the relative error using the SPH nearest neighbours, blue the computation time, and green the number of rays in the IRS case.

A.2 Ray directions

Tracing rays inwards is inefficient, since plenty of SPH particles will be passed several times, especially those close to the star. To avoid this, one can reverse the ray tracing by starting at the stellar surface, and pointing the rays outwards until the edge of the simulation is reached, instead of tracing a ray from each particle to the star. The rays that originate from the star, should be traced uniformly outwards. The directions of these rays can be determined using HEALPIX³ (Górski et al. 2005). HEALPIX subdivides the 2-sphere into iso-laterally distributed equal area pixels, such that the centre of each pixel can be used as the direction of a ray. The scheme starts with 12 pixels, so-called HEALPIX-order 0, and the amount of rays can be increased by subdividing each cell into four, increasing the order o by 1. This results in n_rays = 12 × 4^o. Once the rays are traced, the information along the rays needs to be mapped back to the SPH particles. Using HEALPIX, the closest ray to a particle can be found easily by leveraging its smart pixel positioning. Once the nearest ray is found, the particle is projected on the ray, and the integral is evaluated there by linear interpolation between the two closest known values along the ray.

The performance of this algorithm strongly depends on the HEALPIX order. The relative error and the required computation time are shown as a function of the number of rays in Fig. A.1. There is an improvement in computation time, as long as the number of rays is less than the number of particles in the simulation. This is, however, at the expense of the relative error. For HEALPIX order 5, for instance, the computation time goes down with a factor 250, while inducing a relative error of 2%.

A.3 Ray interpolation

To further reduce the 2% relative error of the algorithm, we investigate the interpolation from the rays onto the particles. Instead of using only the information of the ray nearest to a particle, an interpolation between multiple rays can be used. We consider either the four or nine closest rays, as well as different interpolation exponents, k, in the interpolation formula, given by $$I_i=\frac{1}{{\displaystyle {\sum }_{j=1}{r}_j^{-k}}}{\displaystyle \sum _j\frac{I_j}{r_j^k},}$$(A.2)

where I_i is integral at the location of the i’th particle, I_j is the integral at the projection of the i’th particle on the j’th ray, and r_j the distance between these two. Here, j is a sum over the four or nine closest rays, and k specifies the interpolation exponent. The relative error as a function of the number of rays, used in the interpolation for HEALPIX order 5, is shown in Fig. A.2. There is a significant improvement when moving from one to four rays, as the latter interpolation can account for differences in between rays. However, going up to nine rays does not significantly improve the accuracy, as most of the smoothing already happens in between the four closest rays. The relative error even increases slightly when using k = 1, because the interpolation smooths out large variations. Since in our simulations large variation in density, and thus optical depth, are to be expected, using nine rays and k = 1 will not be ideal. It turns out that for a interpolations exponent of k = 2 and four rays is the best combination for the specific model that we considered.

Fig. A.2 Relative error (Eq. A.1) as a function of the number of rays used in the interpolation (using HEALPIX order 5).

The resulting relative error as a function of computation time for this interpolation is shown in Fig. A.3. For low HEALPIX orders, the interpolation takes a significant amount of time, while only slightly improving the relative error. This is because the initial calculations along the few rays, that are traced, cannot capture all the complexities in between these rays. Going to higher orders (starting from HEALPIX order 3) the interpolation is worthwhile. The increase in computation time becomes ever smaller, as the extra time used for the interpolation only scales with the number of particles, and not with the number of rays. Using HEALPIX order 5, there is a speedup of a factor 250 compared to the IRS, while obtaining a relative error of only 1.5%.

Due to this significant speed-up, the calculation of the integrals using HEALPIX order 5 only take about 10% of the computation time of a normal SPH hydro timestep (or shorter for lower orders). Hence, the extra computation time to include these calculations is practically feasible to perform as on-the-fly calculations.

Fig. A.3 Relative error (Eq. A.1) as a function of computation time for calculating τ. Different dots represent different HEALPIX orders, where the blue represents no interpolation, and the orange represents the interpolation using four rays and k = 2.

Appendix B Additional figures

Fig. B.1 Lucy optical depth, τL (Eq. 23), as a function of radial distance in the single-star simulation, for the high mass-loss rate (3 × 10−6 M⊙ yr−1) model.

Fig. B.2 Optical depth, τ (Eq. 28), in the attenuation prescription as a function of radial distance, for the high mass-loss rate (3 × 10−6 M⊙ yr−1) model.

Fig. B.3 Density distribution in a slice through the meridional plane for the four prescriptions: free-wind (top left), geometrical (top right), Lucy (bottom left), and attenuation (bottom right) prescription, for the high mass-loss rate case with a binary companion. Both stars are on the x-axis, where the primary AGB star is on the left, and the companion on the right.

Fig. B.4 Density distribution in a slice through the orbital plane close to the companion in the Lucy approximation, revealing the presence of an accretion disk. The small arrows indicate the direction of the velocity field.

References

All Tables

All Figures

	Fig. 1 Visual representation of the ray-tracing algorithm. Starting from the upper right corner at point P1 and following the direction of the ray, the subsequent point P2 is determined within the sphere of influence of radius Rkernh. The procedure is then repeated until the boundary of the domain is reached (see text for explanations).
In the text

Fig. 2 Velocity profiles from the SPH simulation (points) and 1D semi-analytical wind solution (lines) in a single-star configuration for the low (ṀAGB = 10−8 M⊙ yr−1, upper panel) and the high (ṀAGB = 3 × 10−6 M⊙ yr−1, bottom panel) mass-loss rate models. The free-wind approximation (red SPH points and 1D solid line profile), the geometrical approximation (light-green and dashed line), the Lucy approximation (orange and dash-dotted line), and the attenuation approximation (blue and dashed-double dotted line) are shown. The geometrical, Lucy, and attenuation approximation overlap in the upper panel. The vertical lines in the lower panel indicate the corresponding dust condensation radius for the geometrical, Lucy, and attenuation approximation. The dotted line shows the escape velocity.

In the text

Fig. 3 Radiative equilibrium temperature of the SPH particles (points) and 1D semi-analytical wind solutions (lines) in a single-star configuration for the low mass-loss (ṀAGB = 10−8 M⊙ yr−1, upper panel) and the high (ṀAGB = 3 × 10−6 M⊙ yr−1, bottom panel) mass loss-rates models. The free-wind approximation (red SPH points and 1D solid line profile), the geometrical approximation (light-green and dashed line), the Lucy approximation (orange and dash-dotted line), and the attenuation approximation (blue and dashed-double dotted line) are shown. The geometrical, Lucy, and attenuation approximation overlap in the upper plot. The vertical lines in the lower panel indicate the corresponding dust condensation radius for the geometrical, Lucy, and attenuation approximation.

In the text

Fig. 4 Density distributions in a slice through the orbital plane for the four simulations, each of which using different radiative transfer prescription: free-wind (top left), geometrical (top right), Lucy (bottom left), and attenuation (bottom right) approximation, for the high mass-loss rate case with a binary companion. Both stars are on the x-axis, where the primary AGB star is on the left, and the companion on the right.

In the text

	Fig. 5 Same as Fig. 4, but for the velocity distribution.
In the text

Fig. 6 Relevant properties of the high mass-loss rate binary simulation using the Lucy approximation. The density is plotted in the upper left panel, the Lucy optical depth in the upper middle panel, the dust temperature in the top right panel, the velocity in the lower left panel, the opacity in the lower middle panel, and the Eddington factor in the lower right panel. The thin solid black contour indicates the location of the dust condensation surface.

In the text

	Fig. 7 Same as Fig. 6 but for the attenuation approximation. τL is replaced by τ in the upper middle panel.
In the text

Fig. 8 Radiative equilibrium temperature as a function of distance from the primary AGB star for the simulation using the Lucy approximation. Magenta represents the temperature calculated using the geometrical approximation, orange using the Lucy approximation, and green using MAGRITTE. The lower envelope of the Lucy simulation (orange points) which is not visible as it falls behind the green points, closely follows the Teq,Lucy dot-dashed line.

In the text

	Fig. 9 Difference in the orbital plane between the radiative equilibrium temperature calculated with the Lucy approximation and MAGRITTE. The temperature is similar in the two cases except in the directions of the companion where the Lucy approximation yields higher temperatures.
In the text

	Fig. 10 2D histogram of the relative non-radial component of the radiation force, as a function of distance from the AGB star, given by MAGRITTE for the snapshot of the model using the Lucy prescription. The colourbar is in log-scale.
In the text

	Fig. 11 2D histogram of the non-radial component of the radiation force, relative to the gravitational attraction of the AGB star (see Eq. (35)) as a function of distance from the AGB star, given by MAGRITTE for the snapshot of the model using the Lucy prescription. The colourbar is in log-scale.
In the text

	Fig. 12 Eddington factor Γ (Eq. (5)), calculated using the magnitude of the flux obtained with MAGRITTE, in a slice through the orbital plane for the snapshot of the model, using the Lucy prescription.
In the text

	Fig. 13 Radiative equilibrium temperature as a function of distance from the primary AGB star, for the simulation using the attenuation prescription. Magenta represents the temperature calculated, using the geometrical prescription, orange using the attenuation prescription, and green using MAGRITTE.
In the text

	Fig. 14 Difference in the orbital plane between the radiative equilibrium temperature, calculated with the attenuation prescription MAGRITTE. Significant differences persist over the entire domain.
In the text

	Fig. 15 2D histogram of the relative non-radial component of the radiation force, as a function of the distance from the AGB star, given by MAGRITTE, for the snapshot of the simulation using the attenuation prescription. The colourbar is in log-scale.
In the text

	Fig. 16 2D histogram of the non-radial component of the radiation force, relative to the gravitational attraction of the AGB star (see Eq. (35)), as a function of distance from the AGB star, given by by MAGRITTE, for the snapshot of the simulation using the attenuation prescription. The colourbar is in log-scale.
In the text

	Fig. 17 Eddington factor Γ (Eq. (5)), calculated using the magnitude of the flux vector, in a slice through the orbital plane, using MAGRITTE for the snapshot of the simulation, using the attenuation prescription.
In the text

	Fig. A.1 Relative error (Eq. A.1) and computation time, required for calculating τ as a function of the number of rays traced in the reference model. The solid red line represents the relative error using the SPH nearest neighbours, blue the computation time, and green the number of rays in the IRS case.
In the text

	Fig. A.2 Relative error (Eq. A.1) as a function of the number of rays used in the interpolation (using HEALPIX order 5).
In the text

	Fig. A.3 Relative error (Eq. A.1) as a function of computation time for calculating τ. Different dots represent different HEALPIX orders, where the blue represents no interpolation, and the orange represents the interpolation using four rays and k = 2.
In the text

	Fig. B.1 Lucy optical depth, τL (Eq. 23), as a function of radial distance in the single-star simulation, for the high mass-loss rate (3 × 10−6 M⊙ yr−1) model.
In the text

	Fig. B.2 Optical depth, τ (Eq. 28), in the attenuation prescription as a function of radial distance, for the high mass-loss rate (3 × 10−6 M⊙ yr−1) model.
In the text

	Fig. B.3 Density distribution in a slice through the meridional plane for the four prescriptions: free-wind (top left), geometrical (top right), Lucy (bottom left), and attenuation (bottom right) prescription, for the high mass-loss rate case with a binary companion. Both stars are on the x-axis, where the primary AGB star is on the left, and the companion on the right.
In the text

	Fig. B.4 Density distribution in a slice through the orbital plane close to the companion in the Lucy approximation, revealing the presence of an accretion disk. The small arrows indicate the direction of the velocity field.
In the text