© The Authors 2024

1 Introduction

At the end of the core-He burning phase, evolved low- and intermediate-mass stars enter the asymptotic giant branch (AGB) phase. During this phase, the stars expand and cool down, and lose their outer envelopes through pulsation-enhanced dust-driven winds. Afterwards, they evolve through the post-AGB and planetary nebulae (PNe) evolutionary phases to become white dwarfs. Many of these AGB successors host binary or higher-order (more than two stars) stellar systems (see e.g. De Marco et al. 2013; Oomen et al. 2018; Escorza & De Rosa 2023), and have very complex morphologies that are believed to be related to this multiple nature (Jones & Boffin 2017; Akashi & Soker 2017; De Marco 2009; García-Segura et al. 2018; Corporaal et al. 2023; Verhamme et al. 2024). Similarly, AGB outflows have recently been observed with high angular and spectral resolution, revealing non-spherically symmetric circumstellar environments in which complex structures such as arcs, spirals, discs, and bipolar outflows have been discovered, with the main shaping mechanism of these complexities believed to be binarity (Mauron & Huggins 2006; Maercker et al. 2012; Decin et al. 2020; Danilovich et al. 2024). Moreover, both in PNe and in AGB outflows, the complexity of some observed morphological structures suggests that they may be created by the impact of more than one companion. An example is the AGB star R Aql, which has been observed within the ATOMIUM ALMA (Atacama Large Millimeter Array) large programme, and which has a circumstellar envelope (CSE) with a lot of intricate spiral- and arc-like structures (Decin et al. 2020). Recent theoretical and observational studies have strengthened the hypothesis that several PNe are shaped by triple stellar systems (Akashi & Soker 2017; Jones et al. 2019; Schreier et al. 2019; Miszalski et al. 2019). According to population synthesis studies, ~50% of low-mass stellar systems have at least two stellar components, and 13% have three or more (Offner et al. 2023; Tokovinin 2014; Gao et al. 2022). Triple systems often tend to be in a stable hierarchical configuration, meaning that they consist of an inner binary and a distant third component (Hut & Bahcall 1983; van den Berk et al. 2007; Tokovinin 2021). Although in a hierarchical setup, many triple systems have very complex orbital mechanics, with for example inclined and eccentric orbits in which Kozai-Lidov cycles are at play (Lidov 1962; Kozai 1962; Mardling & Aarseth 2001; Naoz 2016). In this work, we take a first step in understanding the effect of a triple stellar system on the structure formation of AGB outflows. We investigate this with the use of smoothed-particle hydrodynamic (SPH) models of triple systems including a wind-launching AGB star. For this first study, we limit ourselves to the most elementary triple systems that are co-planar and in stable hierarchical orbits, such that complex mechanisms like unstable N-body orbits and Kozai–Lidov cycles are not relevant.

This paper builds on previous studies in which the wind morphology in various binary AGB systems has been investigated in detail with the use of different hydrodynamic models, including Kim et al. (2017); Liu et al. (2017); Chen et al. (2017); Saladino et al. (2018); El Mellah et al. (2020); Maes et al. (2021); Malfait et al. (2021); Esseldeurs et al. (2023); and Malfait et al. (2024). We introduce the input physics of our SPH model, and the set-up of our grid, consisting of triple and binary systems, in Sect. 2. We first investigate the wind morphology in our set of binary systems in Sect. 3, before analysing the impact of adding a third body to these simulations in Sect. 4. This analysis includes a study of the impact of varying the wind velocity and orbital eccentricity on the wind morphology. The source of inspiration to study the morphology of triple systems is the AGB star R Aql. We therefore present a short study of ALMA observations of R Aql and compare these to a synthetic simulation of one of our triple hydrodynamic simulations in Sect. 5. Finally, our conclusions are provided in Sect. 6.

2 Method

2.1 Input physics SPH model

We performed hydrodynamic simulations of an AGB outflow in a hierarchical triple system with the SPH code PHANTOM (Price et al. 2018). We used the same numerical input physics as is described in Malfait et al. (2024), where more details can be found. In short, the AGB star and companions were modelled as gravity-only sink particles that can accrete the mass and momentum of gas particles that fall within 0.8 to 1 times the set-up accretion radius (see Price et al. 2018; Malfait et al. 2024). In wind simulations, SPH particles are launched with an initial velocity v_ini, on concentric spheres around the primary AGB star (Siess et al. 2022). They are accelerated by a pressure gradient and with the use of the free-wind approximation, in which the gravity of the AGB star is neglected by the wind particles, as such artificially balancing the gravity by the radiation force (Theuns & Jorissen 1993; Esseldeurs et al. 2023). H I cooling, following Spitzer (1978), is included in the energy equation, as is explained in Malfait et al. (2024). We adopt constant values of the polytropic index γ = 1.2 and the mean molecular weight µ = 2.38, modelling a molecular wind with a typical, radially decreasing, power-law temperature profile for AGB winds that strongly depends on the heating and cooling mechanisms included in the energy equation (Millar 2004; Maes et al. 2021).

Fig. 1 Set-up of hierarchical triple with inner binary M1 –M2 orbiting around ${\text{CoM}}_{M_1+M_2}$ and outer companion M3 orbiting around the system’s centre of mass ${\text{CoM}}_{M_1+M_2+M_3}$. M3 induces an orbital motion of ${\text{CoM}}_{M_1+M_2}$, and thereby also of M1 and M2, around the system’s centre of mass.

2.2 Grid set-up triples

Fig. 1 illustrates the hierarchical set-up used in our models. It consists of an inner binary given by M₁−M₂ and a third companion M₃. The third star and the centre of mass (CoM) of the inner binary $\left({\text{CoM}}_{M_1+M_2}\right)$ form the outer binary, orbiting the CoM of the entire system $\left({\text{CoM}}_{M_1+M_2+M_3}\right)$. Table 1 provides the set-up parameters that are the same for all models. The inner binary, M₁–M₂, is set up with a mass ratio q₁ = M₂/M₁ = 0.25 in which M₁ = 1.6 M_⊙ is the mass of the AGB star and M₂ = 0.4 M_⊙ is the mass of the close companion. The primary has a mass-loss rate of Ṁ = 1.1 × 10⁻⁶M_⊙ yr⁻¹ and surface temperature of T₁ = 2800 K. These parameters are chosen to correspond to those of the AGB star R Aql for which high-resolution observations reveal a very intricate wind morphology that might be created by the presence of multiple companions (Decin et al. 2020, see also Sect. 5.1).

Table 2 summarises the properties of all our triple simulations. We start by varying the initial wind velocity, v_ini, between 8 and 15 km s⁻¹, as this parameter is known to strongly impact the outflow morphology (Maes et al. 2021; Malfait et al. 2021, 2024). Secondly, we vary the mass of the third companion (M₃) to 0.1, 0.6, and 1.2 M_⊙ to analyse its gravitational influence on the wind. Finally, we investigate the impact of eccentricity on the morphologies in triple systems by introducing an orbital eccentricity of e_i = 0.4 and e_o = 0.4 once in the inner and the outer orbit, respectively. We do this for the triple system v15m06 (v_ini = 15 km s⁻¹, M₃ = 0.6 M_⊙), which will be used as a reference model for comparison. We select this model because it has a rather regular, not too complex morphology due to the moderate tertiary mass and high wind velocity.

