Pinwheel Outflow Induced by Stellar Mass Loss in a Coplanar Triple System

We develop a physical framework for interpreting complex circumstellar patterns whorled around asymptotic giant branch (AGB) stars by investigating stable, coplanar triple systems using hydrodynamic and particle simulations. The introduction of a close tertiary body causes an additional periodic variation in the orbital velocity and trajectory of the AGB star. As a result, the circumstellar outflow builds a fine non-Archimedean spiral pattern superimposed upon the Archimedean spiral produced by the outer binary alone. This fine spiral can be approximated by off-centered circular rings that become tangent to each other at the location of the Archimedean spiral. The superimposed fine pattern fades out relatively quickly as a function of distance from the center of the system, in contrast to the dominant Archimedean spiral pattern, which presents a much slower fractional density decrease with radius. The different rates of radial decrease of the density contrast in the two superimposed patterns, coupled with their different time and spatial scales, lead to an apparent, but illusory radial change in the observed pattern interval, as has been reported, for example, in CW Leo. The function describing the detailed radial dependence of the expansion velocity is different in the two patterns, which may be used to distinguish them. The shape of the circumstellar whorled pattern is further explored as a function of the orbital eccentricity and the inner companion’s mass. Although this study is confined to stable, coplanar triple systems, the results are likely applicable to moderately noncoplanar systems and open interesting avenues for studying noncoplanar systems.

