© The Authors 2025

1 Introduction

Low- to intermediate-mass stars are stars with typical masses ranging between ~0.8 and ~8 M_⊙. After the depletion of hydrogen and helium in their cores, these stars enter the asymptotic giant branch (AGB) phase. In this stage of their evolution, stars are characterised by an extended dusty envelope, and a stellar wind with typical mass-loss rates ranging between 10⁻⁸ and 10⁻⁵ M_⊙ yr⁻¹, and terminal wind velocities between 5 and 30 km s⁻¹ (Höfner & Olofsson 2018). The mechanism driving these winds is generally believed to be a combination of pulsations and radiation pressure on dust grains. Pulsations levitate material high up in the circumstellar envelope, where the cooler temperature allows the formation of dust grains that can efficiently absorb the infrared stellar radiation, thereby accelerating outwards and dragging along the surrounding gas (Lamers & Cassinelli 1999). Depending on the mass, stars in the AGB phase can also experience third dredge-up events, where the outer convective envelope penetrates deeply into the star, enriching the envelope with material originating from nuclear burning and the s-process (Merrill 1952; Schwarzschild & Härm 1965; Kippenhahn et al. 2013). All together, AGB stars contribute ~85% of all the gas and ~35% of all the dust to the interstellar medium (Tielens 2005). Galactic chemical evolution models predict that AGB stars are the main source of carbon and nitrogen, and about half of the elements heavier than iron in the Universe (Kobayashi et al. 2020).

The high-resolution observations of these outflows, obtained with instruments such as the Atacama Large Millimeter/Submillimeter Array (ALMA) interferometer, reveal a variety of complex morphologies that contrast with the longheld assumption of spherical symmetry (see e.g. Maercker et al. 2016; Kervella et al. 2016; Homan et al. 2018; Decin et al. 2020). The leading hypothesis to explain these structures is a previously unrecognised binary companion that shapes the outflow as the mass-losing AGB star and its companion orbit around the common centre of mass (CoM), as first numerically predicted by Theuns & Jorissen (1993). The post-AGB stars are also known to harbour complex shapes in their circumstellar envelopes (CSEs), including bipolar nebulae, shells, consecutive rings, and arcs. Similarly, these structures have been linked to the presence of companions (van Winckel 2003; Bollen et al. 2022). Hence, the dominant shaping mechanism of the CSE is believed to be the gravitational interaction of the AGB wind with a (sub-)stellar binary companion orbiting around the AGB star (Decin et al. 2020). This assumption is supported by population synthesis since the binarity fraction of AGB progenitors is found to be above ~50%, even reaching ~100% when planetary companions are also included in the binarity rate (Moe & Di Stefano 2017; Fulton & Petigura 2018).

Unfortunately, high angular resolution interferometric observations, where the impact of a companion is readily observable, remain scarce, and are available for only ~20 systems. More frequently, studies of AGB outflows are limited to the emission lines from molecules present in these outflows. Over 100 molecules have been detected in the outflows of AGB stars (Giesen et al. 2006; Decin 2021), which is roughly 50% of all molecules detected in outer space. The most-studied molecule is CO, whose lower energetic rotational states are easily excited over much of the CSE. By studying these lower transitions, we can probe the entire wind and through modelling, we can derive the massloss rate (see e.g. Kerschbaum & Olofsson 1999; Olofsson et al. 2002; Decin et al. 2006; De Beck et al. 2010). In this modelling procedure, the assumption of a single star is always made. Given the likely presence of companions, it is possible to make systematic errors on the derived parameters, especially the mass-loss rate (Decin et al. 2019). Hence, it is crucial that we model the CO emission lines with the inclusion of the binary companion.

The first steps towards this goal have already been undertaken in Homan et al. (2015, 2016) utilising simpler analytical models to study the effect of different morphologies. Moreover, Kim et al. (2019) made use of hydrodynamical simulations, but did not consider detailed radiative transfer to create synthetic spectral lines. Instead, they considered the column density within velocity bins as a proxy for the spectral line. In contrast, here we use 3D hydrodynamical simulations combined with full non-local thermodynamic equilibrium (NLTE) radiative transfer, thereby including more adequate structures and physics in our models. Our goal is to conduct an initial investigation of how various morphologies arising from these hydro-simulations translate to spectral lines.

Three-dimensional hydrodynamic studies of binary AGB systems have been performed using both grid-based and smoothed particle hydrodynamics (SPH) codes. The early studies of Theuns & Jorissen (1993) and later Mastrodemos & Morris (1998, 1999), already showed that Archimedes-like spiral structures, elongated morphologies, and accretion discs can form depending on the orbital and wind parameters. Since then, higher resolution simulations became possible, facilitating more comprehensive studies of the formation of different structures, including spirals, arc- and ring-shaped morphologies, equatorial density enhancements, and accretion discs, and the effects of these structures on the mass transfer and orbital evolution (see e.g. Chen et al. 2017; Saladino et al. 2018, 2019; Saladino & Pols 2019; Chen et al. 2020; El Mellah et al. 2020; Maes et al. 2021; Malfait et al. 2021; Lee et al. 2022; Aydi & Mohamed 2022; Malfait et al. 2024a,b).

This paper is outlined as follows. In Sect. 2, we describe the numerical set-up and parameter space of the HD simulations, outline the radiative transfer to obtain the spectral lines, and describe how we take into account CO photodissociation. Next, in Sect. 3, we present the hydro-simulations and the resulting spectral lines for the grid of models. In Sect. 4, we discuss the implications and accuracy of our results and outline the necessary steps for the future. We conclude in Sect. 5.

2 Methodology

2.1 Hydrodynamic simulations

We use the open-source smoothed particle hydrodynamics (SPH) code PHANTOM¹ (Price et al. 2018) to model the outflows of AGB stars. SPH codes solve the hydrodynamic equations using a Lagrangian method. These equations are given by $$\begin{array}{cc}\frac{\text{d}\rho }{\text{d}t}=-\rho \left(\mathbf{\nabla }\cdot \mathbf{\upsilon }\right),& \left(\text{continuity\hspace{0.17em}equation}\right)\end{array}$$(1) $$\begin{array}{cc}\frac{\text{d}\upsilon }{\text{d}t}=-\mathbf{\nabla }\left(\frac{P}{\rho }\right)+\mathbf{a},& \left(\text{conservation\hspace{0.17em}of\hspace{0.17em}momentum}\right)\end{array}$$(2) $$\begin{array}{cc}\frac{\text{d}u}{\text{d}t}=-\frac{P}{\rho }\left(\mathbf{\nabla }\cdot \mathbf{\upsilon }\right)+\Lambda ,& \left(\text{conservation\hspace{0.17em}of\hspace{0.17em}energy}\right)\end{array}$$(3)

where ρ is the gas density, u the velocity, P the pressure, u the internal energy, a represents all accelerations from external forces, and Λ denotes additional energy losses or gains, which in the present case includes shock heating and HI cooling.

To close this system of differential equations, we assume a polytropic equation of state (EOS), given by $$P=\left(\gamma -1\right)\rho u$$(4)

with γ the polytropic index. The gas temperature T can be obtained by using Eq. (4) and the ideal gas law $$P=\frac{\rho k_{B}T}{\mu m_{\text{H}}},$$(5)

where k_B is Boltzmann’s constant, m_H the mass of a hydrogen atom, and µ the mean molecular weight. In previous studies, the polytropic index was fixed to a value of 1.2 (Maes et al. 2021), and µ was taken to be constant at a value set to 2.38 or 1.26 to mimic either a molecular or atomic gas, respectively (Malfait et al. 2024a). To make our simulations more physically consistent, we calculate the values for γ and µ as done in Siess et al. (2022). For the calculation of the mean molecular weight, we take into account only the most abundant species (H, He, and H₂), giving $$\mu =\frac{\left(1+4{ϵ}_{\text{He}}\right)n_{\langle \text{H}\rangle }{k}_B{T}_g}{P_{\text{H}}+P_{{\text{H}}_2}+{ϵ}_{\text{He}}{n}_{\langle \text{H}\rangle }{k}_B{T}_g}$$(6) $$\gamma =\frac{5P_{\text{H}}+7P_{{\text{H}}_2}+5{ϵ}_{\text{He}}{n}_{\langle \text{H}\rangle }{k}_B{T}_g}{3P_{\text{H}}+5P_{{\text{H}}_2}+3{ϵ}_{\text{He}}{n}_{\langle \text{H}\rangle }{k}_B{T}_g}.$$(7)

In these expressions, n_〈H〉 is the total number of hydrogen atoms per unit volume in the medium, ϵ_He is the total number of helium atoms per hydrogen atom which is taken to be 1.04 x 10⁻¹. P_H and $P_{{\text{H}}_2}$ represent the partial pressures of atomic and molecular hydrogen, respectively. As for the calculation of μ, contributions from other molecules are neglected because they contribute less than 1% of the total mass fraction (Siess et al. 2022). This more physical treatment of the mean molecular weight and adiabatic index impacts the heating and cooling throughout the simulation, which is needed to generate more realistic temperatures.

Additional cooling is implemented in the form of HI cooling, as done in Malfait et al. (2024a). This is calculated using the relation $${\Lambda }_{\text{cool}\_\text{HI}}=7.3\times {10}^{-19}\text{g\hspace{0.17em}c}{\text{m}}^3\times \frac{n_e{n}_{\text{H}}{e}^{-118400\text{K}/T}}{\rho }\text{erg\hspace{0.17em}}{\text{g}}^{-1}{\text{s}}^{-1},$$(8)

where n_e and n_H are the electron and total hydrogen number densities (Spitzer 1978). This cooling becomes effective at temperatures greater than 3000 K, and is dominant at temperatures higher than 8000 K. As shown in Malfait et al. (2024a), this term is necessary to deal with the high temperatures generated in the high-pressure regions around the companion stars, leading to unphysical instabilities and energetic outflows if not included, as present in Malfait et al. (2021). This term, together with the adiabatic work due to compression and expansion, and shock heating Λ_shock (Price et al. 2018) are included in Eq. (3).