The orbital period of the inner binary, M₁ –M₂, is P_i = 7.91 yr, and the outer binary orbital period is P_o = 143.0, 128.5, and 115.8yr for masses of M₃ = 0.1,0.6, and 1.2 M_⊙, respectively. All systems are modelled during at least four orbital periods of the outer binary. The eccentric systems run until they reach the fourth apastron passage, which occurs at different times in the simulation depending on the initial configuration.

To test if the systems are expected to be dynamically stable, we adopt the commonly used stability criterion of Mardling & Aarseth (2001) for co-planar hierarchical triple systems. This criterion predicts a system to be stable if the ratio of the semi-major axes of the outer to the inner orbit (a_o/a_i) is larger than the critical ratio (a_o/a_i)|_crit, $${\frac{a_{\text{o}}}{a_{\text{i}}}|}_{\text{crit\hspace{0.17em}}}=\frac{2.8}{1-e_{\text{o}}}{\left(\frac{\left(1+q_{\text{o}}\right)\cdot \left(1+e_{\text{o}}\right)}{\sqrt{1-e_{\text{o}}}}\right)}^{2/5},$$(1)

in which the mass ratio of the outer orbit is defined as q_o = M₃/(M₁ + M₂). This critical ratio is given in Table 2 and is lower than the ratio a_o/a_i = 7.0 for all modelled systems, indicating that the systems are expected to be dynamically stable. As an additional test, the orbit of all systems is integrated with the REBOUND N-body code (Rein & Liu 2012) over 1000 outer orbital periods, P_o, to confirm the stability of the orbits. Note that this test does not take into account the effect of mass-loss and transfer of mass and angular momentum to the companions, which might influence the orbits of this system during the AGB phase.

2.3 Grid set-up binaries

To analyse the impact of adding a tertiary component to a binary system, we first modelled a limited set of binary models that are subsystems of the hierarchical triple, as is described in Fig. 1 and Tables 1 and 2. The set-up of the binary models is given in Table 3. Firstly, we investigated the wind morphology of the inner binary (M₁ –M₂), for the two different initial wind velocities of 8 and 15 km s⁻¹ (models v08bin_i and v15bin_i), to which an outer companion, M₃, was added in the triple simulations. Additionally, we modelled one circular and one eccentric (e = 0.4) outer binary system (M_p = M₁ + M₂)–(M₃) to study the effect of splitting up the primary into a close binary (models v15bin_om06 and v15bin_om06e40).

2.4 Classification parameter and flattening ratio

The strength of the wind-companion interaction and the expected deviation from a spherically symmetric morphology can be estimated with the use of a classification parameter, ε, proposed by Maes et al. (2021) and applied by Malfait et al. (2021). In a binary system, ε is defined as the ratio of gravitational energy density around the companion to the kinetic energy of the wind, $$\epsilon =\frac{e_{\text{grav}}}{e_{\text{kin}}}=\frac{{\left(24G^3{M}_{\text{s}}^2{M}_{\text{p}}\right)}^{1/3}}{v_{\text{w}}^{2}a(1-e)},$$(2)

with v_w an estimate of the wind velocity at the location of the companion. This velocity is approximated by the vector sum of the wind velocity at the orbital separation of the companion in a 1D (spherically symmetric) simulation without a companion and the average orbital velocity of the primary star, to account for the effect of the orbital motion on the wind velocity field. As such, v_w is $$v_{\text{W}}=\sqrt{v_{1\text{D}}{(r=a)}^2+v_{\text{orb},\text{p}}^2}.$$(3)

We anticipate that when ε ≪ 1, a quasi-spherical outflow with a thin, low-density contrast spiral feature will result; for ε ~ 1, a regular Archimedes spiral will form; and for ε ≫ 1, a more complex asymmetrical morphology will form (Maes et al. 2021). For our triple systems, we calculated ε on two hierarchical levels, for the inner system, M₁–M₂ (ε_i), and for the outer binary system, (M₁ + M₂)–M₃ (ε_o). The values are provided in Table 2 and are used throughout our analysis.

Further, it is known that the impact of a binary companion results into a flattening, or elongation of the morphology, by the induced orbital motion of AGB star (Maes et al. 2021; Malfait et al. 2021). As is explained by Malfait et al. (2021), we theoretically expect that the AGB star launches wind material with a velocity of v_w = v_ini + v_orb,p, the vector sum of the input wind velocity and AGB orbital velocity. Because the orbital velocity vector lies in the orbital plane, the resulting wind has a larger velocity component in the orbital plane direction, which results in an elongation of the morphology. As is illustrated by Malfait et al. (2021), this flattening can be quantified by the ratio, R_f, of the height to the width of the edge-on wind structures. Theoretically, when neglecting the impact of the direct gravitational attraction of matter by the companion on the shape of the arc-structures, this flattening ratio is expected to be $$R_{\text{f}}=\frac{v_{\text{ini\hspace{0.17em}}}}{v_{\text{ini\hspace{0.17em}}}+v_{\text{orb,p}}},$$(4)

with a lower ratio, R_f, predicting the formation of more elongated wind structures. The ratios for the circular binary models are provided in Table 3, and are used in Sect. 3.

Fig. 2 Density (top) and velocity (bottom) distribution in a slice through the orbital (left) and meridional (right) plane of model v08bini for r ≤ 475 au.

3 Wind morphology in binary systems

Before analysing the impact of a third companion on the wind morphology, we investigate the outflow of the binaries to which our triple systems will be compared. First, we discuss the inner binary without outer companion, for initial wind velocities of v_ini = 8 and 15 km s⁻¹ (models v08bin_i and v15bin_i). Figs. 2 and 3 reveal the wind morphology of these systems, showing the density and velocity distribution in slices through the orbital (left panels) and meridional (right panels) planes. As is expected from the ε values close to 1 (Table 3), in both models an Archimedean spiral structure forms in the orbital plane, and concentric arcs appear in the edge-on plane. For model v15bini shown in Fig. 3, these arcs resemble ‘banana’-like shapes with narrow ends that stretch towards the poles and thicker two-edged regions around the orbital plane (z = 0). The thicker part of the banana-shapes emerges from the two-edged gravity wake of the companion, which is shaped by the orbital motion of the AGB star and the direct gravitational attraction of the companion on the wind, while the narrow ends, reaching up to the poles, result solely from the orbital motion of the stars (Kim & Taam 2012). The transition between the thicker centres and narrower edges is visible as a straight line at an angle θ ~ 50° with respect to the orbital plane. The stronger the relative importance of the direct gravitational impact of the companion (compared to the impact of the orbital motion of the AGB), the larger the cross-section of the gravity wake will be, and the larger the angle of this transition line is expected. For model v08bini, in which the lower velocity causes a stronger gravitational influence of the companion on the wind (larger ε value, Table 3), the distinction between the thin edges and broader centres of the arc-like shapes can only be made for the innermost arcs at r ≲ 50 au. Further out, the cross-section of the gravity wake (so the thicker part of the arcs) expands up to the poles. As a result, the lines resulting from the transition between the thin edges and broader centres, are not linear, but follow an increasing angle, θ, with respect to the orbital plane for increasing r. When the edges reach the poles (θ = π/2), they interact with the structures from the opposite side (at x = 0, y ≈ ± 120 au).