1. INTRODUCTION Planetary nebulae (PNe) and preplanetary nebulae (pPNe) exhibit a vast variety of morphologies, and spherical PNe are extremely rare.Only 3.4% of 119 young PNe and none of 23 pPNe in the survey using the Hubble Space Telescope have a round shape, while the majority (> 60%) of them were classified as bipolar or multipolar (Zuckerman & Aller 1986;Sahai et al. 2007Sahai et al. , 2011)).The nonspherical morphologies, in general, include highly collimated jet-like features, bipolar compact knots, point-symmetric filaments (or strings of knots), bipolar ansae, and equatorially dense disk/toruslike features.Also invoked nowadays are objects dis-Corresponding author: Hyosun Kim hkim@kasi.re.kr playing moderate deviations from spherical symmetry, such as spirals, intertwining off-centered circular rings, and irregularly distributed arcs, often found to be surrounding mass-losing asymptotic giant branch (AGB) stars and their compact descendants at the nuclei of pPNe and PNe.As noted in a review by Balick & Frank (2002) and references therein, no single mechanism for shaping offers a comprehensive explanation of all observational properties of such systems, but it is nowadays commonly accepted that these structures require an interaction with a "close" binary companion (De Marco 2009).
A spiral-shell pattern surrounding a mass-losing giant star is an important tool for identifying binarity (e.g., Soker 1994;Mastrodemos & Morris 1999;Kim & Taam 2012b).But the observed whorled patterns examined to date mostly, if not all, indicate "wide" (longperiod) binaries.By searching ∼ 650 optical and in-frared images of pPNe and PNe, Ramos-Larios et al. (2016) found that 29 sources possess whorled patterns, and the time lapse between their consecutive rings and arcs ranges from 500-1200 yr.The time intervals between successive waves in the circumstellar whorled patterns are also estimated using spatio-kinematic structures revealed in molecular line observations, including 350, 300, and 800 yr for the well-known carbon-rich AGB stars (carbon stars), R Scl, CIT 6, and AFGL 3068, respectively (Maercker et al. 2012;Kim et al. 2013Kim et al. , 2017)).The (sub)millimeter images obtained in a new survey of oxygen-rich AGB stars seem to indicate the time interval to be not less than 30 yr (Decin et al. 2020).Therefore, so far, the time interval of the circumstellar whorled pattern (corresponding to the orbital period of the binary) does not seem to indicate the presence of companions that are close enough to induce the morphological transitions from AGB to PN stages.This observational result is possibly due to a bias by observers who tend to target close, large sources by taking into account the limitations of telescopes in angular resolution and sensitivity.Another possibility, which we focus on in this paper, is that observers may have missed the spatial clues implying the presence of another (inner) companion in a triple system appearing as a minor pattern having less surface brightness (or density) contrast and a smaller pattern interval.
Classical statistical studies of stellar populations indicated that the fraction of triple to binary systems is about 0.11 and the fraction for higher multiplicity systems consisting of n objects is f n = N n /N n−1 ∼ 0.25 (Duquennoy & Mayor 1991;Tokovinin 2001).Thanks to high-resolution imaging, recent discoveries of additional subsystems in known binaries further boost the fraction of hierarchical multiplicity relative to binaries up to 30-50% (Raghavan et al. 2010;Hirsch et al. 2021).This result implies that at least 3 out of 10 observed patterns whorled around evolved stars may present additional features induced by a tertiary component.In addition, close binaries tend to be accompanied by at least one other star, with the probability of finding additional companions increased by up to 96% for spectroscopic binaries in the shortest-period group in Tokovinin et al. (2006).This naturally suggests that the progenitors of strongly bipolar PNe, which presumably originated from close binaries, may also possess an outer whorled pattern characterized by a longer time interval.
Before invoking a third stellar component, a binary model with a highly eccentric orbit can be considered to explain a bipolar outflow in terms of the short, but impactful, pericenter passages of the companion star.As an example, the carbon star V Hya has experienced high-speed bullet-like ejections once every ∼ 8.5 yr, which are hypothesized to be produced during the periapse passages of a binary companion having an orbital period of ∼ 8.5 yr (Sahai et al. 2016).However, Salas et al. (2019) raised an issue with the long-term orbital sta-bility of such a binary system because the eccentric orbit should be circularized within a relatively short time by the dynamical and tidal interactions of the companion with the AGB star.They proposed that a triple system could be a solution for V Hya by invoking an orbital configuration in which the inner companion (≲ 0.01 M ⊙ ) grazes the Roche limit of the mass-losing star in an eccentric orbit, while the eccentricity of that orbit can be maintained by the gravitational influence of an outer companion.
The multiple-shell structure of the closest carbon star, CW Leo (a.k.a.IRC+10216), has been under debate since the discovery of the shells (Mauron & Huggins 1999).The multiple shells appear roughly spherical, but they are not simply explained by equally spaced concentric circles; some of them are incomplete arcs, some are spiral-like, and some even appear to intertwine with other shells.Many authors also noticed that the centers of curvature are offset from each other (e.g., Mauron & Huggins 1999;Guélin et al. 2018, see also Guelin et al. 1993).By combining molecular line data taken with the Atacama Large Millimeter/submillimeter Array (ALMA; at a resolution of ∼ 0. ′′ 3) with those from the Submillimeter Array (SMA; at a resolution of ∼ 3 ′′ ), Guélin et al. (2018) found a regular interval of shells of ∼ 16 ′′ (or 2000 au at the distance of ∼ 123 pc derived by Groenewegen et al. 2012) in the outer envelope (up to a radius of ∼ 110 ′′ ), which is consistent with a binary model having an orbital period of ∼ 700 yr.Notice that Guélin et al. claimed that the orbit of the binary star system is eccentric and is viewed nearly face-on (see also Cernicharo et al. 2015).This geometry associated with a face-on and eccentric orbit binary was supported by a separate study of the position-angle dependence of the transverse wind velocity (Kim et al. 2021;Kim 2023).However, this proposed geometry is contradictory to the nearly edge-on geometry proposed based on the elongated shape of the dust continuum emission and the bipolar-like optical image of the central 1 ′′ region (e.g., Men'shchikov et al. 2001;Jeffers et al. 2014;Decin et al. 2015).Guélin et al. also noted that the spatial and time intervals between shells decrease in the inner (and younger) part of the envelope at radii r < 10 ′′ , where ∆r ∼ 2 ′′ (∆t < 100 yr), compared to larger radii between 10 ′′ and 40 ′′ , where ∆r = 5 ′′ -10 ′′ (∆t ∼ 300 yr).Their speculation that such a reduction of orbital period was caused by mass transfer either to the envelope or to the companion star fails because, as they argue, it requires a much higher mass-loss rate than the current rate of CW Leo.These discrepancies between the inner and outer parts of the circumstellar envelope may imply the presence of an additional (inner) companion.
In this paper, we present hydrodynamic and particle (kinematic) simulations for the circumstellar structure induced by a triple system consisting of a mass-losing star and two companion stars without mass loss.We mostly focus on particle simulations in order to explicitly track the circumstellar density-enhanced pattern and its velocity structure without the complexities due to gas pressure effects.Comparison of the particle simulation with the corresponding hydrodynamic simulation emphasizes the purely hydrodynamic effects, in particular, the density dispersion caused by shocks and gas pressure.Unlike hydrodynamic simulations, particle simulations can employ a constant velocity for the intrinsic wind of the mass-losing star, which facilitates the comparison of models over a wide parameter space for the stellar properties, e.g., the stellar masses.In a hydrodynamic simulation, even in a single star case, the wind velocities throughout the simulation domain are dependent upon the stellar masses (see Section 2.2).Particle simulations also have an advantage in computational efficiency, with which we can trace the pattern out to the radii where the pattern induced by the third star has been repeated on multiple scales.This paper will also inform the merits and limits of particle simulations, which often have been, and will be, utilized for quick modeling of observed whorled patterns in the interpretational framework of binary (and hereafter, multibody) systems.
In Section 2 we describe our numerical methods for determining the orbits of a triple system (Section 2.1), for the hydrodynamic simulations (Section 2.2), and for the pinwheel model based on following the trajectories of wind particles ejected at different moments (Section 2.3).The results are shown in Section 3. The implications of our models are discussed, and our conclusions are presented in Section 4.

Coplanar Triple System
The preservation of dynamical stability within triple systems hinges upon the implementation of a hierarchical configuration, wherein an inner binary is orbited by an outer entity with a much wider orbital trajectory (Salas et al. 2019).Systems with smaller ratios of outer to inner orbital periods are susceptible to being unstable because of the eccentric Kozai-Lidov effect (Kozai 1962;Lidov 1962), unless the system is coplanar and the orbits are near-circular.Among evolved star systems that are suspected to have a third object, the AGB star π 1 Gru could be an example of a hierarchical triple with an orbital period ratio of > 500, suggested by the recent finding of evidence for a close companion with an orbital period of only 10 yr (Homan et al. 2020).However, owing to the very large ratio between the orbital periods, the circumstellar patterns associated with the individual companions could have a large density contrast with respect to each other.It would likely leave a vestige of one of the companion stars on the circumstellar envelope that would be too tenuous to be observable, either because of the density attenuation of the larger pattern or because one companion is too low mass to build a significant pattern.As our aim in this work is to find and track the third object's footprints remaining in the circumstellar medium, as may have been observed in some AGB sources, like CW Leo, a large orbital period ratio is not addressed in our investigations.
Adopting such a framework, we assume a stable triple system, for which we further assume coplanar orbits.Coplanar orbits are relatively less affected by the eccentric Kozai-Lidov mechanism, albeit not completely free of such a mechanism, leading to large-amplitude eccentricity and inclination oscillations in near-coplanar triple systems (Li et al. 2014).The eccentric Kozai-Lidov mechanism is not taken into account in this paper because the timescales for its operation are much longer than the few-orbit timescales that we are considering here, and possibly even longer than the duration of the AGB phase.
Seven parameters that determine the orbits within a triple system include the masses of the three objects (M A , M B , and M C ), the average separation between the center of mass of the inner (A-C) binary system and the outer B object (a AC + a B ), the average separation between the inner objects (a A + a C ), and their eccentricities (e AC−B and e A−C ).Here, a represents the semimajor axis of the orbit of a star relative to the center of mass of the corresponding (AC-B or A-C) binary system.
We first calculate the orbits of a binary system composed of two mass components, M A + M C and M B , separated by a AC + a B , which determines the final orbit of object B as well as the time-dependent position of the center of mass of the inner binary system, objects A and C. The mass ratio, M A /M C , and the A-C separation, a A + a C , then impose the individual orbits of objects A and C with respect to the time-dependent position of the center of mass of the inner binary system.
Figure 1 illustrates the orbits of a triple system (a) in XY coordinates corotating with object B about the center of mass of the three objects (coinciding with the coordinate origin) and (b) in the inertial frame of observers sitting on the +z-axis.In the former frame, the locations of the center of mass of the inner binary system and object B are fixed, and the orbits of A and C are closed and circular (elliptical, in eccentric orbit cases).In the observer's frame, however, the orbits of A and C are not necessarily closed.The small fluctuation in the orbit of the mass-losing star A (red curve in Figure 1(b)), relative to a circle, is the primary cause of the main features described in this paper.
This study is limited to triple systems having coplanar orbits in order to capture the essential features created by the presence of a third object.Finally, the orbital motions are assumed to be stable in their predefined circular or eccentric orbits, with constant orbital parameters.In future studies, stability considerations will be included in the modeling of the observations.

Hydrodynamic Simulations
Three-dimensional hydrodynamic simulations for triple systems were performed using version 4.5 of the code FLASH (Fryxell et al. 2000) by solving the Eulerian hydrodynamic equations in Cartesian coordinates with the origin at the center of mass of the triple system.The adaptive mesh refinement scheme imposes a higher level of refinement toward the instantaneous stellar positions.The refinement is determined by the stellar positions, not by the usual density or velocity criteria.Specifically, the maximum refinement level of 9, with 64 grid points per block length for a total image size of 6000 au, imposes the highest resolution of about 0.37 au near the stellar positions; it assigns 19 grid points to the diameter (7 au) of the hypothesized wind-launching sphere surrounding the mass-losing star, where the initial vector velocities in the observer's frame of the wind material blowing out of the star are reset every simulation time step.The mass-losing star is treated as a point source having a gravitational softening radius (also called the Plummer radius) of 1 au.The mass-loss rate for all models presented in this paper has been set at 3 × 10 −6 M ⊙ yr −1 ; the morphological results are insensitive to the assumed mass-loss rate.The equation of state of an ideal gas is assumed with the ratio of specific heats, γ, chosen to be equal to 1.4.Radiative cooling and heating are not included for simplicity; we only follow the density of the material.
Following a commonly adopted technique for wind acceleration (e.g.Theuns & Jorissen 1993), the gravitational force attributed to the mass-losing star of mass M A located at ⃗ r = ⃗ r A is reduced by a factor of 1 − f , where the wind acceleration factor f is a constant representing the ratio of the outward force due to radiation pressure on dust grains to the inward gravitational force (see the wind velocity profiles for different constant values of f in Kim & Taam 2012a for an isothermal gas; the wind solution for a polytropic gas can be found in Shivamoggi et al. 2021).In a test simulation for a corresponding single mass-losing star located at the center of the simulation domain, the intrinsic wind profile as a function of radius followed a supersonic branch of the above papers as the updated version of Parker's wind solution (Parker 1958), and the terminal velocity was around 13 km s −1 .The intrinsic wind velocity is governed by the momentum equation of hydrodynamics; therefore, it would be scaled down with an increased companion mass.
The gravitational forces attributed to the outer and inner companion stars of mass M B and M C located at ⃗ r = ⃗ r B and ⃗ r = ⃗ r C , respectively, are also implemented in the code.In order to explore the effects of the perturbed orbital path of the mass-losing star and to distinguish those effects from those created by the density wakes formed behind the companion stars, we run simulations that alternatively include and exclude the terms for the companions' gravitational forces by setting the shutoff parameter, W, defined by Kim et al. (2019) and Kim (2023), to 1 and 0, respectively.In any case, the orbit of the mass-losing star is preset according to the dynamics of the triple system, as described in Section 2.1.The results of the hydrodynamic simulations are displayed in Figures 2(a) and (b) for density and in Figures 3(a) and (b) for the expansion velocity; a further description can be found in Section 3.

Particle Simulations
We consider that most of the characteristics of the spiral-shell pattern induced by a binary or triple system can be explained by a simpler model that tracks the trajectories of wind particles ejected from the mass-losing giant star.In order to identify the purely hydrodynamic effects, we compare the density and velocity fields obtained through a hydrodynamic simulation with the corresponding particle model characteristics calculated by the pinwheel code, as described below.
The first version of the pinwheel model accumulates the particles that are freely moving with the instantaneous momenta gained at the moments of their ejection from the mass-losing star as it passes through the predefined positions along its orbit with the orbital velocities as calculated in Section 2.1.Here, the particle velocity is defined as the vector sum of the intrinsically isotropic wind velocity (denoted by V exp ) with the orbital velocity of the star.If N par particles are ejected with equal spacing over all solid angles, from one of the discrete N pos positions along the orbit of the star within one orbital period, then within the observing time t = N turn orbits, a total of N par × N pos × N turn particles are moving within the computational domain.Among them, the number of particles located within each grid cell is used as a measure of the density distribution.The velocity map is also calculated by averaging the velocity vectors of individual particles falling into each grid cell.An example binary model that was calculated by this method is shown in Kim et al. (2017)'s supplementary Figure 4.
The second version of the pinwheel code is written using the concept of a piston based at the systemic center of mass.Once the direction of a piston (θ, ϕ) is chosen, with the usual notations θ for the polar angle and ϕ for the azimuthal angle, the time-variable length and velocity of the piston head are determined by the projection of the position and velocity of the star in its orbit onto the piston direction: where (x A , y A ) and (v x,A , v y,A ) indicate the position and velocity of object A in the orbital plane.The wind particle is ejected with speed V exp + V piston toward the piston direction.The output images of the first and second versions of the pinwheel model are apparently the same under the condition of a far-distance approximation (i.e., r ≫ l piston ).In the third version of the pinwheel code, the particles are made to be sticky.For this, a piston approximation is used for convenience.Whenever the particles ejected at later times overtake previously ejected particles, the velocities of all these particles at the same distance are averaged and updated to a common new value.This assumes that the particles spatially coinciding with each other perfectly stick together.A binary model of this version was exhibited in He (2007).
The fourth version of the code reduces the efficiency of stickiness so that the velocity dispersion of the particles that are coincident with each other is reduced to the predefined speed of sound; if the initial velocity dispersion of the overlapping particles is smaller than the specified speed of sound, their velocities are not updated.
Figures 2(c) and 3(c) present the resulting density and velocity maps of the first (or second) version of the pinwheel model (i.e., without stickiness), while Figures 2(d) and 3(d) show the third version's results for particles with the efficiency of stickiness set to 100%.The example density maps produced via the fourth version of the piston model with the velocity dispersion of 0.5 km s −1 for the coincident particles are exhibited in Sections 3.3 and 3.4.Table 1 summarizes the parameters for the particle simulations displayed in this paper.We note that a binary system consisting of two objects with masses M A + M C and M B induces the densityenhanced structure of an Archimedean spiral in the orbital plane, as denoted by the black line in the figure, satisfying where T AC−B represents the orbital period in this binary system.Figure 2(a) shows that the triple system creates the same spiral structure as in the abovementioned binary system (hereafter referred to as the main spiral), along with finer structures in the inter-ridge region of the main spiral.As shown in the right panel for the r-ϕ polar coordinate map, the finer structures have a spiral shape with the opposite orientation (i.e., the slope in the polar coordinate map is negative).The segments of the finer spiral structure effectively merge into the main spiral ridge, rather than passing continuously across it.The overall features in Figure 2(a) are compared to the correspondences in Figure 2(b), where the gravitational wakes of objects B and C are excluded; we have further compared it to the images (not shown in this paper) made by excluding just one of the gravitational wakes of objects B and C. The nonlinear perturbations appearing along the main spiral, as clearly seen in both the left and right panels in Figure 2(a), can be attributed to the overlay of the wake of the outer companion (object B).The irregularities in the main spiral have a vertically limited extent, as indicated in previous studies (see, e.g., Kim & Taam 2012a,b;Kim et al. 2019).On the other hand, the broadening of the fine spiral segments shown in Figure 2(a), relative to the width of the corresponding features shown in Figure 2(b), occurs with the introduction of the gravitational effect of the inner companion (object C).The density profile across the fine spiral is also modified from two peaks shown in Figure 2(a) to one peak in Figure 2(b).We also note that the slope of the black line in the right panel is greater in Figure 2(b) than in Figure 2(a) because of the absence of the companions' gravitational influences in Figure 2(b), which yields an expansion velocity that is about 0.5 km s −1 greater than in Figure 2(a).The deceleration of the wind is smaller when the companions' gravity is excluded.
With increasing distance from the system's center of gravity, the relative density decrease of the finer pattern is much faster than that of the main spiral.Furthermore, while the radial broadening of the main spiral is insignificant, the width of the finer pattern increases with radius, eventually smoothly filling the region between the density ridges of the main spiral.Therefore, the finer pattern in the outermost part would likely be unidentified in, for example, molecular line observations of the circumstellar medium of an AGB star in a triple system.
In Figure 3 for the distribution of expansion velocity, the velocity gradient is maximized at the location of the main spiral pattern, traced by the black solid line.Overall, the expansion velocity gradually increases from one arm to the next outer arm of the main spiral.The velocity jumps at the locations of fine spiral structures are minor, which provides an important characteristic to distinguish the main and fine spirals of a triple system among the observed intensity peaks.

Particle Simulations: Differentiating Hydrodynamic and Nonhydrodynamic Effects
Panels (c) and (d) of Figure 2 present the results of pinwheel model calculations with stickiness efficiencies of 0% and 100%, respectively.Both images well approximate the locations of the whorled patterns observed in the hydrodynamic model displayed in Figure 2(b).The width of the main spiral pattern in the hydrodynamic simulation is, however, smaller than in the nonsticky model shown in Figure 2(c) and larger than in the extremely sticky model in Figure 2(d).It suggests the necessity of adjusting for the efficiency of the stickiness of the particles when they coincide.
The net velocity distribution, as plotted in Figure 3, is also overall well reproduced by the pinwheel model, in particular with the sticky particles.In the perfectly sticky model shown in Figure 3(d), the expansion velocity sharply drops from > 14 to < 11 km s −1 at the main spiral's ridge, whose width is unresolved in this model.The corresponding hydrodynamic model, drawn in Figure 3(b), has a similar velocity distribution, except for the slightly larger width of the ridge of the main spiral.In the nonsticky model, the relatively smaller velocity between the split edges of the main spiral is distributed over a wider area (see Figure 3(c)), inconsistent with the hydrodynamic result shown in Figure 3(b).Within the region between the split edges of the main spiral, it is also found that the continuation of the sharp velocity structures of the finer pattern having the opposite orientation does not cross the black solid line (see the right panel of Figure 3 Although there is a close similarity in the global distribution of the expansion velocity of fluid, the number of fine spirals in the inter-ridge region of the main spiral is doubled in the hydrodynamic model (Figure 3(b)) compared to that in the pinwheel model with sticky particles (Figure 3(d)).The fine spirals have an extremely small width in the sticky model, while they show broadening in both the nonsticky and hydrodynamic models, accompanied by double-peaked velocity profiles (see Figure 6 of Kim et al. 2019, and below for more details).
In the full three-body hydrodynamic model, the nonlinearities in the density wake of the outer companion are apparent in Figure 3(a).The fine spirals in the full three-body model are broad (compare Figures 2(a) and (b)), within which triple peaks are presented, owing to the overlay of the gravitational wake of the inner companion.Accordingly, the velocity structures of the fine partial spirals are more complex (compare Figures 3(a Figure 4 compares the density and velocity profiles between the 100% sticky particle model and the hydrodynamic model induced by the orbital motion of the mass-losing object.The red-dotted vertical lines indicate the positions of the peaks in the sticky particle model.Note that the typical radial structure of one spiral ridge is characterized by a one-peak density profile, surrounded by double-peaked temperature profiles, and the velocity profile with one inflection point (see Figure 6 of Kim et al. 2019).Each of the radial zones indicated by a horizontal two-headed arrow and shading in Figure 4 shows that the dispersed spiral ridge of the hydrodynamic model coincides with a common shape of the velocity profile.Furthermore, the density peaks in the sticky particle model coincide with inflection points in the velocity profile, as indicated by the red-dotted vertical lines.The darker regions represent the overlaps of these shaded regions.These shaded regions, determined based on the velocity profile, agree well with the enhanced regions in the density profile.
The spirals of the hydrodynamic model are enclosed by two shocks at their inner and outer edges.Their radial variations in density and velocity are very similar to the characteristic density and velocity profiles of an outgoing forward shock and a reverse shock, the latter of which accelerates material inward, relative to the expanding forward shock, as in a supernovae remnant (e.g., Figure 1 in Truelove & McKee 1999), albeit showing smaller shock jumps due to a much slower wind speed.The individual peak structures of the fine spiral, shaded in Figure 4, have such a velocity profile, including one inflection point, superimposed upon the larger-scale velocity variation following the main spiral pattern.As a consequence, the density peaks in the hydrodynamic model are largely dispersed through the regions enclosed by the forward and reverse shocks at the inner and outer edges of the spiral patterns.

Efficiency of the Stickiness of Wind Particles
The efficiency of stickiness is adjusted in the piston model by allowing a certain range of velocity dispersion of the coincident particles up to a predefined constant value.As a result, the width of the main spiral in the hydrodynamic model is reasonably reproduced by applying a velocity dispersion of ∼ 0.5 km s −1 to the pinwheel model (see Figure 5(a)), which corresponds to the adiabatic speed of sound of a gaseous medium at a temperature of ∼ 20 K.In our hydrodynamic simulations, the temperature shows an overall decrease with radius except for the jumps at the major and minor spiral patterns, and the interarm temperature is below 20 K beyond 1 kilo-au.

Eccentric Orbits
In this section, we demonstrate the influences of eccentricities of stellar orbits in the density distribution of the circumstellar spiral patterns.The velocity dispersion up to ∼ 0.5 km s −1 for the coincident particles is adopted.Figure 5(a) presents the case in which the orbital shapes of the three stars are all circular, as described in the previous sections.Compared to this, in the model with an eccentric orbit for the inner companion (e A−C > 0, Figure 5(b)) 1 , the individual patterns are similar in location but become very widened, in particular toward the apocenter of the mass-losing star (ϕ ∼ 0, which we have defined to be in the −x direction).
In the model with an eccentric orbit for the outer companion (e AC−B > 0, Figure 5(c))2 , a one-sided dearth of matter in the inter-ridge regions is clearly seen toward the pericenter of the mass-losing star (ϕ = π, or in the +x direction), just like in a typical eccentric binary system (Kim et al. 2015(Kim et al. , 2019)).The fine spirals accordingly congregate near the main spiral at the position angle corresponding to the pericenter of the mass-losing star, making a bowtie shape in the r-ϕ map.The knot of the bowtie occurs at values of ϕ corresponding to the pericenter of the mass-losing star, and the knot is tighter for larger orbital eccentricities of the outer companion.In this model, the broadening of the fine spiral pattern is nearly independent of the position angle.

Mass of the Inner Companion
The effect of the mass ratio M A /M C at a given M A + M C is also explored.We compared our fiducial model for the inner companion as a low-mass star (M C = 0.1 M ⊙ , Figure 5(a)) with the model for a slightly more massive star companion (M C = 0.3 M ⊙ , Figure 6(a)) and the model for a planetary-mass companion (M C = 0.01 M ⊙ , Figure 6(b)).In the pinwheel model at a constant wind velocity of 13 km s −1 , the fixed M A + M C guarantees the same slope of the main spiral pattern.The mass M A is 1.0, 0.8, and 1.09 M ⊙ , respectively, in Figure 5(a), Figure 6(a), and Figure 6(b).
From the comparison of the three cases for the mass of the inner companion, we find that the peak density of the pattern does not differ much, but the density in the inter-ridge regions is significantly reduced as the mass of the inner companion is increased, thereby increasing the density contrast of the pattern.The fine pattern formed due to the presence of a more massive inner companion is somewhat broader in width, which yields a countereffect in the density contrast because the column density in the finer spiral features becomes more spread out.The latter effect, governed by the velocity dispersion of the pattern, is, however, largely limited by the local speed of sound, which would decrease with increasing distance from the star.Because the speed of sound in the outer circumstellar envelope is likely to be lower than the single value used in the current simple calculations, the latter (widening) effect would be reduced in the outer circumstellar envelope, thereby preserving the high-density contrast.Therefore, the additional pattern established by a relatively massive inner companion tends to have a higher-density contrast, on top of the predominant binary-induced pattern.

Circular Ring Approximation
The fine spirals are formed from the introduction of the inner companion (object C) into the binary system composed of the mass-losing star (object A) and the outer companion (object B).The actual cause of these spirals is the wiggling of the orbital trajectory of the mass-losing star (the red curve in Figure 1(b)) and the consequent variation of its orbital velocity.The wiggling of the orbit around the center of mass of the objects A and C is related to their orbital period, T A−C .
In Figure 7(a), the circular rings illustrated in red very well approximate the fine spirals.These circular rings (or the components of the fine spiral) that has the time interval corresponding to the inner orbital period of T A−C ∼ 85 yr converge on and form the main spiral pattern, which has a time interval corresponding to the outer orbital period of T AC−B ∼ 690 yr.There-fore, the seemingly position-angle-dependent change in pattern interval is an illusion caused by the offsets of the centers of curvature.Also, note that the centers of the rings follow a spiral that can be expressed by Equation (3) but for T A−C instead of T AC−B (see Figure 7(b)).The radii of the rings increase by a constant value of 260 au during each of the inner orbits.
On the other hand, as shown in Figures 5(c) and (d), in the models with the outer companion moving along a highly eccentric orbit, the fine spirals can be approximated by rings centered along the x-axis, which is the line connecting the pericenter and apocenter of the orbit.
4. SUMMARY AND DISCUSSION Three-body systems are being recognized as important for explaining complex circumstellar structures and as a possible driver for the observed morphologies of a subset of objects undergoing the transition between the late stellar evolutionary phases from AGB to PN.In this paper, we have investigated the hydrodynamic and kinematic influences on the morphology of the expanding circumstellar envelope of a mass-losing star when it is being orbited not only by a relatively distant companion, but also by a third star that is closely orbiting the mass-losing star.We propose a triple system for the AGB star, CW Leo, which accounts for its complex circumstellar shell pattern, such as the off-centering of the ring-like pattern and the abrupt radial change in the interval of this pattern.
We have first demonstrated the circumstellar pattern induced by a triple system by performing a hydrodynamic simulation in which the outflowing gas ejected by the AGB star is affected only by its initial velocity and the mass of the AGB star.We then compared the resulting density and velocity fields of the gas to the full three-body hydrodynamic simulation in which the gas response to the mass of all three bodies was computed, revealing incidental substructures originating from the gravitational wakes of companions.
The hydrodynamical effects are elucidated through comparison with the pinwheel model, which is nonhydrodynamic as simply following particles ejected isotropically from the mass-losing star having a predefined orbital motion.The pinwheel model mimics the circumstellar spiral pattern of the hydrodynamic model, exactly coinciding in location.However, it does not reproduce the density and width of the pattern, thereby highlighting the hydrodynamical effect upon the fluid as it encounters the shocks at the inner and outer edges of the dense ridge of the spiral pattern.
As a result, the density and density contrast of the fine spiral pattern (newly revealed to occur in a triple system) quickly decrease as the radius increases, while those along the main spiral pattern decline more slowly, just as in a two-body system.Therefore, the fine pattern may be observable only in the inner portions of the circumstellar envelope, making the observed image appear to have undergone an abrupt change in the pattern interval, as has been claimed for CW Leo (e.g., Guélin et al. 2018).
Furthermore, the fine spiral pattern can be very well represented by off-centered circular rings having a regular radius increment.The center positions of the rings are located along a spiral, with the slope of that spiral in polar coordinates being related to the orbital period of the inner binary.The approximating off-centered rings coincide tangentially and therefore stack up along the ridge of the main spiral.
The detailed structure of the expansion velocity of the outflowing envelope is closely tied to the shape of the main spiral; after crossing the main spiral outward, the velocity quickly drops and slowly recovers its maximum value before reaching the next winding of the main spiral.Superimposed on that large-scale velocity pattern, the velocity jumps at the locations of the fine spiral pattern are much smaller.If the propagation speed of the matter is measurable, the detailed variations in the expansion velocity could contribute to the discrimination between the main and fine spirals in a triple system.
The dependence of the pattern on orbital eccentricities and inner companion mass has also been explored.A larger orbital eccentricity of the outer companion creates a density pattern that tightens a "bowtie" with the "knot" at the position angle toward the pericenter of the mass-losing star.And the density outside the bowtie significantly drops, causing a one-sided gap in the interpattern density.A larger orbital eccentricity of the inner companion broadens the width of the fine spiral, in particular at angles corresponding to the apocenter of the mass-losing star.Finally, a more massive inner companion reduces the density in the regions between the ridges of both the main and fine spiral patterns, thereby raising the density contrast of the pattern.
According to the results of our model calculations, CW Leo's abrupt radial change in the interval of the pattern and the off-centering of the pattern with intertwining ridges may hint at its nature as a triple system.As CW Leo has been getting brighter for already two decades, possibly indicating its evolutionary phase at the critical AGB-pPN transition moment (Kim et al. 2021), continuous monitoring of it would be beneficial to understanding the roles of the companions, in particular of the inner companion, in the morphological transformation.
In the case of another example, the AGB star π 1 Gru, recent ALMA observations reveal an HCN spiral having an expansion time interval between successive turns of ∼ 10 yr, which is consistent with the orbital period of a potential companion located at the secondary continuum peak, assuming that it has about one solar mass (Homan et al. 2020); the implied companion differs from the known G0V-type companion located 10 times further from the mass-losing star (Feast 1953;Ake & John-son 1992), plausibly suggesting a triple system with an orbital period ratio of > 500.However, the CO spiral found by Homan et al. (2020) has an interval of ∼ 80 yr, which is inconsistent with either of the above orbital timescales.It remains unclear whether this implies a fourth object in this system, or whether this is a triple system having a large eccentricity and/or noncoplanar orbits.
Only the midplane structures are examined in this paper.The three-dimensional structure of the pattern induced by a triple stellar system and the position-velocity diagnostics for interpretation of spectro-imaging observations are deferred to a future study.In this paper, the orbits of three stellar objects are assumed to be coplanar for simplicity.It would be interesting to study the complexity of the spiral-shell patterns by relaxing this orbital configuration.Dynamical stability consideration for the orbits will also need to be pursued in the future, especially for the eccentric orbit cases, since only a subset of triple systems will be stable on timescales long enough to produce a mass-losing AGB star.In a different subset of triple systems, the radius of the AGB star will ultimately swell to a size comparable to the radius of the inner orbit, creating a common-envelope binary, thereby leaving the original triple system in a well-separated binary configuration.
all models, the orbital periods are set the same as TAC−B = 690 yr and TA−C = 85 yr, yielding a ratio of ∼ 8.The gravitational wakes of the companion objects, B and C, are not considered in particle simulations, corresponding to the shutoff parameter W = 0 in the hydrodynamic simulations.

Figure 2
Figure 2(a) shows the density distribution in the (common) orbital plane of the triple star system with the mass-losing giant star of mass M A = 1 M ⊙ , the outer companion of mass M B = 1 M ⊙ , and the inner companion of mass M C = 0.1 M ⊙ .The distances from the giant star are 100 au and 20 au for the outer and in- ) and 3(b)).

Figure 5
Figure 5(d) shows eccentric orbits for both inner and outer companions, in which the fine spirals are threaded within the bowtie (as in Figure 5(c)), and the individual fine spirals are widened around ϕ ∼ 0 (as in Figure 5(b)).

Figure 1 .
Figure 1.Orbits of three objects in a triple system used for the calculations displayed in Figures 2 and 3.All three stars are initially aligned along the x-axis in the order of A-C-B from left to right, starting from their apocenters in eccentric orbit cases.(a) Orbits of objects A, B, and C with individual masses of MA, MB, and MC in the corotating XY frame in which the locations of object B and the center of mass of the A-C binary system are fixed.The center of mass of the triple system is located at the origin.The term aAC denotes the (average) distance of the center of mass of the A-C binary system from the center of mass of the whole system, while aB indicates that of object B. Dashed curves schematically indicate the directions of the motions of the center of mass of the A-C binary system (red) and of object B (blue) in the observer's frame.The notations aA and aC demonstrate the (average) distances of objects A and C, respectively, with respect to the center of mass of the A-C binary system.Red and green circles present the orbits of A and C in this rest frame.Red and green solid straight lines display the velocity vectors of objects A and C at their initial positions, and are scaled to each other.(b) Orbits of objects A (red), B (blue), and C (green) in the frame of nonrotating observers located on the +z-axis.

Figure 2 .
Figure 2. (a) Hydrodynamic simulation for the circumstellar matter distribution governed by orbital motions of three objects (with orbital parameters of MA = 1.0 M⊙, MB = 1.0 M⊙, MC = 0.1 M⊙, aAC + aB = 100 au, aA + aC = 20 au, eAC−B = 0, and eA−C = 0).The density in the orbital plane is displayed in a logarithmic scale in units of g cm −3 .The black solid line represents an Archimedean spiral with a pattern speed of 12.5 km s −1 .The angle ϕ is measured from the −x-axis in the clockwise direction.(b) Same as (a), but with the gravitational density wakes of objects B and C excluded by turning off their gravitational effects on circumstellar gas by setting the shutoff parameter, W, to 0. The black solid line corresponds to an Archimedean spiral with a pattern speed of 13 km s −1 .(c) A pinwheel model with a constant wind velocity of Vexp = 13 km s −1 , adopted to match the major spiral pattern with that in panel (b), and (d) the corresponding piston model with sticky particles.The color bar in (c) and (d) panels shows the number counts of particles within each grid cell, which is a function of the total number of particles present in the computational domain after the computation has reached a steady state.

Figure 3 .
Figure 3. Same as Figure 2, but illustrating the expansion velocity of the fluid in the orbital plane.

Figure 4 .
Figure 4. Comparison between hydrodynamic and sticky particle models in density and velocity profiles.The density profiles along the azimuthal angle ϕ = π/4, passing (x, y) = (−1, 1), in Figures 2(b) and (d) are drawn in black and red curves, respectively, in the top panel.The corresponding velocity profiles from Figures 3(b) and (d) are plotted in the bottom panel.Red-dotted vertical lines mark the peak positions of the sticky particle model.Shaded regions, along with the horizontal twoheaded arrows, identify the characteristic spiral ridges (see Figure 6 of Kim et al. 2019, and the relevant text therein).

Figure 5 .
Figure 5. Midplane density in pinwheel models of a three-body system in which the mass-losing star is object A. The stickiness of particles is adjusted to provide a speed of sound of 0.5 km s −1 .The employed parameters are Vexp = 13 km s −1 , MA = 1.0 M⊙, MB = 1.0 M⊙, MC = 0.1 M⊙, aAC + aB = 100 au, and aA + aC = 20 au.The orbital eccentricity of object B (eAC−B) is 0 in the upper two models ((a) and (b)) and 0.8 in the bottom two models ((c) and (d)).The orbital eccentricity of objects A and C (eA−C) with respect to their center of mass is 0 in (a) and (c) but 0.8 in (b) and (d).

Figure 6 .
Figure 6.Same as in Figure 5(a), but for (a) MA = 0.8 M⊙ and MC = 0.3 M⊙, while in (b) MA = 1.09M⊙ and MC = 0.01 M⊙.Notice that MA + MC is equivalent to the value for the models in Figure5; therefore, the orbital motions of object B and of the center of mass of objects A and C remain the same as in Figure5(a).

Figure 7 .
Figure 7. (a) Ring approximation for the locations of the fine spirals (red), overlaid upon Figure 2(d), the sticky particle model, in grayscale.(b) The positions of the centers of the rings illustrated in panel (a) with respect to the systemic center of mass.The increment of ring radius in each of the inner orbits is 0.26 kilo-au.

Table 1 .
Parameters for particle simulations