Both the primary and secondary star are modelled as sink particles. The primary and secondary have an accretion radius R_accr, beyond which particles are accreted if their energy is low enough (Price et al. 2018). The actual wind is driven in a free-wind approximation, where the gravity of the primary AGB star is artificially balanced by the radiation force (Theuns & Jorissen 1993). To drive the wind, SPH particles are distributed on spheres around the primary sink particle, and injected into the simulation with a fixed initial velocity (see Siess et al. 2022).

2.2 Grid set-up

We computed a grid of nine models with varying injection velocities and orbital separations. In previous studies (e.g. Malfait et al. 2024a), the wind was launched from the stellar surface. However, the wind is launched at the location where dust grains can start forming, i.e. the dust condensation radius. Injecting wind particles from the stellar surface gives rise to unrealistic temperature profiles in the inner region of the simulation. To estimate the wind injection radius, we assume a simple power law for the temperature starting from the surface of the star (Decin et al. 2006). This is given by $$T\left(r\right)=T_{\ast }{\left(\frac{R_{★}}{r}\right)}^{\zeta }$$(9)

with ζ ≈ 0.5 for optically thin absorption and emission, and T_★ and R_★ the temperature and radius of the star, respectively. We choose values for T_★ and R_★ of 2750 K and 1.267 au, respectively. These values are selected arbitrarily, and are in line with typical values of AGB stars (Habing & Olofsson 2004). Assuming a dust condensation temperature of 1500 K (Whittet 2003), the dust condensation radius is located at R_d ≈ 4.26 au, which is where we inject the free wind.

In the current paper, we limit ourselves to a low massloss rate of 10⁻⁷ M_⊙ yr⁻¹ for all our models. There are two main reasons for this. First of all, since we do not yet include radiative acceleration and dust formation in our simulations, increasing the mass-loss rate only increases the density. If radiation pressure on dust grains is taken into account, the larger opacities expected when using higher mass-loss rates can alter the resulting morphology. Including more adequate descriptions of radiation shows that the discrepancy between the free wind approach and detailed descriptions increase with the mass-loss rate, and is small at lower mass-loss rates. This shows that the free wind approach is most viable for these lower mass-loss rates (Esseldeurs et al. 2023). Secondly, the photodissociation radius strongly depends on the density, and in turn on the mass-loss rate. Increasing the mass-loss rate therefore also leads to larger models, and thus to a strong increase in the computation time of our models (see also Section 2.4).

In earlier studies of hydro-simulations of AGB stars, for example Malfait et al. (2024a), particles reaching a critical radius were deleted from the simulation to increase computational efficiency. This allowed the authors to use higher resolutions, whilst keeping computation time moderate. Since we are interested in low-J spectral lines, the size of the emitting envelope is of extreme importance, and we cannot use an outer boundary. In order to keep computational cost moderate, we therefore have to lower the resolution of our simulations. We have run test simulations with high resolution as well, and did not find significant differences in the inner wind morphology, as all the characteristic features of the high-resolution run are present in the low-resolution runs.

For all our models, the primary AGB star has an initial mass of M_p = 1.5 M_⊙, an effective radius of R_p,_eff = 1.267 au, an accretion radius of R_p,_accr = 1.220 au and a temperature of T_p,_eff = 2750 K. The secondary has a mass of M_s = 1 M_⊙ and an accretion radius of R_p,_accr = 0.2 au. This accretion radius is rather large for the companion, but this was necessary to remove numerical instabilities as a result of the lower resolution. A lower value for the accretion radius of the secondary, in combination with a higher resolution, is especially needed when interested in resolving the accretion disc (see Malfait et al. 2024a), at the cost of a higher computation time. Since we are looking at the simulations on a larger scale, we are not interested in resolving the accretion disc in great detail. This will not impact the resulting line profiles, as the region inside the accretion radius is negligible compared to the total region probed by the low-J CO spectral lines. Additionally, if thermal dissociation were to be taken into account, we do not expect any CO to be present in the disc. The set-up parameters are summarised in Table 1.

We consider three different orbital separations of 9, 15, and 25 au, corresponding to orbital periods of approximately 17, 37, and 72 yrs, respectively. All models are assumed to be circular (i.e. e = 0). Lastly, we take three different initial velocities v_ini for the free wind, of 5, 10, and 20 km s−1. The list of models are summarised in Table 2. It should be noted that for the higher outflow velocities, observations suggest a higher mass-loss rate (Ramstedt et al. 2020). However, to explore the parameter space these models still remain interesting to study and understand.

We run the simulations for the v05 and v10 models until ≈2560 yr, and the v20 models until ≈ 1910 yr (see Section 2.4). The number of particles in each simulation varies with the initial velocity, as the time between consecutive shells of injected particles can be written as (Siess et al. 2022) $$\delta t=w_{ss}\frac{d_{\perp }}{{\upsilon }_{\text{ini}}}$$(10)

with d_⊥ the distance between neighbouring particles on an injection shell, w_ss a value close to 1, and v_ini the initial wind velocity.

This means that increasing the initial velocity of the wind, decreases the time between consecutive shell ejections, leading to more particles. At the end, the models with initial wind velocities of 5, 10, and 20 km s−1 have approximately 2 x 10⁶,4 x 10⁶, and 6 x 10⁶ particles, respectively.

2.3 Radiative transfer

To post-process our hydrodynamic simulations, we use MAGRITTE² (De Ceuster et al. 2019, 2020, 2022; Ceulemans et al. 2024), an open-source 3D NLTE line radiative transfer code. MAGRITTE solves the radiative transfer equation along a fixed set of rays using a second-order Feautrier scheme. The ray tracing algorithm is grid-agnostic, and can easily cope with SPH models as input. For the CO abundance, we assume a constant fractional number density of CO to hydrogen of 1.5 x 10⁻⁴, i.e. $$\frac{n_{\text{CO}}}{n_{\text{H}}+2n_{{\text{H}}_2}}=1.5\times {10}^{-4},$$(11)

such that when all hydrogen is in H₂, we have a constant fractional abundance of 3 x 10⁻⁴, which is a commonly assumed value for oxygen-rich AGB stars (Teyssier et al. 2006), but this value can be as high as 8 x 10⁻⁴ for carbon-rich AGB stars (Li et al. 2014). The number densities of H and H₂ are calculated using equilibrium chemistry, such that the following condition holds: $$P_{{\text{H}}_2}=P_{\text{H}}^2{K}_{{\text{H}}_2}.$$(12)

Here $K_{{\text{H}}_2}$ is the dissociation constant of the reaction given by $$K_{{\text{H}}_2}=\exp \left(\frac{-\Delta G_{{\text{H}}_2}}{R{T}_g}\right)$$(13)

with R the universal gas constant, and ∆GH₂ the change in Gibbs energy. It follows that $$n_{\langle \text{H}\rangle }k{T}_g=P_{\text{H}}+2P_{\text{H}}^2{K}_{{\text{H}}_2}.$$(14)

Once the partial pressures are known, the number densities are found using $n_i=\frac{P_i}{k{T}_g}$. After this, CO photodissociation is applied to obtain a realistic CO profile (see Sect. 2.4). Additionally, we include a turbulent velocity field of 1 km s⁻¹ throughout the domain in our MAGRITTE simulations, which impacts the line broadening. This value is in line with previous studies, although it can be taken as high as 2 km s⁻¹ (see e.g. Ramstedt et al. 2006, 2008; De Beck et al. 2012).

In our NLTE calculations, we take into account the first 40 rotational transitions of CO in the ground vibrational state. All the molecular data are taken from the Leiden atomic and Molecular Database³ (Schöier et al. 2005, 2010). The collisional rate coefficients are adopted from Yang et al. (2010). For all v05 and v10 models, we add a spherical outer boundary inside MAGRITTE at 5000 au, and for the v20 models this is located at 6000 au. This is done in order to remove the non-physical outer regions of the SPH model where the velocities are very high due to the expansion into vacuum. These boundaries are, however, well outside the emitting CO region (see Section 2.4 and Fig. 2), and thus do not impact the line profiles. We add the cosmic microwave background radiation at the outer boundary of the model, and add a spherical inner boundary on the surface of the AGB star (i.e. on R_p,eff), which is emitting as a blackbody with a temperature of T_★. We always put the observer at a distance of 500 pc, which is a realistic distance for the observed AGB stars (Andriantsaralaza et al. 2022), changing it would uniformly scale the observed flux proportional to the ratio of the distances squared. Finally, we always use velocity bin sizes of 0.1 km s⁻¹ for the spectral lines.

The region between the injection radius and the surface of the AGB star is left empty, but can still contribute to the spectral lines, albeit towards the higher J lines. We sample this region before calculating the spectral lines. Since many uncertainties are present regarding the complex structure, and a detailed treatment is out of the scope of this paper, we use simple prescriptions for the quantities in this region, similar to Decin et al. (2006). As stated in Section 2.2, we assume a power law for the temperature according to Eq. (9). The density is given by the mass conservation equation $$\rho \left(r\right)=\frac{\dot{M}}{4\pi r^2\upsilon \left(r\right)},R_{★}<r<R_{\text{inject}}$$(15)

and the velocity is assumed to follow a classical β law of the form $$\upsilon \left(r\right)={\upsilon }_{\infty }{\left(1-\frac{R_{★}}{r}\right)}^{\beta },R_{★}<r<R_{\text{inject}}$$(16)

with β = 0.5, and v_∞ such that v(r) equals the initial wind velocity of the model at r = R_inject. In the inner region, the CO number density is calculated as explained above.

Because of the computational cost that comes with full NLTE calculations, MAGRITTE can recursively remesh the grid, thereby reducing the computational load (see Ceulemans et al. 2024, for details). This method has two parameters, the maximum variation of the density in a box r_max, and maximum recursion level l. The exact number of points that are left after reduction depends on these two parameters, as well as the original model. We opt for a value of l = 15, and increase the r_max parameter until we are left with one quarter of the original number of points. This was found to be optimal for reducing the computational cost without sacrificing accuracy. Additionally, we benchmarked our results to the Monte-Carlo radiative transfer code SKIRT (Baes et al. 2003), for which the results can be found in Appendix A.