Comparing the morphology in the meridional plane of both models, the arc structures in model v08bini are more elongated than in model v15bin_i. In these models, the primary’s orbital velocity, which causes a flattening of the wind morphology (see Sect. 2.4), is v_brb,p = 3.77 km s⁻¹, which is relatively low due to the low secondary mass M₂ = 0.4 M_⊙. The edge-on velocity distribution of model v08bin_i in Fig. 2 (lower right panel) illustrates that there is a flattening due this orbital motion, as the velocity |v| within the arc-structures decreases with increasing angle away from the orbital plane, with |v| ≈ 12 km s⁻¹ in the orbital plane, to |v| ≈ 6 km s⁻¹ close to the poles. The theoretical flattening of these models, from Eq. (4), resulting solely from the impact of the orbital velocity on the velocity distribution in the wind, are 0.68 and 0.80, for models v08bin_i and v15bin_i, respectively. The actual flattening ratio in the simulations is estimated by measuring the height-to-width ratio of the edge-on arc-like structures in Figs. 2 and 3, and are ~0.6 and ~0.8 for v08bin_i and v15bin_i. This shows that the theoretical prediction of the flattening (which does not take into account the direct gravitational effect of the companion) is somewhat underestimating the flattening in model v08bin_i, but is a good approximation for model v15bin_i.

Next, we present the outer binaries v15bin_om06 and v15bin_om06e40, which consist of a primary AGB star with M_p = M₁ + M₂ and outer companion M_s = M₃ = 0.6, an input wind velocity of v_ini = 15 km s⁻¹, and eccentricities of 0 and 0.4, respectively. Because of the large orbital separation, and thereby lower density, ρ, at the location of the companion, its gravitational impact is expected to be smaller, which is in agreement with the smaller ε < 1 (Table 3). Fig. 4 shows the morphology of model v15bin_om06. In the orbital plane, it consists of an extended face-on two-edged Archimedean spiral density enhancement on top of a spherically symmetric wind. In the meridional plane, there are edge-on arcs with thin ends reaching towards the polar regions and thicker two-edged structures in the centre that reach up to only θ ~ 25° with respect to the orbital plane, due to the reduced gravitational pull of the companion. The orbital velocity of the primary star, v_orb,p = 1.87 km s⁻¹, is lower than in the inner binary models due to the large orbital separation. As a consequence, the morphology is less elongated, with a theoretical flattening ratio R_f = 0.89 (Eq. (4)) and a measured ratio R_f,sim ≈ 1. Additionally, due to the large orbital separation and relatively high wind velocity, the spiral structure travels a significant distance during one orbital period, such that the distance between two consecutive arcs or spiral arms is larger compared to shorter period binary models. Unlike in the inner binary models, the structures do not interact and overlap in the spatial range of the simulation, and the quasi-spherical wind is still observable between the 3D structures. In the snapshot of Fig. 4, the inner edge of the spiral wake crosses the x axis at x ≈ 450, 900 and 1350 au, indicating that the structure travels ~450 au radially outward during one orbital period of 128.5 yr, thereby moving with a velocity of approximately 16.6 km s⁻¹. In observational studies, the distance between two recurrent high-density structures can in this way be used to estimate the orbital period or wind velocity, under the condition that a good estimate of the other quantity is available.

Fig. 5 exhibits the density distribution in the orbital and an edge-on plane of the eccentric system v15bin_om06e40. The inclusion of eccentricity results in a globally more asymmetric morphology. At the apastron side of the companion (x > 0), the minimal orbital velocity results in a widening of the gravity wake, and at periastron side (x < 0), the higher orbital velocity makes the wake thinner. As a consequence, the edge-on arcs at apastron side consist mainly of broad two-edged structures centred around the orbital plane, and at periastron side of large, thin concentric arcs.

Fig. 3 Same as Fig. 2, but for model v15bini.

Fig. 4 Density distribution in a slice in (top) and perpendicular to (bottom) the orbital plane of model v15binom06.

4 Wind morphology in triple systems

We study the wind morphologies in hierarchical triple systems, using the set of models as defined in Fig. 1 and Tables 1 and 2. We focus on the impact of changing the triple companion mass, M₃, the wind velocity, v_ini, and the eccentricity, e. We expect from theory that a lower wind velocity, increased companion mass, and higher eccentricity will make the wind deviate stronger from spherical symmetry, which is reflected in the values of ε, given in Table 2. Figs. 6 and 7 display the density distribution of the non-eccentric models with outer companion masses M₃ = 0.1, 0.6, 1.2 M_⊙, in the orbital plane (top) and a meridional plane (bottom) for input velocity of v_ini = 15 km s⁻¹ and v_ini = 8 km s⁻¹, respectively. Additionally, a more zoomedin view of the density and velocity distribution of these models are shown in Figs. 8 (v15) and 9 (v08). The eccentric models are represented in Figs. 10, 11, 12, and 13.

Fig. 5 Density distribution in a slice in (top) and perpendicular to (bottom) the orbital plane of model v15binom06e40. The apastron passage of the companion is on the y = 0, x > 0 axis

4.1 Impact of a hierarchical third body

We study the impact of hierarchically adding a third body to the binary systems, described in Sect. 3, from two different points-of-view: (i) considering the close binary model v15bin_i where a third, outer companion is added, and (ii) considering the wide binary model v15bin_om06, where the primary is split up into a close binary. First, we consider the effect of adding a third, outer body. Comparing the density distribution in the triple models in Fig. 6 to the inner binary v15bin_i in Fig. 3, we see that the face-on spiral and edge-on circular arcs, created by the inner binary M₁ –M₂, are still present and resolved, independent of the outer companion mass M₃. On top of this compact underlying structure, there is an additional gravity wake in the orbital plane originating from the companion M₃, which appears as large concentric arcs in the meridional plane. We find that the more massive M₃ and thus the larger ε_o (from left to right in Fig. 6), (i) the larger the density enhancement of its gravity wake and arcs with respect to the underlying wind and compact spiral pattern of M₁ –M₂, and (ii) the more radially compressed its spiral and arcs are, because its orbital period is shorter, and the radial velocity within its gravity wake is lower. Although ε_o is very small for all outer binary systems, the outer companion M₃ is strongly altering the global 3D morphology.

Fig. 8 reveals in more detail the onset of the gravity wake spiral originating from M₃, through a more zoomed-in view of the density and velocity distribution in the orbital and edge-on plane of the triple v15 models. Because the outer companion on its orbit continuously encounters outward moving density- and velocity-enhanced spiral shocks from the inner binary, the outer edge of the M₃ gravity wake has a waved pattern. The orbital plane velocity and density distribution of model v15m06 (Fig. 8, central column) illustrate more clearly how shocks propagate through this gravity wake, resulting in the wavy outer edge. This waved pattern and the shocks within the gravity wake become less prominent further out in the wind. This is less noticeable in the edge-on view, but still visible, mainly in the velocity distribution of model v15m12 (lower right plot in Fig. 8). The companion M₃ also affects the shape and properties of the underlying spiral created by M₁–M₂. Because it induces an orbital motion of M₁–M₂ around the system’s centre of mass CoM_{M1+ M2+ M3}, the orbital velocity of the AGB star and inner companion have an additional component, which affects the velocity and morphology of the M₁–M₂ spiral. The effect is stronger for a more massive outer companion M₃, and is thereby best visible in model v15m12. In the orbital plane and edge-on velocity plots of this model in Fig. 8 (right column), material in front of the inner binary (x ≲ 0, y ≲ 0) has significantly larger velocities than the material directly behind it and further back on its orbit (y > 0 and x < 0). Due to this orbital motion, the M₁–M₂ spiral appears as a snail-shell pattern.

Next, we discuss the effect of splitting the primary star of the wide binary v15bin_om06 into a close binary, resulting in the triple model v15m06. Comparing the M₃ gravity wake in the middle panels in Fig. 6, with the gravity wake in the wide binary in Fig. 4, we see that the interaction of the companion M3 with the inner spiral of M₁–M₂ has several effects. It not only (i) generates the wave-pattern in the outer edge of the M₃ gravity wake, but it also (ii) produces a broadening of this wake, and furthermore, (iii) the presence of the inner binary and its spiral pattern reduce the radial velocity of the M₃ gravity wake. This can be seen by comparing the radial distances travelled by the recurrent arcs and spiral arms in both models, and from the velocity profiles in Fig. 8 (model v15m06, second row, middle column) and Fig. A.1. For example, the inner and outer edge, respectively, of the third spiral and arc on the x < 0 side have travelled about 850 and 1100 au in model v15m06, and about 1100 and 1250 au in the binary model v15bin_om06.

Fig. 6 Density distribution in a slice in the orbital plane (top) and in a perpendicular meridional plane (bottom) for models v15m01 (left), v15m06 (middle), and v15m12 (right).

Fig. 7 Same as Fig. 6, but for models v08m01 (left), v08m06 (middle), and v08m12 (right).

Fig. 8 Density and velocity distribution in a slice in the orbital plane (top) and in a perpendicular meridional plane (bottom) for models v15m01 (left), v15m06 (middle), and v15m12 (right), for r ≤ 475 au.

Fig. 9 Same as Fig. 8, but for models v08m01 (left), v08m06 (middle), and v08m12 (right).

4.2 impact of a different wind velocity

Decreasing the wind velocity in triple systems increases the gravitational impact of both companions (see also lower ε_i,ε_o in Table 2). The effect on the global morphology can be seen by comparing the density distribution of the three models with dif-ferent mass M₃ in Fig. 6 (v15) and Fig. 7 (v08), and at smaller scales in the density and velocity distribution in Fig. 8 (v15) and Fig. 9 (v08). In the edge-on view, there is a stronger flattening in case of lower wind velocity. The cross-section of the spirals produced by M₃ can no longer be described as concentric banana-shaped arcs, but resemble a bit more bicentric elongated structures that extend up to the polar axes. Similar to the binary model v08bin_i (Fig. 2), the edge-on velocity distribution in Fig. 9 reveals the characteristic decrease in velocity with increasing angle with respect to the orbital plane, that appears in flattened morphologies. The flattening ratio of the inner binary model v08bin_i was measured as R_f,sim ≈ 0.6. In the triple v08 mod-els, these ratios are determined by the edge-on arcs originating from the triple companion M₃, as long as the density contrast of these arcs with the underlying wind is large enough to be visible and measurable. This measured R_f is approximately 0.5 for all masses of M₃. This indicates that compared to the binary models, the triple companion slightly increases the elongation of the morphology, but varying its mass does not significantly affect the degree of flattening, at least in models where the density contrast is still visible.

As in the v15 models, the orbital plane morphology of the v08 models is again characterised by a large gravity wake spiral on top of the underlying compact binary spiral pattern. The den-sity distribution in Fig. 9 reveals that the gravity wake of M₃ in models v08m01 and v08m06 is a two-edged spiral, with a wave-pattern in the outer edge. The structure in model v08m12 is different, with a bow shock originating in front of M₃, instead of a two-edged gravity wake behind it. We first focus on comparing models v08m01 and v08m06 in Fig. 9 with their correspond-ing v15 triple models in Fig. 8, before going into more detail about model v08m12. Firstly, the waves in the outer edge of the M₃ gravity wake are less pronounced in the v08 models because the interaction of the inner M₁–M₂ spiral shocks with M₃ and its gravity wake is less strong. This is, firstly, because the inner binary spiral shocks have lower velocities when they collide with the M₃ wake, as can be seen by comparing the orbital plane velocity distribution in Figs. 9 and 8. We note that the colour bars in these figures have different limits. Secondly, the density contrast encountered by M₃ as it passes a M₁–M₂ spiral shock is much smaller than in case of the higher wind velocity models, because in the v08 models, the inner edge of the M₁–M₂ spiral has already caught up with the outer edge at the orbital radius of the companion, as can be seen in the zoomed in density dis-tribution of models v08m06 and v15m06 in Fig. A.2. Thirdly, when the wind has a lower velocity, the M₃ gravity wake is also wider (the distance between the inner and outer edge is larger), and is wrapped more closely around the system CoM due to the lower radial velocity of the material within the wake. As a result, the outer edge of the M₃ gravity wake catches up and merges with the inner edge that originated one orbital period earlier, after, respectively, ~2 and 1 orbital periods in models v08m01 and v08m06. Beyond this interaction radius, the spiral pattern of the M₁ –M₂ binary is no longer visible. From this, we conclude that the lower the wind velocity and the higher the mass M₃, the lower the radius up to where the M₁ –M₂ spiral pattern can be resolved.

Of all the systems, the highest ε_i and ε_o values are found in model v08m12, indicating that the impact of the compan-ions on the wind morphology is strongest. Due to the cumulative effect of the relatively large companion mass M₃ and low wind velocity, the companion M₃ is able to gravitationally attract enough material around it to form an accretion disc and a bow shock. Accretion discs are expected to form around the AGB companion(s), but to resolve them numerically requires very small accretion radii and a very high number of SPH particles (Malfait et al. 2024). These conditions cannot be reached in our triples set-ups because of the large radial extent of the systems. We find an accumulation of wind material around and behind the companion M₃ of model v08m12, such that after ~300 yr of simulation time the wind morphology starts to change. This effect is visible through the high-density region behind the outer companion in the orbital plane density distribution in the upper right panel of Fig. 9. We performed a test simulation of this model with µ = 1.26, which is more suited for atomic gas, and find no accumulation of material. Additionally, this change of µ increases the wind velocity throughout the entire wind, which strongly affects the entire morphology, as was concluded as well by (Malfait et al. 2024). This indicates that a more robust treatment of the equation of state that takes into account the variations in µ and γ is needed.