2.4 CO Photodissociation

When modelling low-J CO transitions, the size of the envelope is of high importance, as it greatly influences the observed flux coming from the lowest transitions. This size is set by external ultraviolet (UV) radiation, that destroys CO from outside inwards, determining the CO distribution. When modelling CO spectral lines using spherically symmetric models, one can use the tabulated values of CO density profiles, as calculated by e.g. Mamon et al. (1988), and more recently Groenewegen (2017) and Saberi et al. (2019).

When dealing with more complex morphologies, the spherically symmetric description is insufficient, since the shielding from UV radiation by CO itself, and by H₂, varies depending on the direction of the incoming UV photons due to the flattening of the model and local density variations in a spiral structure. Therefore, to calculate the CO density profiles, we implement photodissociation in 3D, which closely follows the method described by Groenewegen (2017), which is recapitulated and generalised below.

Assuming no molecules are created, and UV photons are the only destruction mechanism, the CO number density, n, along a ray can be written as (Jura & Morris 1981) $$\frac{1}{r^2}\frac{\partial }{\partial r}\left(r^{2}n\upsilon \left(r\right)\right)=-I\left(r\right)n,$$(17)

where v(r) is the expansion velocity, and I(r) is the photodissociation rate at a distance r to the central AGB star. In order to find a solution, we define x as r²nv(r), leading to $$x\left(r_i\right)=x\left(r_{i-1}\right)\exp \left(-{\displaystyle {\int }_{r_{i-1}}^{r_i}\frac{I\left(r^{\prime }\right)}{\upsilon \left(r^{\prime }\right)}d{r}^{\prime }}\right),$$(18)

where one assumes x(r₀) = 1. The dissociation rate I(r) is given by $$I\left(r\right)=\frac{1}{4\pi }{\displaystyle \int k\left(r,\theta ,\varphi \right)d\Omega ,}$$(19)

where k(r, θ, ϕ) is the dissociation rate at a distance r from interstellar photons that are coming at a direction defined by θ and ϕ towards point r. In the spherically symmetric case, the integral only needs to run over θ. The dissociation rate k(r, θ, ϕ) can be written as $$k\left(r,\theta ,\varphi \right)=\chi I_0{\Theta }_{dust}\left(r,\theta ,\varphi \right){\Theta }_{{\text{H}}_2,CO}\left(r,\theta ,\varphi \right),$$(20)

where I₀ = 2.6 x 10⁻¹⁰ s⁻¹ is the unshielded dissociation rate (Visser et al. 2009), χ indicates the strength of the interstellar radiation field (ISRF) relative to the Draine field (Draine 1978) used in Visser et al. (2009), Θ_dust(r, θ, ϕ) represents the amount of shielding by dust, and ${\Theta }_{{\text{H2}}_{}}{}_{\text{,CO}}\text{(}r,\theta ,\varphi \text{)}$ corresponds to the CO self shielding and H₂ shielding. For the latter, we use the tabulated shielding functions of Visser et al. (2009), which depend on the CO and H₂ column densities, and the excitation temperature T_ex, and are interpolated accordingly. We assume the excitation temperature equals the kinetic temperature. The shielding by dust is given by Θ_dust = exp(-τ_dust), with τ_dust the dust optical depth. In our PHANTOM models, we do not yet include dust; therefore, we opted to consider a constant dust extinction at 1000 Å, given by $${\tau }_{\text{dust}}\left(r,1000\AA \right)=\frac{4.65\times 2\times N_{{\text{H}}_2}\left(r\right)}{1.87\times {10}^{21}}$$(21)

with $N_{{\text{H}}_2}$ (r) the H₂ column density (Morris & Jura 1983). Finally, we assume χ = 1.

Since the dissociation rate depends on the CO column density, which in turn depends on the CO density profile, we need to solve this iteratively. Since we now work in 3D, this would entail tracing a ray through every point, and for each point along that ray, calculate the column densities along all directions, and finally solve Eq. (18) along all the rays. As this is computationally demanding, we instead opt for an approximation, where we apply the process described above for a predetermined amount of rays, which are subsequently interpolated to obtain values of x(r) for the full domain, similar to Esseldeurs et al. (2023). The full implementation is outlined in Appendix B.

One also has to keep in mind the initial size of the model, as the H₂ dissociation is not being accounted for, and thus, the H₂ and dust shielding is always increasing. In Groenewegen (2017), the size of the model is taken to be 15$R_{\frac{1}{2}}$ (i.e. 15 times the radius where the fractional abundance of CO is half its central value), but running our hydro-simulations to such a size is computationally unfeasible. Instead, we run the model sufficiently large, such that the value of $R_{\frac{1}{2}}$ does not change by more than one percent. We find that running the v05 and v10 models to 2560 years suffices for this, whilst for the v20 models, this state is reached sooner, and we stop the simulations at 1910 years. Running the model even further could still lead to slightly different profiles in the outer regions, but this should only marginally impact the J = 1→0 line, and this is prone to larger uncertainties such as the exact value of χ.

3 Results

3.1 Morphologies

As a novel part of the wind parameter space is explored, the different resulting morphologies are briefly explained. Still, the general structures that are observed are similar to those found in earlier studies. Thus, we omit a detailed description of the formation of these structures, and for a far more in-depth explanation of how these specific morphologies are created, we refer to Maes et al. (2021) and Malfait et al. (2021, 2024a) and references therein. We focus here on the inner region of the models first, and discuss the overall global morphology (after we have applied CO photodissociation) in the subsequent section.

3.1.1 Inner wind morphology

A companion influences the outflow in two main ways. First, it interacts with the outflowing wind particles due to gravitational attraction. Secondly, it induces an orbital motion of the AGB star around the CoM. We show the resulting mass density for the v10 models in slices through both the orbital (x − y) and meridional (x − z) plane for the inner regions in Fig. 1, with increasing orbital separation from left to right. Analogous plots for the v05 and v20 models can be found in Figs. C.1 and C.2.

For the initial velocity of 10 km s⁻¹ (Fig. 1) a clear spiral structures arises in the orbital plane, regardless of the separation. In all three cases, the spiral consists of a more dense inner edge, the ‘backward spiral edge’ (BSE), and a less dense outer edge, the ‘frontal spiral edge’ (FSE) as defined in Malfait et al. (2021). This FSE has a greater radial velocity, and therefore the spiral wake will widen, eventually catching up with the BSE of the spiral sent out earlier, and form a single spiral shock. This interaction is more pronounced at the smaller separations.

Although mostly similar in the orbital plane, the morphologies in the meridional plane are very distinct. Looking at the largest separation of 25 au (right column), the density distribution consists of concentric arcs, where the inner part of each arc is the BSE, and the outer edge is the FSE. At the ends of the arcs, there are thinner regions, resulting from particles that are not gravitationally attracted by the companion. At a separation of 15 au (middle column), we still see the same concentric shells, but the catching up of the FSE and BSE occur sooner, removing also the thin ends of the arcs. At the smallest separation of 9 au (left column), the arcs have a very distinct shape, being elongated in the orbital plane, indicating the spiral does not reach very high above the orbital plane.

Compared to the lower initial outflow velocity of 5 km s⁻¹ (Fig. C.1), the spiral structures look very similar to the v10a09 model, where the FSE quickly catches up with the BSE. One notable difference is the appearance of a bow-shock spiral for the v05a15 and v05a25 models. This can be seen by the presence of a second spiral structure arising in front of the companion in the outer region of the regular spiral. Interesting to note is the absence of this feature in the v05a09 model, as it is expected to be more present when the wind velocity at the location of the companion is lowest. In the meridional plane, we see structures that are reminiscent of those observed for the v10a09 model, that are elongated along the orbital axis, although the exact shape differs between the models.

For the highest initial outflow velocity of 20 km s⁻¹ (Fig. C.2) we see a very well-defined narrow spiral structure in the orbital plane for all three models. Still there is a FSE and BSE present; however, only at the smallest separation can we observe the widening of the spiral and the FSE catching up with the BSE at this scale. In the meridional plane, the concentric shells reach nearly up to the polar axis, and are not too dissimilar from those observed in the v10a25 model, although the shells appear more circular, whilst those in the v10a25 model are more triangular, because of a stronger interaction of the companion with the wind.

Fig. 1 Density distributions in slices through the orbital plane (upper row) and the meridional plane (lower row) for the models with initial wind velocities of 10 km s−1. From left to right, the orbital separation is 9, 15, and 25 au.

3.1.2 Large scale morphology and CO photodissociation

The global morphologies of all nine models, after the calculation of the CO number density, are displayed in Fig. 2, where the solid and dotted lines represent the location where fractional abundance of CO with respect to H₂, compared to the initial fractional abundance, is 0.5 and 0.01, respectively. We always show a slice through the x − z plane, as in the orbital plane the CO map is circular in all models. In the meridional plane, some interesting features are visible. Due to the motion of the AGB star about the CoM, the wind particles are injected into the simulation with a velocity $${\mathit{\upsilon }}_{\text{wind}}={\mathit{\upsilon }}_{\text{ini}}+{\mathit{\upsilon }}_{\text{orb},\text{AGB}}$$(22)

with v_ini the injection velocity, and v_orb,AGB the orbital velocity of the AGB star around to CoM. As noted in Malfait et al. (2021), the wind particles ejected perpendicular to the orbital plane follow a radial direction with an angle $$\theta ={\tan }^{-1}\left(\frac{{\mathit{\upsilon }}_{\text{ini}}}{{\mathit{\upsilon }}_{\text{orb},\text{AGB}}}\right)$$(23)

with respect to the orbital plane. This deflection of the particles’ trajectory causes the global morphology to be flattened.

This is readily observed in Fig. 2, where the v05 models all show very flattened morphologies, whilst the v20 models are more spherical. For the v10 models, flattening is strong for the v10a09 model, only marginal for the v10a15 model, and very small for the v10a25 model, which is in line with expectations as the orbital velocity around the CoM is lowered for larger separations. In the meridional plane, the arcs also reach lower latitudes and are elongated along the orbital axis with decreasing separation. The gravitational attraction of the companion can further focus particles to the orbital plane, building up pressure gradients that can cause particles to be accelerated towards lower density regions. At larger scales, we also see that the FSE still catches up with the BSE for the larger separations for the v20 models.