Fig. 10 Density distribution in a slice in the orbital plane (top) and in a perpendicular meridional plane (bottom) for model v15m06ei40, with apastron passage of the inner companion at x > 0, y = 0, z = 0.

4.3 Impact of eccentricity

To study the characteristics of eccentricity in triple system mor-phologies, an eccentricity of 0.4 is introduced to the inner and outer orbit in models v15m06e_i40 and v15m06e_o40, respec-tively. For the model with an inner eccentric orbit, the global density distribution (Fig. 10) is very similar to the one of the non-eccentric model v15m06 of Fig. 6 (middle row). The largest difference is that at apastron side (x > 0), the spiral wake of M₃ is slightly distorted. This distortion is revealed as an asymmetry in the meridional plane slice (through y = 0, Fig. 10). The edge-on arcs have a different shape at apastron side (x > 0), with a thicker central bulge around the orbital plane, than at periastron side (x < 0). The effect of the eccentric inner orbit is more prominent on smaller scales, where it induces an asymmetry in the compact spiral of M₁–M₂ that is visible between the windings of the outer M₃ spiral wake (see Fig. 11). This asymmetric compact M₁–M₂ spiral deviates from the one in the non-eccentric model v15m06 (Fig. 8), especially at the upper right side of the orbital plane (x > 0, y > 0), in a very similar way as model v20e50 in Malfait et al. (2021, 2024), where the origin of this type of asymmetry is explained in detail. We conclude that the eccentricity of the inner orbit has a limited effect on the large-scale wind structures, but is revealed on a smaller-scale, indicating that this type of eccen-tricity can probably only be detected observationally if the inner wind structures are resolved well.

The impact of the eccentricity on the large-scale wind morphology is more pronounced if it is induced in the outer orbit, as can be seen from the density distribution in the orbital and edge-on planes in Fig. 12 of model v15m06e_o40. In the orbital plane, the outer spiral wake is thicker at the apastron side (x > 0) and more narrow at the periastron side (x < 0). In the meridional plane, the periastron side contains extended arcs reaching up to the polar regions, whereas at apastron side, the structures resem-ble more banana-shapes with thicker two-edged bulges more confined around the orbital plane. The asymmetries induced by eccentricity are less obvious than in the eccentric close binary systems studied by Malfait et al. (2021) and Malfait et al. (2024). To investigate whether this is purely due to the larger extent of the orbit, or due to the added complexities induced by the triple nature of the system, we compare this model (v15m06eo40) to the eccentric binary model v15bin_om06e40 in Fig. 5. We conclude that the previously mentioned asymmetries are more pronounced in the binary model, and are thereby indeed affected by the triple nature of the system. When we split up the primary into a close binary, similar effects result as in the comparison between the circular triple model v15m06 and the wide binary v15binom06 in Sect. 4.1. Due to the triple nature of the sys-tem, there is (i) a waved pattern on the outer edge of the gravity wake of M₃, this time particularly around periastron side, (ii) a broadening of the gravity wake, resulting in thicker arcs in the edge-on view, and (iii) a reduced radial velocity of the wake, making it travelling a shorter radial distance in the same num-ber of orbital periods. The aforementioned effects render it more challenging to discern the eccentric nature of a triple system when examining its wind morphology, in comparison to a binary system.

To study in more detail how the eccentricity-induced asym-metries arise in model v15m06e_i40, Fig. 13 shows the onset of the gravity wake of M₃ at four different orbital phases. Under-standing the wind-companion interaction throughout the orbit is very challenging due to the many interferences and interacting components. Ignoring the inner binary spiral shock of M₁–M₂, the changing orbital velocity of the outer companion M₃ (both in absolute value and direction) strongly affects the opening angle and width of the M₃ gravity wake, resulting in the asymmetries revealed as well in the outer binary model v15bin_om06e40 in Fig. 5. Further, the interaction of the companion M₃ with the spiral shocks created by the inner binary M₁–M₂ strongly deter-mines the shape of the resulting gravity wake of M₃. The time in-between two encounters of M₃ with a spiral shock becomes phase-dependent because of the varying orbital velocity and sep-aration along the orbit, and the varying distance between two M₁–M₂ spiral shocks. Further, the impact of an encounter of M₃ with a passing spiral shock depends on the density con-trast of the shock with the wind. This contrast is larger when the outer edge of the M₁–M₂ compact binary spiral has not yet caught up and merged with the previous inner edge. In this model, because of the varying orbital separation the M₃ compan-ion sometimes resides within and sometimes beyond the region where these edges merge. Finally, the shape of the M₁–M₂ binary spiral is affected by the phase-dependent orbital velocity of the ${\text{CoM}}_{M_1+M_2}$, and by the direct gravitational attraction of wind material towards companion M₃.

Taking into account these complexities, we attempt to understand the different shapes of the M₃ gravity wake at different orbital phases in Fig. 13 as follows. After apastron passage (first plot, t = 508 yr), the companion M₃ is accelerating, which results in a decrease in the width of its gravity wake. Because M₃ has a radially inward orbital velocity component, it encounters successive outward moving M₁–M₂ spiral shocks at a faster rate. Additionally, between apastron passage (x > 0, y = 0) and this orbital phase, the orbital radius of M3 is relatively large, such that it resides in a region where the inner and outer edge of the M₁–M₂ spiral already merged, and the density contrast of the passing spiral shocks with the wind is relatively low. As a consequence of these two effects, rather small waves result on the outer edge of the gravity wake, which fade out at larger radii (Fig. 12). In contrast, around periastron passage (second plot, t = 536 yr), M₃ resides in a region where the edges of the M₁–M₂ spiral have not yet merged, such that the density con-trast of the passing spiral shocks and the wind is larger, which results in more impact-full collisions. After periastron, the radi-ally outward directed orbital velocity reduces the rate at which the companion encounters spiral shocks. Further, it encounters strong spiral shocks, because it takes some time before it enters the region where the M₁–M₂ spiral edges have merged. This results in the large waves in the outer edge of the gravity wake that can be seen in the third plot of Fig. 13 (e.g. the wave at −120 au ≲ x ≲ 0 au and −100 au ≲ y ≲ −50 au at t = 575 yr). As the companion approaches apastron, its orbital velocity direc-tion becomes perpendicular to the radial velocity of the spiral shocks, which increases the rate at which spiral shocks pass by. Additionally, the encounters are less violent due to the larger orbital separation, such that the outer gravity wake edge becomes smoother again. This explains why the gravity wake in the global orbital plane morphology in Fig. 12 is more wavy at periastron side (x < 0), and smoother at apastron side (x > 0).

Fig. 11 Density distribution in a slice in the orbital plane at two different orbital phases for model v15m06ei40. The apastron passage of the inner companion star M2 is on the x > 0, y = 0, z = 0 axis.