Due to UV photodissociation, the fractional CO abundance in the orbital plane, with respect to the central fractional abundance, drops to half for all models between 1500 and 2000 au, and for the flattened morphologies, this is already at ~1000 au in the z direction. This inner part is important for most CO transitions as they are only excited within this region, and the outer regions only impact the J = 1→0 and marginally impact the J = 2→1 lines. The CO fractional abundance, relative to the central value, drops to 1 percent at ~4000 au for the v05 and v10 models in the orbital plane, and 2000 au in the z direction. For the v20 models this occurs around 5000 au.

Overall, the global CO distribution closely follows the inner wind morphology, as the shape of the R_½ contour (or where the fractional abundance of CO is half of its initial value) closely resembles the shape of the arcs in the meridional plane, which is best visible for the v05 models. Even for the v20 models, where no flattening occurs, the contours show deviations from spherical symmetry with small bumps towards the polar regions. This also shows the importance of accounting for the CO photodissociation in a detailed manner. Finally, we note that when looking at the global morphology, the interaction of the binary companion with the wind (besides the global flattening), is hardly visible for the smallest separations, and becomes more visible with increased separations.

Fig. 2 Number density of CO in slices through the meridional plane after accounting for photodissociation. The solid and dotted lines represent the location where x(r) is 0.5 and 0.01, respectively.

Fig. 3 Synthetic CO J = 2→1 spectral lines for all models. The columns correspond, from left to right, to initial outflow velocities of 5, 10, and 20 km s−1. The solid, dashed, and dotted lines represent orbital separations of 9, 15, and 25 au. The spectral lines are creating viewing face-on at an inclination of 0° (upper row) and edge-on at 90° (lower row).

3.2 Synthetic spectral lines

For the CO spectral lines, not only the density structures, but also the velocity distribution is important, as the latter dictates the frequency at which certain material will be emitting. Additionally, the temperature profile is of great relevance, since this greatly impacts the level populations. Given the asymmetry of our simulations, we create spectral lines at six different inclinations between 0° and 90°, where 0° is defined as looking straight at the orbital plane, or face-on, whilst for 90° we are looking at the meridional plane, or edge-on. When viewing at an angle, the position angle (PA) of the mass-losing AGB star also can play a role, as it changes the appearance of the spiral structure. We define a PA of 0° when both the AGB star and the companion are on the x-axis, with the AGB star on the left (i.e., as they are presented in Fig. 1), and anti-clockwise is positive. Our focus here lies on the low-J lines, starting from the J = 1→0 up to the J = 6→5 transition, as these lines are often used to derive mass-loss rates from observations. In the following sections, any telescope effects are omitted. These will be discussed in Sects. 3.2.5 and 4.2.

3.2.1 Face-on and edge-on line profiles

Fig. 3 shows the CO J = 2→1 spectral lines for all nine models viewed both face and edge-on, always at a PA of 0°. As noted before, this line is excited throughout the wind, and thus should trace the overall morphology. For the v05 models (left column) very defined double peaks emerge at the maximal radial velocities when viewed face-on, and a strong central bump appears when viewed edge-on, regardless of the orbital separation. These characteristics are are still present for the v10a09 and v10a15 models (middle column), with the central bump for the v10a09 model being very prominent. The two peaks still show up for the v10a25 model when viewing face-on, but the central bump has completely disappeared for this model, and the line profile approaches a simpler parabolic shape. Looking at the v20 models (right column), the profiles show less pronounced features. Some peaks are still present when viewed face-on, where, except for the v20a09 model, the peaks are no longer present at the outer region of the spectral line, but have moved closer to the central velocities. Viewing edge-on, the line profiles are approaching flat topped line profiles, characteristic of optically thin spherically symmetric outflows (see also Fig. 5).

Additionally, we see an increase in the width of the spectral lines, when going from the face-on to the edge-on view. For the v05 models, the line width is increased by approximately a factor of two for the three orbital separations. For the v10 models, it changes between the models, with again a factor of two for the v10a09 model, 1.5 for the v10a15 model, and for the v10a25 the width has remained similar between the two inclinations.

When comparing the different spectral lines more closely, we see that the flux changes significantly between the different models, especially when comparing different outflow velocities. This can be attributed to several reasons. First of all, increasing the wind velocity naturally widens the line profiles, but also decreases the density given a constant mass-loss rate (Eq. (15)), thereby lowering the amount of flux received per velocity bin. Secondly, the interaction with the companion also drastically alters the temperature of the wind, which in turn greatly affects the level populations, and therefore the received flux. As described in Section 3.1, this modification strongly depends on both the local wind velocity as well as the orbital separation.

To gain insight into the shapes of these spectral lines, we show the velocity distribution of the hydro simulations in the z and x direction (i.e., v_z and v_x) in slices through the meridional plane. For the v10 models, we show this in Fig. 4. Analogous plots for v05 and v20 are presented in Figs. C.3 and C.4, respectively. In these plots, the observer is located above the plot in the upper row, and to the right for the lower row. Additionally, two sets of contours are indicated to visualise velocity bins corresponding to interesting features of the spectral lines. For reference, in a spherically symmetric wind with a constant outflow velocity, the volume of equally sized velocity bins is equal and are conically shaped, leading to flat topped and parabolic line profiles in optically thin and thick regimes, respectively. Spectral lines of the first six transitions from a spherically symmetric PHANTOM model, with a mass-loss rate of 10⁻⁷ M_⊙ yr⁻¹ and initial wind velocity of 10 km s⁻¹ are shown in Fig. 5. Note that it is also possible to obtain two peaked profiles from single star systems, depending on whether the outflow is resolved or not (see also Sect. 3.2.5).

Instead of the conically shaped velocity bins, the interaction with the companion for the v10a09 model causes material in the velocity bin around 2.2 km s⁻¹ to deviate strongly from the conical shape (as visible in the upper left plot of Fig. 4). In this case the velocity bin has an arc-like shape. The velocity bin centred around 0 km s⁻¹ does remain conical, and thus the relative contribution of the velocity bin centred around 2.2 km s⁻¹ will be greater, leading to an increase in the flux at ±2.2 km s⁻¹, and two peaks in the spectral line (Fig. 3 top middle plot, solid line). As seen in Fig. C.3, this also produces the two peaks for all the v05 models, although the exact shape of the velocity bin structures strongly differs between the orbital separations. All the double peaks in Fig. 3 also show a stronger peak at redshifted velocities, in contrast with the blue shifted peak, which is due to the opacity being generated in slightly hotter regions for the red peak (Morris et al. 1985). As we are looking face-on, this cannot be an effect of the position angle.

The central bump when viewing edge-on has a similar explanation as the appearance of two peaks. Although in this case all velocity bins are roughly conical, the volume of the bins towards the central velocity is larger compared to that of higher velocities. This is more outspoken for the v05 models (see Fig. C.3). For these models, the velocity bins deviate more from the conical shape as well, as the contours show more curved edges, due to a stronger interaction with the companion. Note also the different scales of the colour bars in the top and bottom row. In the orbital plane, higher speeds are reached compared to the meridional plane, giving rise to the increase in the width of the spectral lines described above.

The two peaks that are present in the v10a15 model are not due to an arc-like structure, but rather due to two separate regions having the same projected velocity in the z direction, thereby creating a peak in the line profile. The double peaks from the v10a25 model are yet again different in origin, this time they are created akin to the central bump feature. The conical regions are increased in size with respect to the central velocity bin (top right plot of Fig. 4), creating two distinct peaks. We also note that the individual contours show some coves, indicating that the features of the spiral are influencing the line shapes stronger than for the v05 models, where the contours are more smooth (see Fig. C.3).

Finally, the cause of the deviations from flat topped or parabolic line profiles for the v20 models are similar in origin, but as these deviations are less pronounced, they are also less visible in Fig. C.4. For the v20 models, however, the spiral features in the contours are very pronounced, which can impact the line profiles (see also Sects. 3.2.3 and 3.2.4).

Fig. 4 Velocities in the z direction (upper row) and in the x direction (lower row) in slices through the meridional plane for the models with initial wind velocities of 10 km s−1. Contours are shown on the plots around velocities interesting for the line profiles, indicated by the horizontal lines on the colour bar. From left to right, the orbital separation is 9, 15, and 25 au.

Fig. 5 Synthetic CO spectral lines of the first six transitions from a single-star PHANTOM model, with a mass-loss rate of 10−7 M⊙ yr and an initial wind velocity of 10 km s−1.

3.2.2 Intermediate inclinations

In Fig. 6 we show the CO J = 2→1 spectral lines for all v10 models viewed at intermediate inclinations of 18°, 36°, 54°, and 72°, still at a PA of 0°. Analogous plots for the v05 and v20 models are presented in Figs. D.1 and D.2, respectively. In general, there is a smooth transition between the face-on and edge-on cases, as at ever-increasing inclinations, the two peaks become less pronounced, and the central bump starts appearing. At the two most central inclinations of 36° and 54°, both features can be distinguished. This is most visible at 36° for the v10a09 model and at 54° for the v10a15 model, where two small peaks are present on a larger central emission. The spectral lines at the smallest and largest inclinations of 18° and 72° are virtually indistinguishable from the face and edge-on spectral lines presented in Fig. 3. The same findings apply to the v05 and v20 models.

3.2.3 Different CO transitions

The CO J = 2→1 line, which we have so far solely focussed on, probes mainly the outer part of the CSE. Higher transitions probe hotter regions due to a higher excitation temperature. In Fig. 7 we show how the spectral lines evolve as the transition increases from J = 1→0 up to J = 6→5 for the v05 models viewed face-on, and in Fig. 8 for the v20 models viewed edge-on. For the remaining models, we show only the J = 1→0 and the J = 6→5 transitions viewed both face-on and edge-on in Figs. 9 and 10, as the results are similar to the cases described here.