Fig. 12 Density distribution in a slice in the orbital plane (top) and in a perpendicular meridional plane (bottom) for model v15m06eo40, with apastron passage of the outer companion at x > 0, y = 0, z = 0.

Fig. 13 Density distribution in a slice in the orbital plane at four different orbital phases within the outer orbital period of model v15m06eo40. The companion star M3 is at periastron passage in the second image (x < 0, y = 0) and at apastron passage in the fourth image (x > 0, y = 0, z = 0).

5 Comparison to R Aql

Observations of the circumstellar envelope of the AGB star R Aql reveal a very high degree of complexity, indicating that one or multiple, still undetected, companions are probably shaping the wind (Decin et al. 2020). This inspired us to examine the possibility that a triple stellar system is embedded and hidden within this outflow, and to explore which type of wind structures and global morphologies are expected to form around such triple systems. We therefore modelled the set of systems, described and analysed in Sects. 2 and 4, in such a way that the primary AGB star has similar properties as R Aql.

The exact nature of any companions of R Aql and the orbital properties are still unknown. Therefore, the parameter space that could be explored to model R Aql is extremely large and degenerate. However, in this section we perform a prelim-inary comparison of the observations of R Aql to synthetic observations of one or our triple models that lies within this degenerate parameter space. A thorough study of the observa-tional data, including different molecular lines, which all provide different insights on the possible nature of its central system and on the complex envelope at different spatial scales within the wind, is beyond the scope of this work and is needed before a more detailed comparison with a simulation can be made. To describe the observed morphology, we solely focus on the ¹²CO J = 2→1 rotational line emission, which is an optimal tracer of the density in the entire envelope, due to its large fractional abundance with respect to molecular hydrogen, and because it is predominantly collisionally excited (De Beck et al. 2010). We first shortly describe this observational data, and then examine what we can learn from a comparison to a synthetic observation of one or our models.

5.1 ALMA observation