Starting with the v05 models viewed face-on in Fig. 7, we see that the peaks remain strongly visible for the models with orbital separations of 9 and 15 au, although the central dip diminishes in strength for the 15 au case towards higher transitions. For the orbital separation of 25 au, at the highest considered transition the two peaks are replaced by a single off-centre peak, as the dip has disappeared completely. Remarkable is also that at the J = 1→0 transition, the strength of both peaks is equal, due to a lower optical depth, and the effect described in Sect. 3.2.1 does not take place. In general, for the v05 models, the line shape remains largely similar when viewing face-on.

For the v20 models viewed edge-on (Fig. 8), we see a distinct behaviour. For the lower transitions (J = 1→0 and 2→1), the line shapes are quasi symmetric. However, towards higher transitions, these lines become asymmetric, with a distinct double peak forming for the v20a25 model. These peaks differ in strength, similar to the face-on line profiles of the v05 and v10 models, but at the J = 3→2 transition the red peak dominates, whilst the blue peak is stronger at the J = 6→5 transition. Hence, this feature cannot be explained by the effect described in Sect. 3.2.1, but rather by the fact that the higher J lines probe different regions for different models. To highlight this effect, we show the temperature distribution in a slice through the orbital and meridional plane in Fig. 11, for the v05a09 (left) and v20a25 (right) models, to show two extreme cases. In the v05a09 model, the temperature distribution is quite smooth, with the spiral only barely visible in the innermost region. Higher J lines probe hotter regions, but as the temperature profile is quite uniform, the overall line profile is mainly determined by the global wind morphology. In contrast, the temperature distribution of the v20a25 model is strongly dominated by the spiral structure. This implies that, while the lower J transitions are set by the global morphology, the higher J transitions probe the spiral structure where the temperature is higher. This can produce the strong asymmetries visible in Fig. 8, as the spiral structure at blue and red shifted velocities differ by half a period.

In Fig. 9 we show the J = 1→0 and 6→5 transitions for the v10 and v20 models, viewed face-on. For the v10 models, the same general results described above apply, as the overall line shape remains largely the same, with some differences in the dip between the two peaks. For the v20 models, the line shapes change significantly, due to the transition from the global wind to the spiral structure. For the J = 1→0 transition, the lines are similar to those in Fig. 3, with a more pronounced double peak for the v20a09 model. At the J = 6→5 transition, all lines show two peaks, and most notably, the v20a25 model now has a central bump, not too dissimilar from the edge-on line profiles of the v05 and v10 models. This central bump comes from the fact that the concentric arcs are slightly thicker towards the orbital plane, than towards the polar regions, as visible in the lower right plot of Fig. C.2. In the face-on view the spiral structure is symmetric at the blue and red shifted velocities, and as a result, the spectral lines are also quasi symmetric, in contrast with the edge-on view.

Finally, in Fig. 10 we show the J = 1→0 and 6→5 transitions for the v05 and v10 models viewed edge-on. Here, the main features of the line, i.e. the strong central bump, are still present, but the general shape has still changed significantly with J. For both v05a09 and v10a09 at the J = 1→0 transition, two small peaks at the outer edges of the line become visible. At the J = 6→5, the central bump has increased in strength relative to the outer wings for all models. The v10a25 model, which does not possess a central bump, shows a very small dip at the centre.

Fig. 6 Synthetic CO J = 2→1 spectral lines for the models with initial wind velocities of 10 km s−1, viewed under inclinations of 18°, 36°, 54°, and 72°. The solid, dashed, and dotted lines represent orbital separations of 9, 15, and 25 au.

Fig. 7 Synthetic CO spectral lines of the first six rotational transitions, for the models with initial outflow velocities of 5 km s−1, viewed under an inclination of 0°. The solid, dashed, and dotted lines represent orbital separations of 9, 15, and 25 au.

Fig. 8 Synthetic CO spectral lines of the first six rotational transitions, for the models with initial outflow velocities of 20 km s−1, viewed under an inclination of 90°. The solid, dashed, and dotted lines represent orbital separations of 9, 15, and 25 au.

Fig. 9 Synthetic CO J = 1→0 (top) and J = 6→5 (bottom) spectral lines for the v10 (left) and v20 (right) models, viewed under an inclination of 0°. The solid, dashed, and dotted lines represent orbital separations of 9, 15, and 25 au.

Fig. 10 Synthetic CO J = 1→0 (top) and J = 6→5 (bottom) spectral lines for the v05 (left) and v10 (right) models viewed under an inclination of 90°. The solid, dashed, and dotted lines represent orbital separations of 9, 15, and 25 au.

Fig. 11 Temperature distributions in slices through the orbital plane (upper row) and the meridional plane (lower row) for the v05a09 (left) and v20a25 (right) models.

3.2.4 Effect of PA

So far, all the spectral lines have been created at a PA of 0 degrees. Given the results of Sect. 3.2.3, we can make a prediction whether the PA will play a role in the line shape as well. As long as the shapes of the line profiles are determined by the global wind, as is the case for the v05 and part of the v10 models, the PA will have a negligible effect. In contrast, as described above the higher J lines for the v20a25 model directly probe the spiral structure. Hence, by changing the PA, the projected velocities of the different parts of the spiral are altered, thereby also changing the spectral line.

This is illustrated in Fig. 12, where the J = 6→5 line viewed edge-on is shown, for four different values of the PA, differing by a quarter of a period for the v05a09 (left) and v20a25 (right) models. Changing the PA for the v05a09 has a very small effect on the line profiles, contrary to the v20a25 model where the location and the strength of the stronger peak strongly depend on the PA. This effect is strongest when viewing perfectly edge-on, and decreases in strength towards smaller inclinations, disappearing when viewing face-on. This phase dependency of the spectral line could be a diagnostic for locating the companion. Due to the self-similarity of the models, and low optical depth, spectral lines that are viewed at a PA that differ by 180° are near perfectly mirrored. This would likely no longer be the case if we introduce eccentricity in our models.

Fig. 12 Synthetic CO J = 6→5 spectral lines for the v05a09 (left) and v20a25 (right) models, viewed edge-on and at four different values of the PA (ϕ). The different colours correspond to the different values of the PA.

3.2.5 Single dish telescopes

When observing AGB outflows with single dish telescopes, the spectral features can be heavily influenced by the beam profile, depending on whether the source is resolved or not. To see how the spectral lines are influenced by the beam, we apply Gaussian beams with different widths on the CO J = 2→1 line of the three v10 models. Since the model has an outer boundary at 5000 au, and is at a distance of 500 pc, the model spans 20". We apply beams, centred on the CoM, with five different beam sizes of θ_b ∈ [20″, 5″, 3″, 2″, 1.5″, 1″]. These beam sizes are selected purely for visual purposes, and represent different ratios between the beam size and size of the model. This is done for the v10 models, both on the face and edge-on line profiles.

The resulting line profiles (Fig. 13) show that when viewing models face-on (top row), the overall line features are preserved, with two distinct peaks present regardless of beam size. This is to be expected given the origin of these peaks. Applying a Gaussian beam generally affects the flux at the central velocities more, as these have a greater extent on the sky. The peaks originate above the CoM (see the top row of Fig. 3) at the greatest velocities towards and away from us. Thus, they subtend a smaller area on the sky, which is near the centre of the beam, and are less affected. Viewing edge-on, the line shapes change quite drastically when varying the beam size. The central peak remains visible for the v10a09 and v10a15 model, although the strength compared to the wings decreases, and at the smallest beam sizes, two small peaks start to form around the central peak for the v10a09 model. This is in line with expectations, as the beam removes more flux from the central velocities. When the beam is very small, this effect can create a double peaked line profile even in spherically symmetric outflows.

Finally, for the v10a25 model, we observe a similar behaviour, where the two peaks in the face-on view remain very pronounced, even at the smallest considered beam. In the edge-on view, where the line profile starts out as quasi parabolic, the line slowly becomes flat topped, and two peaks start forming at the smallest beam. This is very similar to a spherically symmetric outflow, as only the region in front and behind the star is captured within the beam, and gas moving at larger angles with respect to the line of sight are excluded more.

Fig. 13 Synthetic CO J = 2→1 spectral lines for the models with initial wind velocities of 10 km s−1, where a beam with five different values of θb has been applied. The spectral lines are viewed face-on at an inclination of 0° and edge-on at 90° in the upper and lower row, respectively. From left to right, the orbital separation is 9, 15, and 25 au.

4 Discussion and future work

4.1 Implications

Given the strong deviations from flat-topped and parabolic line profiles, direct observations of line profiles as shown in Sect. 3 can be indirect evidence of a hidden binary companion. However, the line shape can be very degenerate with the orbital and wind parameters, as attested by the strong similarities between the v05 models and the v10a09 and v10a15 models. Additionally, all the v20 models show strong similarities for all spectral lines, making it very hard to retrieve the orbital parameters solely using a spectral line. Not only can the line shape change with the orbital parameters, as clearly shown when varying the separation of the v10 models (Fig. 3, middle column), but also its strength. This is most visible for the v05 models, whose line shape is similar between the different orbital separations, but the line strength of the J = 2→1 transition of the v05a25 model is almost twice as high as those of the v05a09 and v05a15 models (see Fig. 3). In addition, the behaviour of the spectral lines when moving to higher J transitions can also differ between models, both in line shape and strength. Hence, observing multiple lines simultaneously can alleviate this degeneracy.

Closely related to this, is the terminal wind velocity one retrieves from the lines, whose value strongly depends on both the orbital parameters, as well as the inclination. For the v05 models viewed face-on, the terminal velocity is approximately the input velocity, though it has slightly decreased for the v05a25 model. When looking edge-on, however, all three considered orbital separations have terminal velocities around 10 km s⁻¹. As we do not include radiative acceleration in our models, the acceleration of the material can only be achieved by the interaction with the secondary star. In contrast, all v20 models have terminal velocities close to the input velocity, regardless of inclination and orbital separation, as the interaction with the companion is less strong. For the v10 models, we see a strong dependence on the orbital separation in the face-on view, as the 9, 15, and 25 au models have terminal velocities of approximately 5, 7.5, and 10 km s−1, respectively. For the v10a09 and v10a15 models, this is even lower than the initial velocity. In the edge-on view, the difference between the models is smaller, as the v10a15 and v10a25 models have very similar velocities around 12 km s−1, and the v10a09 model has 10 km s⁻¹.

De Beck et al. (2010) derived a fitting function (their Eq. (9)) for the mass-loss rate, that depends on the integrated line intensity of low-J lines, the terminal wind velocity v_∞, fractional CO abundance, and parameters linked to the spherically symmetric models. Using this expression with typical values relevant for spherically symmetric models, we calculate the mass-loss rate of our models, viewed both face and edge-on, by combining the result of the analytical expression for the J = 1→0,2→1,3→2,4→3,6→5, and 7→6 transitions. It should be noted that this equation is the result of a minimisation over a grid of spherically symmetric models, where detailed cooling and heating mechanisms have been taken into account. In our hydro models, we have only included HI cooling, which strongly affects the integrated line strengths. Therefore, the exact value retrieved from this formula is not of greatest importance, and we instead wish to see the difference in retrieved mass-loss rate between the models and between inclinations. The input massloss rate for all models is 10⁻⁷ M_⊙ yr⁻¹. For reference, we ran a single star simulation in PHANTOM using an initial wind velocity of 10 km s−1 and a mass-loss rate of 10⁻⁷ M_⊙ yr⁻¹. Using the formula of De Beck et al. (2010) we obtained a mass-loss rate of ~5 x 10⁻⁸ M_⊙ yr⁻¹, which can be attributed to omitting the cooling and heating mechanisms.

The derived values for the mass-loss rate of each model, which represent the mean value and standard deviation of all the considered lines, viewed face-on and edge-on are shown in Table 3. First of all, the derived values from the binary star simulations deviate from the single star PHANTOM model, which again can be contributed to a change in temperature profile, in this case caused by heating due to compression by the binary companion. Therefore, the good agreement between input and output mass-loss rates is likely a coincidence. More interestingly, the estimates of the mass-loss rate when viewing face or edge-on can be quite different. For the v05 and v10 models, this difference between the inclinations is typically a factor of two. For the v20 models, the derived mass-loss rates agree within 1σ between the face and edge-on views. This difference can be attributed to a change in the estimate for the terminal wind velocity v_∞. When using spherically symmetric models, v_∞ is strongly connected to the density, and hence to the mass-loss rate. However, when viewing the models at different angles, v_∞ changes as described above, causing a difference in the derived mass-loss rates. The difference between the retrieved mass-loss rates of the single star and binary star simulations shows that even in the case where the line profiles approach a simple flat topped or parabolic shape (i.e. the v10a25 and v20 models), the retrieved mass-loss rate can still differ by at least a factor of two, due to impact the binary companion has on the underlying temperature profile.

4.2 Observability of line shapes

When observing CO spectral lines, one also has to worry about noise present on the observations. Additionally, we have worked at a spectral resolution of 0.1 km s⁻¹ so far; however, this can be significantly larger for observations, although single dish telescopes can reach high spectral resolutions. One can then wonder which of the features described above remain observable. A detailed description of these effects is out of the scope of this paper, and instead we increase the velocity bin size to 0.5 km s⁻¹, and add white noise on top with a 1σ uncertainty of 10% of the peak flux, corresponding to a signal-to-noise (S/N) of 10. This is shown in Fig. 14, where we show the CO J = 2→1,4→3 and 6→5 lines for the v05a09, v10a25 and v20a25 models, respectively, both viewed face-on and edge-on. The underlying real spectral line is always shown as a dashed black line.

First, looking at the v05a09 model, generally the two distinct features of the spectral lines remain distinguishable. Especially the double peaks in the face-on view remain well visible. The central bump feature in the edge-on view is also visible, although the black dotted line helps significantly in identifying it, and it can be easily misinterpreted as a triangular line profile, as the transition from the outer region of the spectral line, to the central bump is easily hidden by the noise. The double peaks of the v10a25 model also remain visible, though they have become less pronounced, and in some case could disappear. The small features on the quasi parabolic profile in the edge on view have been entirely encapsulated within the noise. Similarly, all the small features that are visible on the lower J lines of the v20 models are also likely to be hidden within the noise. At the highest transitions, some features become visible for the v20a25 model. In the face-on view, one can still distinguish the central bump. The two peaks present on this central bump are still present, although without the knowledge of the underlying spectral line, this goes unnoticed. The double peaks, whose strength depends on the PA, remain visible to a certain extent, though it is unlikely that any changes in the shape due to the PA will remain significant for our models. Fig. D.3 shows the same spectral lines, but now assuming a S/N of 20. Immediately, all the characteristic features become more pronounced. Most noteworthy is that for this S/N, the asymmetric peak for the J = 6→5 transition of the v20a25 model, likely will remain observable.

On top of the noise, we need to consider the beam profile as done in Sect. 3.2.5. In the case where the emitting region is larger than the beam, the change in line shape combined with the noise will make it very difficult to retrieve the true underlying line shape, and hence infer any information about the orbital parameters. This is especially true for any double peaked patterns, as these can originate from a spherical outflow if the outflow is resolved, and thus can be easily miscategorised. This degeneracy could be lifted by including additional observables (e.g. channel maps), but this is out of the scope of the current study, and barely available for a large sample of stars.

4.3 Comparison with earlier studies

Given the difference betweens the methods used in similar studies, it is difficult to directly compare the results. Still, we can see how our overarching results fit in. Although all our models possess a clear spiral structure, our results deviate strongly from the analytical spiral models presented in Homan et al. (2015). In this paper, the authors distinguish between narrow and shell spirals (i.e. small and large opening angles), which are embedded in a spherical outflow. Therefore, any flattening effects, as present in the v05 and v10 models, are omitted. Furthermore, they impose a temperature profile on the spiral, and thus neglect the heating by compression caused by the companion. For the narrow spiral models, they obtain spectral lines with a central bump feature, and double peaked line profiles, similar to the spectral lines presented here, although both features occur at different inclinations compared to ours, with the bump feature occurring in the face-on view, and the double peaks in the edge-on view. The origin of these line shapes is very distinct, easily explaining this discrepancy. For the shell spiral, the deviation of their spectral lines from a simpler flat topped or parabolic line profile is less pronounced but still present for some of the models. Our v20 models are representative of the shell spirals, and the lines show a similar deviation from the spherically symmetric case.

Kim et al. (2019) study four HD simulations which cover a distinct part of the parameter space compared to our models. The models have an orbital separation of 68 au, two values for the eccentricity are considered (0 and 0.8), and the gravitational influence of the companion for two of their models is turned off to isolate the effect of the motion of the AGB star around the CoM. Additionally, the spectral lines are created by calculating the column density within velocity bins, thereby omitting any optical depth or temperature dependence effect. In general, the authors find double peaked profiles when viewing face-on, not too dissimilar from the low-J lines of the v20 models. These are, however, also strongly present when viewing edge-on. We do recover double peaks in the edge-on view for the v20 models when looking at higher J lines (Fig. 8), but the lines of Kim et al. (2019) are representative of low-J optically thin lines given method they use, and so do not correspond to any of our results. A more in depth comparison will not be instructive, given the difference in method and orbital parameters.

Fig. 14 Synthetic CO spectral lines for three selected models. From left to right these are v05a09, v10a25, and v20a25 with the J =2→1, 4→3, 6→5 transitions, respectively. The spectral lines are created viewing face-on (upper row) and edge-on (lower row). The black dashed lines show the original spectral line, and the red line has the addition of white noise assuming a S/N of 10, and an increased velocity bin size to 500 m s−1.

4.4 Comparison with observations

Although a one to one comparison between observations and our simulations is not within the scope of this paper, we can still compare the global trends observed. First, the strength of the lines is comparable with observations, leading to a general agreement of the mass-loss rates in Sect. 4.1. A direct comparison can only be made with oxygen rich targets whose mass-loss rates are close to our input mass-loss rate of 10⁻⁷ M_⊙ yr⁻¹. To this end, we can compare to W Hya, T Mic, Y Scl, V1943 Sgr and BK Vir from Ramstedt et al. (2020), who have mass-loss rates ranging between 8 x 10⁻⁸ and 1.5 x 10⁻⁷ M_⊙ yr⁻¹, with terminal velocities around 5 km s−1. After adjusting the observations for the difference in their distance compared to our models (our models are at 500 pc), we find peak flux values for the J = 2→ 1 and 3→2 transitions in the range of 0.5 to 2 Jy, and 1 to 3.5 Jy, respectively. Given the terminal velocity of the observations, we can compare them to the peak fluxes of the v05 models, which fall in the same range as the observations (see Figs. 9 and D.1). The observations in Ramstedt et al. (2020) were made using ALMA, and so there will be resolved out flux, although this should remain minimal in the used configuration.

Moving to the line shapes, most striking is the appearance of the central bump feature, which is observed in targets such as RS Cnc, X Her, AFGL 292, V370 And, and EP Aqr (Danilovich et al. 2015; Díaz-Luis et al. 2019). We note, however, that EP Aqr has been identified as a face-on narrow spiral in Homan et al. (2018, 2020). This is in contrast with our results, where the central bump feature is observed when viewing edge-on. The authors derive a separation of ~65 au, and a companion mass of 0.1 M_⊙, which is well outside the parameter space covered by us. Additionally, these types of spectral lines have been modelled in Díaz-Luis et al. (2019) using an oblate spheroid attached to a biconical structure, viewed nearly face-on. Our models where global flattening occurs do resemble oblate spheroids, but we recover the central feature in the edge-on view. This as a result of the impact of the binary companion on the velocity field, as more material has built up in higher density regions near the central velocities.