R Aql is an O-rich AGB star that is located at a distance of 230 pc from Earth. Its mass is estimated to 1.6 M_⊙ and its mass-loss rate is $~\dot{M}=1.1\times {10}^{-6}{\text{M}}_{\odot }{\text{yr}}^{-1}$ (Decin et al. 2020). R Aql was observed as part of the ATOMIUM ALMA large programme, with three different ALMA antenna array configurations: compact, mid and extended (with baseline ranges of 15– 500, 15–1398, and 111–13894m, respectively) (Decin et al. 2020; Gottlieb et al. 2022). To optimally reveal both small-and larger-scale structures, a combination is made of the data from these three configurations (with a uv taper of FWHM 0.04”, which is equivalent to a baseline of ~9 km s⁻¹ at 1.3mm wavelength, with restoring beams (major axis [mas], minor axis [mas], angle [degrees]), of (40.7, 36.0, 22.38) and (58.3,49.3,25.12) for the continuum and CO cube, respectively. We refer to Gottlieb et al. (2022) for more details on the observations and data processing.

In the observational images, 1 arcsec corresponds to a spa-tial extent of 230 au. The continuum emission of the combined ATOMIUM ALMA dataset is shown in Fig. 14. It consists of one centrally condensed patch of emission with a specific intensity peak of 18.9 mJy beam⁻¹. The total flux density within the 3 × σ_rms contour with a radius of 0.05″, is 19.5 mJy beam⁻¹. Except for a slight elongation of the emission patch due to the beam shape, no deviations from spherical symmetry are resolved, and no close companions can be identified in this way. Outside the 3 × σ_rms contour, the emission is dominated by Gaussian noise with a mean close to zero, which is typical for an image obtained by interferometry.

Emission maps of the combined ¹²CO J = 2→1 emission, with blue- and red-shifted projected velocities ranging from −12.3 to 11.8 km s⁻¹ with respect to the stellar velocity of 47.2 km s⁻¹ are shown in Fig. 15. In most channels, the bulk of emission is confined within a roughly circular region cen-tred around the star, so there is no clear sign of rotation or deviation from a roughly spherical, radial velocity profile. The emission has a maximal radial extent of about 14″ (corresponding to 3200 au) around the central velocity, and decreases for increasing projected velocity. Around the maximal velocity of ≈ ± 9–10 km s⁻¹, at the red-shifted side (8.0 and 9.3 km s⁻¹ chan-nel maps) the emission is located more to the north, whereas at the blue-shifted side (−8.5 and −9.8 channel maps) there is an asymmetry towards the south. This indicates that the global morphology could be slightly oblate instead of spheri-cal, and inclined with the side towards us pointing slightly to the south, corresponding to a small positive inclination, i. Fur-ther, at smaller scales, there are many intricate inhomogeneities including non-circular concentric shells, arcs, and spiral-like structures. For example, the outermost emission structure, which extends up to 7″ at the central velocities, appears as a triangular morphology in the red-shifted channels v = 0.4 to v = 4.2 km s⁻¹, and it is rather circular in the blue-shifted channels. At smaller radii, there are many more intricate structures that are hard to disentangle and identify as shells, arcs, or spirals.

Fig. 14 Continuum emission of the R Aql system. Contours are drawn at 3,12,48, and 96 times the continuum rms noise value (1 × 10−4 Jy beam−1).

5.2 Comparison to synthetic observations

Now that we understand the general impact of a triple compan-ion on the resulting wind structures of mass-losing stars from Sect. 4, we compare these with the observations of R Aql. The typical features appearing in the simulations of triple systems, being an inner spiral (or small arcs) with an additional larger scale spiral (or larger scale arcs), and a waved pattern in the outer spiral, are not distinguishable in the observed ¹²CO J = 2 → 1 emission of R Aql in Fig. 15. However, it is not straightforward to compare the 2D slices of our 3D simulations (in Sect. 4) to these channel maps. Firstly, the channel maps are not slices in the x-y-z spatial dimensions, but 2D maps with material with a certain projected velocity towards the observer. To locate the material on these channel maps in the 3D space, one needs to know the entire wind velocity field, which is generally unknown and can be very complicated (see e.g. De Ceuster et al. 2024). Secondly, the den-sities in the simulations do not represent exactly the observed ¹²CO J = 2 → 1 emission. Therefore, we need to post-process the output of our models with a radiative transfer code, to calculate the level populations and flux densities, and in this way provide synthetic channel maps and line profiles of certain observed molecules. Additionally, to mimic the process of the interfer-ometric observations and the resulting noise and lost flux, an ALMA-simulator needs to be applied to these synthetic observ-ables. Thirdly, in the hydrodynamic simulations in Sect. 4, the inclination has a big impact on how the different structures and global morphologies appear. In observations, the inclination and orientation at which the system is viewed is generally unknown, since the orbit and location of the possible companion(s) are unconstrained, but indications for the inclination can be found from analysing specific observed structures or global asymme-tries. The observations of R Aql might have a slightly positive inclination and oblate global morphology, given the asymmetric location of the emission around the terminal projected velocities (see Sect. 5.1).

We estimate with large uncertainty that the average radial distance between two observed arc- and spiral-like structures in the outer part of the outflow morphology in the central velocity channels in Fig. 15 is ~1.15 arcsec, indicating that these wind structures travel about ~265 au during one orbital period. Assuming that the radial wind velocity in this region is ~10 km s⁻¹, the time between two structures (corresponding to the orbital period), is ~126 yr. Model v08m06 (with char-acteristics given in Tables 1 and 2) consists of an inner binary M₁ + M₂ = 2 M_⊙ and outer companion M₃ = 0.6 M_⊙ at a separation of 35 au and initial wind velocity below 10 km s⁻¹. Given that its outer orbital period is P_o = 128.5 yr, this is a model that could represent the R Aql system.

We therefore compare the morphology in the observation and simulation of model v08m06, to investigate if the outflow of R Aql might be shaped by a similar triple system. To do this, we post-process the PHANTOM output of model v08m06 with the MAGRITTE¹ radiative transfer code, which is a ray-tracer that solves the radiative transfer equation using a second-order Feautrier scheme along each ray (for more details see De Ceuster et al. 2020a,b, 2022). Since the hydro-simulations do not model the chemistry, we assume that all mass is in H₂, and the fractional abundance of CO is [CO/H₂]= 3 × 10⁻⁴, which is a reasonable assumption for an O-rich AGB star (Van de Sande et al. 2019). MAGRITTE uses the atomic and molecular database LAMDA (Schöier et al. 2005), and the calculations are done in non-local thermodynamical equilibrium (NLTE), taking into account the first 40 rotational transitions. The micro-turbulence is set to 0.5 km s⁻¹ (Ramstedt et al. 2008), and the cosmic microwave background radiation is taken as outer boundary condition.

Assuming the system is observed at the same distance as R Aql (230 pc) and under an inclination of 30°, the simulated ¹²CO J = 2 → 1 intensity at similar projected velocities as in the R Aql channel maps (Fig. 15) is shown in Fig. 16. The morphol-ogy consists of bright emission patches that move from north to south of the star between the projected velocities of 8 to 0 km s⁻¹ and 0 to −8 km s⁻¹. This corresponds to wind material residing close to the polar regions, and producing a two-peak pattern on the emission line. The bright emission patches lie on top of the underlying large-scale spiral pattern (gravity wake of outer companion M₃), which is deformed compared to the orbital plane structure in the model, due to the inclination and the uncertain deprojected 3D coordinates of the wind material in each channel map. The underlying close binary spiral pattern and the wave pattern in the outer spiral edge are more difficult to distinguish in these channel maps, which also ignore the noise simulations by ALMA. Figure B.1 shows in more detail the central channel map with v_proj = 0.56 km s⁻¹. The compact spiral pattern is visible within the central region (within ~1 arcsec), and the waved pattern of the outer edge of the M₃ spiral wake is distinguishable at small positive relative RA. We addition-ally performed a test in which we post-process a simulation of the outer binary of this triple system, with M_p = 2.0 M_⊙, M_s = 0.6 M_⊙, v_ini = 8 km s⁻¹, a = 35 au. Assuming the same inclination and distance as the triple model in Fig. 16, we found that with the binary simulation, a clear spiral pattern and less asymmetry appear in the channel maps, making the two cases distinguishable.

Comparing the synthetic observations of the triple model (Fig. 16) to the actual R Aql data (Fig. 15), the global morphol-ogy has a similar shape. It is quasi-spherical around the central velocity channels, and decreases in radial extent with increasing velocities. In both the actual and synthetic observations, at the maximal projected velocities, an emission patch is visible above and below the stars, which is due to the inclination and elongated morphology (see edge-on density profile of model v08m06 in Fig. 7). When looking at the smaller-scale structures, in both the observation and simulations, there are rather intricate spiral-like structures. A detailed one-to-one comparison is, however, not feasible, due to the large uncertainty on the selected parameters of the model and the inclination, and because the SPH models lack some physics such as the impact of the stellar pulsations, which may affect the wind morphology. One has to be very careful with interpreting any more detailed (non-global) resem-blances between the model and observation. Given the before mentioned various uncertainties, in this case, such resemblances may be rather coincidental and be created by entirely different physical processes. Further, a more direct comparison of the MAGRITTE synthetic observables to the observations requires to post-process the data with an ALMA-simulator, to model the noise. This is beyond the scope of this work, as it would further complicate the comparison. R Aql might be shaped by a simi-lar triple system, under a slight positive inclination, but further detailed studies on the observational data are required to discover the actual nature of this system.

Fig. 15 Channel maps of the combined12 CO J = 2→1 emission of R Aql, with blue- and red-shifted projected velocities (annotated in upper right corner in km s−1), ranging from −12.3 to 11.8 km s−1 with respect to the stellar velocity of 47.2 km s−1. Contours are drawn at 3,6 and 12 times the rms noise value of 6.6 ⋅ 10−4 Jy/beam.

Fig. 16 Synthetic channel maps of the 12CO J = 2→1 emission of model v08m06, post-processed with the radiative transfer code MAGRITTE (De Ceuster et al. 2020a,b, 2022), assuming an inclination of 30° and distance of 230 pc. A yellow star is added indicating the location of the central star.

6 Conclusion

With the use of high-resolution 3D SPH PHANTOM models, we investigated the effect of a co-planar hierarchical triple system on the morphology of AGB winds. In general, adding an outer companion, M₃, to an inner binary, M₁–M₂, results in an addi-tional, large-scale spiral structure, on top of the more compact spiral structure created by the inner binary. This additional large-scale spiral is a gravity wake behind the outer companion, and has a wavy outer edge, due to its interaction with the under-lying compact M₁–M₂ spiral structure. The added companion, M₃, also induces an orbital motion of the inner binary, M₁–M₂, around the system’s centre of mass. This makes the inner M₁–M₂ spiral deviate from an Archimedes spiral and resemble a snail-shell pattern. When decreasing the wind velocity, we find that, similar to the binary case, the global morphology significantly changes and becomes more flattened, and the impact of the companions on the morphology increases.

We also explored the impact of eccentric orbits on these morphologies. We find that the effect of introducing eccentric-ity in the inner orbit on the large-scale structures is limited. On smaller scales, this introduces asymmetries that are char-acteristic of eccentric systems, as was studied in Malfait et al. (2021). If the outer orbit is eccentric, the effect is larger, but still less pronounced then in the binary model, indicating that the triple nature reduces the effect of eccentricity. The signatures of eccentricity are clearest when comparing the large-scale edge-on structures at the apastron and periastron side of the companion’s orbit. In an eccentric triple system, very complicated, phase-dependent wind-companion interactions take place, resulting in asymmetric spiral and arc structures.

Finally, we included a brief study of the wind structures emerging in the complex outflow of the AGB star R Aql, observed within the ATOMIUM ALMA large programme. Due to its chaotic morphology, this star could host a triple system. Although a detailed one-to-one comparison with one of our simulations is unfeasible, we used the radiative transfer code MAGRITTE to post-process model v08m06, which could resemble the R Aql system observed at an inclination of ~30°. From the modelled and observed ¹²CO J = 2→1 emission, it is not straightforward to identify the triple nature of the system from the emission maps, because the characteristics identified in the triple models are not resolved in the channel maps. The morphology of R Aql could be the result of a similar triple system to the ones studied in this work, but no clear direct implications of this are detected, and a more detailed study of the observational data of R Aql is required.

A direct comparison between observations and simulations of AGB outflows is very challenging, due to the large uncer-tainties in the wind and orbital properties, and the uncertainties that come with comparing a 3D simulation to 2D projected velocity maps of a wind with an uncertain velocity field. Additionally, more relevant physical, chemical (Maes et al. 2024), and approximate radiation (Esseldeurs et al. 2023) processes, which may play a role in shaping the observed outflow, are still being added to the SPH wind model. Alternatively, in contrast to the bottom-up approach presented in this paper, model recon-struction methods might help in directly obtaining probable 3D models for our most complex observations (see e.g. De Ceuster et al. 2024).

This work provides an important first step in modelling the wind structures resulting from the interaction with companions beyond the binary paradigm.

Data availability

The input files and output dumps to reproduce the Phantom simulations and figures in this paper are available online (https://zenodo.org/records/12795466).

Acknowledgements

The computational resources and services used in this work were provided by the VSC (Flemish Supercomputer Center), funded by the Research Foundation Flanders (FWO) and the Flemish Government – depart-ment EWI. We thank Daniel Price for his help, especially concerning the development of PHANTOM. The animations and figures in this paper are made using the visualisation tools PLONS² and SPLASH (Price 2007). ALMA is a partnership of ESO (representing its member states), NSF (USA) and NINS (Japan), together with NRC (Canada) and NSC and ASIAA (Taiwan), in cooperation with the Republic of Chile. The Joint ALMA Observatory is operated by ESO, AUI/NRAO and NAOJ. JM, ME, and LD acknowledge support from the KU Leu-ven C1 excellence grant C16/17/007 MAESTRO and the FWO research project G099720N. L.S. is a senior F.N.R.S. research associate. FDC is a postdoctoral research fellow of the Research Foundation – Flanders (FWO), grant 1253223N. LD acknowledges support from the KU Leuven C1 excellence grant BRAVE C16/23/009 and KU Leuven Methusalem grant SOUL METH/24/012. SHJW acknowledges support from the Research Foundation Flanders (FWO) through grant 1285221N.

Appendix A Additional figures hydrodynamic simulations

This section contains additional figures of the wind morphology in binary and hierarchical triple systems.

Fig. A.1 Velocity distribution in a slice through the orbital plane of model v15binom06 for r < 475 au.

Fig. A.2 Density distribution in a slice through the orbital plane of models v08m06 (left) and v15m06 (right) for r < 90 au.

Appendix B Additional figure synthetic observation

This section contains an additional figure of the synthetic observation of model v08m06.

Fig. B.1 Synthetic channel map of the 12CO J = 2→1 emission with v = 0.56 km s−1 of model v08m06 presented in Fig. 16.

References

All Tables

All Figures

	Fig. 1 Set-up of hierarchical triple with inner binary M1 –M2 orbiting around ${\text{CoM}}_{M_1+M_2}$ and outer companion M3 orbiting around the system’s centre of mass ${\text{CoM}}_{M_1+M_2+M_3}$. M3 induces an orbital motion of ${\text{CoM}}_{M_1+M_2}$, and thereby also of M1 and M2, around the system’s centre of mass.
In the text

	Fig. 2 Density (top) and velocity (bottom) distribution in a slice through the orbital (left) and meridional (right) plane of model v08bini for r ≤ 475 au.
In the text

	Fig. 3 Same as Fig. 2, but for model v15bini.
In the text

	Fig. 4 Density distribution in a slice in (top) and perpendicular to (bottom) the orbital plane of model v15binom06.
In the text

	Fig. 5 Density distribution in a slice in (top) and perpendicular to (bottom) the orbital plane of model v15binom06e40. The apastron passage of the companion is on the y = 0, x > 0 axis
In the text

	Fig. 6 Density distribution in a slice in the orbital plane (top) and in a perpendicular meridional plane (bottom) for models v15m01 (left), v15m06 (middle), and v15m12 (right).
In the text

	Fig. 7 Same as Fig. 6, but for models v08m01 (left), v08m06 (middle), and v08m12 (right).
In the text

	Fig. 8 Density and velocity distribution in a slice in the orbital plane (top) and in a perpendicular meridional plane (bottom) for models v15m01 (left), v15m06 (middle), and v15m12 (right), for r ≤ 475 au.
In the text

	Fig. 9 Same as Fig. 8, but for models v08m01 (left), v08m06 (middle), and v08m12 (right).
In the text

	Fig. 10 Density distribution in a slice in the orbital plane (top) and in a perpendicular meridional plane (bottom) for model v15m06ei40, with apastron passage of the inner companion at x > 0, y = 0, z = 0.
In the text

	Fig. 11 Density distribution in a slice in the orbital plane at two different orbital phases for model v15m06ei40. The apastron passage of the inner companion star M2 is on the x > 0, y = 0, z = 0 axis.
In the text

	Fig. 12 Density distribution in a slice in the orbital plane (top) and in a perpendicular meridional plane (bottom) for model v15m06eo40, with apastron passage of the outer companion at x > 0, y = 0, z = 0.
In the text

	Fig. 13 Density distribution in a slice in the orbital plane at four different orbital phases within the outer orbital period of model v15m06eo40. The companion star M3 is at periastron passage in the second image (x < 0, y = 0) and at apastron passage in the fourth image (x > 0, y = 0, z = 0).
In the text

	Fig. 14 Continuum emission of the R Aql system. Contours are drawn at 3,12,48, and 96 times the continuum rms noise value (1 × 10−4 Jy beam−1).
In the text

	Fig. 15 Channel maps of the combined12 CO J = 2→1 emission of R Aql, with blue- and red-shifted projected velocities (annotated in upper right corner in km s−1), ranging from −12.3 to 11.8 km s−1 with respect to the stellar velocity of 47.2 km s−1. Contours are drawn at 3,6 and 12 times the rms noise value of 6.6 ⋅ 10−4 Jy/beam.
In the text

	Fig. 16 Synthetic channel maps of the 12CO J = 2→1 emission of model v08m06, post-processed with the radiative transfer code MAGRITTE (De Ceuster et al. 2020a,b, 2022), assuming an inclination of 30° and distance of 230 pc. A yellow star is added indicating the location of the central star.
In the text

	Fig. A.1 Velocity distribution in a slice through the orbital plane of model v15binom06 for r < 475 au.
In the text

	Fig. A.2 Density distribution in a slice through the orbital plane of models v08m06 (left) and v15m06 (right) for r < 90 au.
In the text

	Fig. B.1 Synthetic channel map of the 12CO J = 2→1 emission with v = 0.56 km s−1 of model v08m06 presented in Fig. 16.
In the text