A second similarity are the double peaks formed when observing some models face-on. Care has to be taken, since as discussed in Sect. 3.2.5, double peaks can also form as a result of the beam profile. In Díaz-Luis et al. (2019), such double peaked line profiles were again modelled using an oblate spheroid in Díaz-Luis et al. (2019), now viewing edge-on. This is again in contrast with our models, as we recover peaks when viewing at a small inclination, as the origin of the peaks is the additional region towards the z-axis (see the arc-like structures in Fig. C.3). Other targets show also very pronounced double peaks, although in some cases these are known to originate from detached shells, for example R Scl (Maercker et al. 2012, 2024). We cannot yet produce such shells in our simulations, as these are believed to be caused by a brief increase in the mass-loss rate during a thermal pulse.

Besides these stronger features being observed in some systems, some additional subtle results are also retrieved. The emergence of a double-peaked pattern at higher J transitions is observed in V Aql, V CrB and W And (see Figs. A.1 and A.2 in Danilovich et al. 2015). These line profiles also become more asymmetric towards the higher transitions, similar to our results of Sect. 3.2.3. A single peak to the side of the central velocity, similar to the J = 6→5 transitions of the v05a25 model as seen in Fig. 7, can be seen in R Cas (see Fig. 2 De Beck et al. 2010). However, to be similar in nature, we expect the lines to be reasonably optically thick. Smaller scale deviations from flat topped or parabolic lines, similar to the v20 models (see Fig. 3) are also often observed, for example in SW Vir, RX Boo, R And, and RT Vir (Danilovich et al. 2015).

4.5 Limitations and future work

It is important to keep in mind the assumptions made in both the PHANTOM models and the post-processing with MAGRITTE, and the impact these have on our results. As noted in Sect. 2, we neglect the radiative acceleration in our models, and instead opt for a free wind approach. This implies that the velocity profile in our simulations is also unrealistic. As shown in Maes et al. (2021), the wind velocity at the location of the companion strongly impacts the complexity of the resulting morphology, and so obtaining a realistic velocity profile is crucial. In order to further refine the wind launching mechanism, we need to consider detailed dust nucleation (as done in Siess et al. 2022) and radiation pressure (as done in Esseldeurs et al. 2023), though this still needs to be combined with pulsations to obtain physically correct density structures in the inner region, which is also tied to the dust formation. Therefore, we aim to implement a more accurate radiative pressure formalism, and pulsations into PHANTOM, similar to Aydi & Mohamed (2022).

Currently, all our models show very pronounced spiral like structures. Although some observations indicate a well-defined spiral structure (see e.g. AFGL 3068 in Mauron & Huggins 2006; Morris et al. 2006), the majority reveal more chaotic patterns can exist. Including eccentricity is one pathway to introduce more asymmetry (Malfait et al. 2024a). However, detailed 3D RHD simulations of AGB stars, where pulsations and convection emerge self-consistently (see e.g. Freytag & Höfner 2023; Ahmad et al. 2023) show that mass loss can be non-spherically symmetric due to large convection cells and non-radial pulsations. Implementing anisotropies in the particle ejection could have an appreciable effect on the complexity of the overall morphology.

Additionally, we do not yet include detailed cooling mechanisms (see e.g. the mechanisms described in Decin et al. 2006). The effect on the line shapes are expected to be minor, although the line strengths can deviate very strongly with changing temperature, and greatly impact any mass-loss rate retrieval, as discussed above. Hence, we aim to incorporate more cooling mechanisms in our simulations (Dionese et al. In prep). To obtain the most accurate cooling rates, we need detailed chemistry. The first step has been achieved in Maes et al. (2024), who created MACE, a chemistry emulator, which will be coupled to PHANTOM in the future.

Neglecting dust has a second effect on our results, as it also acts as an additional source term in the radiative transfer. Most importantly, dust can affect the level populations by exciting CO to its first vibrational transition, leading to IR pumping (Lamers & Cassinelli 1999). By including the first vibrational level of CO, the amount of transitions that need to be considered increases by a factor of four, leading to a significant increase in computational cost of the radiative transfer. Furthermore, as shown in Matsumoto et al. (2023) dust can also attenuate the CO emission, though they find this plays a larger role for higher transitions. Taking into account dust in the radiative transfer is again not expected to significantly alter the line shapes, but mainly impact the line strengths.

It is important to stress that a very small part of the possible parameter space has been covered, and that a great number of relevant parameters are taken constant over the grid, most importantly the mass-loss rate, which we justified in Sect. 2.2. Increasing the mass-loss rate will greatly alter the photodissociation radius and line strengths, but given the lack of additional cooling and a detailed description of radiation pressure in the current set-up, the impact on the underlying morphology is expected to remain small. Besides that, we only considered oxygen-rich AGB stars by imposing the corresponding fractional abundance of CO. Increasing this value will strongly alter the photodissociation radius, optical depth for the radiative transfer, and line strengths.

This study focusses on the spectral lines, which is the most often observed feature from AGB outflows. However, a great deal of information can easily be concealed, either in the intrinsic shape of the line or by the noise and beam profile. With the emergence of high-resolution interferometric observations, we can also study different observables, such as channel maps and position-velocity diagrams. This has been studied in Homan et al. (2015, 2016), but again with the more simplistic analytical models.

5 Summary and conclusions

Using the SPH code PHANTOM, we created a grid of nine 3D models of AGB stars with a stellar companion, with different orbital separations and initial wind outflow velocities. The simulations shows distinct AGB outflows exhibiting deviations from spherical symmetry, both in the inner region, where a spiral structure arises in the orbital plane, and in the global morphology where flattening can occur. We implemented a novel method to take into account CO photodissociation in complex 3D morphologies that utilises a ray tracing scheme combined with nearest ray interpolation.

Utilising the 3D NLTE line radiative transfer code MAGRITTE, we computed the level populations of CO considering only rotational transitions. We created synthetic spectral lines of the six lowest transitions of CO for all our models viewed at different inclinations and position angles. A variety of line shapes can emerge, depending on the different input parameters. These includes lines that deviate strongly from the single star scenario, with two strong peaks near the terminal wind velocity when viewed face-on, or pronounced bumps near the central velocity when viewed edge-on. In contrast, quasi parabolic or flat-topped profiles can emerge, where the impact of a companion goes unnoticed. The line strengths, and the evolution of the line when looking at different transitions, or at different values for the PA, can provide useful information on the underlying morphology.

We mimicked the effect of observing the models with single-dish telescopes in two ways, by applying a beam profile, and by adding noise to the spectral lines. This showed that many of the features of the spectral lines that are characteristic of a binary companion can be easily hidden by these instrumental effects, with the possible danger of misinterpreting the line profiles as originating from spherically symmetric outflows. This shows the difficulty of obtaining information about the system when only considering the spectral lines. Even so, the companion can have a significant influence on the temperature profile, and hence on the level populations of the model, thereby impacting the retrieved mass-loss rate.

Given the sensitivity of the spectral lines on the stellar and orbital input parameters, the lack of detailed cooling and heating mechanisms, and insufficient wind driving physics, a direct comparison of our synthetic spectral lines to observations is not yet possible. We aim to improve upon these shortcomings in the future and to move towards this goal.

Acknowledgements

OV acknowledges funding from the Research Foundation – Flanders (FWO), grant 1173025N. ME acknowledges funding from the FWO research grants G099720N and G0B3823N. TC acknowledges funding from the Research Foundation – Flanders (FWO), grant 1166724N. LS is a senior researcher for the FNRS. KM acknowledges funding from the Research Foundation – Flanders (FWO), grant 1169822N. FDC is a Post-doctoral Research Fellow of the Research Foundation – Flanders (FWO), grant 1253223N. TD is supported in part by the Australian Research Council through a Discovery Early Career Researcher Award (DE230100183). CL acknowledges support from the KU Leuven C1 excellence grant BRAVE C16/23/009. LD acknowledges support from the KU Leuven C1 excellence grant BRAVE C16/23/009, KU Leuven Methusalem grant SOUL METH/24/012, and the FWO research grants G099720N and G0B3823N.

Appendix A Comparison with SKIRT

To further validate and check the accuracy of our results, we compare them by generating spectral lines with SKIRT⁴ (Baes et al. 2003), an open-source 3D Monte Carlo radiative transfer code. Although SKIRT is mainly used for dust radiative transfer, it allows for the detailed NLTE calculations necessary for line radiative transfer (Baes et al. 2022; Matsumoto et al. 2023).

We carry out the comparison using the vl0a09 model. The results are presented in Fig. A.1, where we show both face-on and edge-on line profiles of the CO J = 1→0,3→2 and 5→4 lines using MAGRITTE (different shades of blue for different transitions) and SKIRT (red). In the bottom plot we show the relative difference between the line profiles for the different transitions, calculated with Eq. (B.2), where the spectral line from MAGRITTE is taken to be exact. We see that for both the J = 3→2 and 5→4 lines, we find nearly perfect agreement, as the relative difference is around one percent for the entire spectral line. Near the wings the difference increases due to the lower values for the flux in that region.

For the J = 1→0 line, SKIRT has a systematically higher flux than MAGRITTE for both the face and edge-on views, where the line profiles differ by ~5−10%, which is still comparable with the expected uncertainty on the lines from observations. Again, the relative error increases towards the wings, due to the lower flux. The line shape still remains the same between the two codes even for the J = 1→0 line. It should be noted that SKIRT was unable to converge completely, which is believed to be linked to this discrepancy. A difference is also present at the J = 2→l line, but to a lesser extent. However, the agreement between the two codes is very good in general, especially towards higher transitions.

Fig. A.1 Comparison of the CO J = 1→0,3→2, and 5→4 lines of the vl0a09 model, between MAGRITTE (different shades of blue) and SKIRT (red) for the face-on (left) and edge-on (right) views. The bottom plot shows the relative difference between the spectral lines.

Appendix B CO photodissociation

In this section we explain the implementation of CO photodissociation, which is inspired by Groenewegen (2017), and comment on the accuracy of the interpolation scheme.

B.1 Implementation

Our implementation works as follows:

• Starting from the surface of the star, we trace a predetermined amount of rays outwards, until the edge of the simulation is reached. These rays are chosen such that they are uniformly distributed in 3D. This is accomplished using the HEALPix⁵ package (Górski et al. 2005), which discretises a sphere into pixels of equal area. The unit vectors, or directions of the rays, are such that they point towards the centre of these pixels. We use HEALPix order of 4, corresponding to 3072 rays.

• Along each ray, at all points along that ray (which are determined in a similar way as the ray tracing scheme of De Ceuster et al. 2019), we calculate the H₂ and CO column densities by tracing rays in all directions, to the edge of the simulation. These rays are again distributed using HEALPix. In this case, we use HEALPix order 1, corresponding to 48 rays for each point.

• Using the calculated column densities, we derive the photodissociation rate using Eqs. (19) and (20). When the value of I(r) is known at all the points along a ray, we compute x(r) along that ray using Eq. (18).

• Once this is known for all the initially traced rays, we can interpolate between these to obtain the value of x_i for any other point i of the simulation. This is done by finding the eight closest rays to the point i, and interpolating the value of x on each neighbouring ray at the distance r_i. The final value of x_i at position r_i is then calculated by using $$x_i=\frac{1}{{\displaystyle {\sum }_{j=1}^8{d}_j^{-2}}}{\displaystyle \sum _{j=1}^8\frac{x_j}{d_j^2},}$$(B.1)

  where x_i is the interpolated value along each ray j, and d_j represents the perpendicular distance between the point i and the ray j, i.e. closer rays have a larger weight.

• We repeat this process, until the relative difference between subsequent iterations becomes less than 0.1 percent, which is typically achieved in five iterations.

B.2 Validation of interpolation

In order to justify the interpolation approach, we trace rays in four random directions in the simulation, calculate x(r) explicitly using Eqs. (18)–(20), and also compute x(r) using the interpolation scheme outlined above. This test has been carried out for the vl0al5 model, and the result is shown in Fig. B.1, where the top plot shows x(r) for the four rays considered, with the blue dots corresponding to the explicit calculation, with different shades showing different directions, and the red dots shows the values of the interpolation. The bottom plot shows the relative difference between the explicit and interpolated values at each point, calculated as $${\delta }_{\text{rel}}=\frac{\left|x_{\text{ex}}\left(r\right)-x_{\text{in}}\left(r\right)\right|}{x_{\text{ex}}\left(r\right)},$$(B.2)

with x_ex(r) the explicitly calculated values and x_in(r) the interpolated values. We see that in general, the red dots overlap with the blue dots, showing the interpolation is accurate. The bottom plot shows that the relative error remains under 1 percent for the majority of the points. Only towards the outer regions does the relative difference start increasing, which is due to the decrease of x(r) itself. This region is also less important as it does not contribute to the line profiles, since x(r) is close to zero.

Increasing the number of rays that are traced, either in the amount of rays that can be used for interpolation, or in the number of rays used to calculate the column densities, improves these results marginally, but at an additional computational cost. We did not find sufficient improvement in the accuracy to justify this increase in time.

Fig. B.1 Top: x(r) for rays in four different directions in the vl0al5 model. The blue dots show the explicit calculation of x(r), and red the interpolation between the eight nearest rays. Bottom: Relative difference between the explicitly calculated and interpolated values. The red dashed line corresponds to a relative difference of 1%.

Appendix C Morphologies and velocities

This section contains additional slices through the orbital and meridional plane for the v05 and v20 models, for both the density and velocity structures, analogous to Figs. 1 and 4 for the v10 models.

Fig. C.1 Same as Fig. 1, but for the models with initial wind velocities of 5 km s−1.

Fig. C.2 Same as Fig. 1, but for the models with initial wind velocities of 20 km s−1.

Fig. C.3 Same as Fig. 4, but for the models with initial wind velocities of 5 km s−1.

Fig. C.4 Same as Fig. 4, but for the models with initial wind velocities of 20 km s−1.

Appendix D Additional line profiles

D.1 Intermediate inclinations

This section contains additional plots, showing the synthetic CO J = 2→1 line under intermediate inclinations for the v05 and v20 models, analogous to Fig. 6 for the v10 models.

Fig. D.1 Same as Fig. 6, but for the models with initial wind velocities of 5 km s−1.

Fig. D.2 Same as Fig. 6, for the models with initial wind velocities of 20 km s−1.

D.2 Higher S/N

This section contains an additional figure, visualising the observability of features when assuming a higher S/N value of 20.

Fig. D.3 Synthetic CO spectral lines for three selected models. From top to bottom these are v05a09, v10a25, and v20a25 with the J = 2→1, 4→3, 6→5 transitions, respectively. The spectral lines are created viewing face-on (left column) and edge-on (right column). The black dashed lines show the original spectral line, and the red line has the addition of white noise assuming a S/N of 20, and an increased velocity bin size to 500 m s−1.

References

All Tables

All Figures

	Fig. 1 Density distributions in slices through the orbital plane (upper row) and the meridional plane (lower row) for the models with initial wind velocities of 10 km s−1. From left to right, the orbital separation is 9, 15, and 25 au.
In the text

	Fig. 2 Number density of CO in slices through the meridional plane after accounting for photodissociation. The solid and dotted lines represent the location where x(r) is 0.5 and 0.01, respectively.
In the text

	Fig. 3 Synthetic CO J = 2→1 spectral lines for all models. The columns correspond, from left to right, to initial outflow velocities of 5, 10, and 20 km s−1. The solid, dashed, and dotted lines represent orbital separations of 9, 15, and 25 au. The spectral lines are creating viewing face-on at an inclination of 0° (upper row) and edge-on at 90° (lower row).
In the text

	Fig. 4 Velocities in the z direction (upper row) and in the x direction (lower row) in slices through the meridional plane for the models with initial wind velocities of 10 km s−1. Contours are shown on the plots around velocities interesting for the line profiles, indicated by the horizontal lines on the colour bar. From left to right, the orbital separation is 9, 15, and 25 au.
In the text

	Fig. 5 Synthetic CO spectral lines of the first six transitions from a single-star PHANTOM model, with a mass-loss rate of 10−7 M⊙ yr and an initial wind velocity of 10 km s−1.
In the text

	Fig. 6 Synthetic CO J = 2→1 spectral lines for the models with initial wind velocities of 10 km s−1, viewed under inclinations of 18°, 36°, 54°, and 72°. The solid, dashed, and dotted lines represent orbital separations of 9, 15, and 25 au.
In the text

	Fig. 7 Synthetic CO spectral lines of the first six rotational transitions, for the models with initial outflow velocities of 5 km s−1, viewed under an inclination of 0°. The solid, dashed, and dotted lines represent orbital separations of 9, 15, and 25 au.
In the text

	Fig. 8 Synthetic CO spectral lines of the first six rotational transitions, for the models with initial outflow velocities of 20 km s−1, viewed under an inclination of 90°. The solid, dashed, and dotted lines represent orbital separations of 9, 15, and 25 au.
In the text

	Fig. 9 Synthetic CO J = 1→0 (top) and J = 6→5 (bottom) spectral lines for the v10 (left) and v20 (right) models, viewed under an inclination of 0°. The solid, dashed, and dotted lines represent orbital separations of 9, 15, and 25 au.
In the text

	Fig. 10 Synthetic CO J = 1→0 (top) and J = 6→5 (bottom) spectral lines for the v05 (left) and v10 (right) models viewed under an inclination of 90°. The solid, dashed, and dotted lines represent orbital separations of 9, 15, and 25 au.
In the text

	Fig. 11 Temperature distributions in slices through the orbital plane (upper row) and the meridional plane (lower row) for the v05a09 (left) and v20a25 (right) models.
In the text

	Fig. 12 Synthetic CO J = 6→5 spectral lines for the v05a09 (left) and v20a25 (right) models, viewed edge-on and at four different values of the PA (ϕ). The different colours correspond to the different values of the PA.
In the text

	Fig. 13 Synthetic CO J = 2→1 spectral lines for the models with initial wind velocities of 10 km s−1, where a beam with five different values of θb has been applied. The spectral lines are viewed face-on at an inclination of 0° and edge-on at 90° in the upper and lower row, respectively. From left to right, the orbital separation is 9, 15, and 25 au.
In the text

Fig. 14 Synthetic CO spectral lines for three selected models. From left to right these are v05a09, v10a25, and v20a25 with the J =2→1, 4→3, 6→5 transitions, respectively. The spectral lines are created viewing face-on (upper row) and edge-on (lower row). The black dashed lines show the original spectral line, and the red line has the addition of white noise assuming a S/N of 10, and an increased velocity bin size to 500 m s−1.

In the text

	Fig. A.1 Comparison of the CO J = 1→0,3→2, and 5→4 lines of the vl0a09 model, between MAGRITTE (different shades of blue) and SKIRT (red) for the face-on (left) and edge-on (right) views. The bottom plot shows the relative difference between the spectral lines.
In the text

	Fig. B.1 Top: x(r) for rays in four different directions in the vl0al5 model. The blue dots show the explicit calculation of x(r), and red the interpolation between the eight nearest rays. Bottom: Relative difference between the explicitly calculated and interpolated values. The red dashed line corresponds to a relative difference of 1%.
In the text

	Fig. C.1 Same as Fig. 1, but for the models with initial wind velocities of 5 km s−1.
In the text

	Fig. C.2 Same as Fig. 1, but for the models with initial wind velocities of 20 km s−1.
In the text

	Fig. C.3 Same as Fig. 4, but for the models with initial wind velocities of 5 km s−1.
In the text

	Fig. C.4 Same as Fig. 4, but for the models with initial wind velocities of 20 km s−1.
In the text

	Fig. D.1 Same as Fig. 6, but for the models with initial wind velocities of 5 km s−1.
In the text

	Fig. D.2 Same as Fig. 6, for the models with initial wind velocities of 20 km s−1.
In the text

Fig. D.3 Synthetic CO spectral lines for three selected models. From top to bottom these are v05a09, v10a25, and v20a25 with the J = 2→1, 4→3, 6→5 transitions, respectively. The spectral lines are created viewing face-on (left column) and edge-on (right column). The black dashed lines show the original spectral line, and the red line has the addition of white noise assuming a S/N of 20, and an increased velocity bin size to 500 m s−1.

In the text