Shaping planetary nebulae with jets in inclined triple stellar systems

TL;DR

This paper uses 3D hydrodynamical simulations to show that jets launched in inclined triple stellar systems can produce highly irregular, asymmetric planetary nebulae, supporting the idea that such systems shape messy PNe.

Contribution

It demonstrates, through simulations, that inclined jets in triple systems can create irregular planetary nebulae, a novel explanation for their complex shapes.

Findings

Jets in inclined triple systems produce asymmetric PNe.

Simulations show filamentary, irregular nebula structures.

Supports triple systems as a shaping mechanism for messy PNe.

Abstract

We conduct three-dimensional hydrodynamical simulations of two opposite jets launched obliquely to the orbital plane around an asymptotic giant branch (AGB) star and within its dense wind, and demonstrate the formation of a `messy' planetary nebula (PN), namely, a PN lacking any type of symmetry (highly irregular). In building the initial conditions we assume that a tight binary system orbits the AGB star, and that the orbital plane of the tight binary system is inclined to the orbital plane of the binary system and the AGB star (the triple system plane). We further assume that the accreted mass on to the tight binary system forms an accretion disk around one of the stars, and that the plane of the disk is tilted to the orbital plane of the triple system. The highly asymmetrical and filamentary structure that we obtain support the notion that messy PNe might be shaped by triple stellar systems.

Figures

Full text

binaries: close $-$ planetary nebulae $-$ jets

1 INTRODUCTION

The notion that binary stellar systems (e.g, Bond & Livio 1990; Bond 2000; De Marco 2015; Zijlstra 2015) or planetary systems (e.g, De Marco & Soker 2011) shape planetary nebulae (PNe) has gained a huge support over the years, both from theoretical considerations that point to the difficulties of models based on single asymptotic giant branch (AGB) stellar progenitors to shape PNe (e.g., Soker & Harpaz 1992; Nordhaus & Blackman 2006; García-Segura et al. 2014), and from a large number of observations and their analysis (e.g., from 2015 on, Akras et al. 2015; Aller et al. 2015a, b; Boffin 2015; Corradi et al. 2015; Decin et al. 2015; De Marco et al. 2015; Douchin et al. 2015; Fang et al. 2015; Gorlova et al. 2015; Hillwig et al. 2015; Jones 2015; Jones et al. 2015; Manick et al. 2015; Martínez González et al. 2015; Miszalski et al. 2015; Močnik et al. 2015; Montez et al. 2015; Jones et al. 2016; Chiotellis et al. 2016; Akras et al. 2016; García-Rojas et al. 2016; Jones 2016; Hillwig et al. 2016a, b; Bond et al. 2016; Chen et al. 2016; Madappatt et al. 2016; Ali et al. 2016; Hillwig et al. 2017). In most cases of shaping by binary interaction, jets are involved (e.g., Morris 1987; Soker 1990; Sahai & Trauger 1998; Boffin et al. 2012; Huarte-Espinosa et al. 2012; Balick et al. 2013; Miszalski et al. 2013; Tocknell et al. 2014; Huang et al. 2016; Sahai et al. 2016; Rechy-García et al. 2016; García-Segura et al. 2016, to list a small number out of hundreds of papers). A relatively small number of papers raise the possibility that triple stellar systems are behind the morphologies of some PNe (e.g., Soker et al. 1992; Soker 1994, 2016; Bear & Soker 2017; Jones 2017). A claim for a triple stellar system in the PN SuWt 2 (Bond et al., 2002; Exter et al., 2010) was rejected recently by Jones & Boffin (2017).

In general there are two orbital planes. One is that of the more tight binary system, and the second one is that of the triple stellar system, i.e., of the tight binary system motion around the center of mass with the third star. If the two orbital planes coincides, the interaction leads to a mass loss geometry that has a plane of symmetry, even if the mass loss geometry departs from axi-symmetry. Binary systems, e.g., if the interaction time scale is shorter than the orbital period or if the orbit is eccentric, can also cause departure from axisymmetry (e.g., Soker & Hadar 2002).

If the two orbital planes of the triple system are inclined to each other, the morphology of the descendant PN is likely to possess no symmetry plane. The descendant PN is actually likely to lack any symmetry; nor point-symmetry, nor axial-symmetry, and nor mirror symmetry. We term this a ‘messy PN’.

As part of the increasing interest in mass transfer in evolving triple stellar systems (e.g., Michaely & Perets 2014; Portegies Zwart & van den Heuvel 2016), we extend our study of the shaping of PNe by jets, now with properties as expected in some triple stellar systems. Specifically, we consider the flow structure described by Soker (2004), where a tight binary system orbits an AGB star and accretes mass from the AGB wind. Because the orbital plane of the tight binary system is not parallel to the orbital plane around the AGB star, the jets’ axis is not perpendicular to the orbital plane of the triple system. We describe the initial setting and the numerical code in section 2. We then describe the numerical results in section 3. We summarize our results in section 4.

2 NUMERICAL SET-UP

2.1 Initial conditions

Soker (2004) studied triple stellar systems where a tight binary system orbits an AGB star that blows a dense wind. The two orbital planes, that of the tight binary system and that of the triple stellar system are inclined to each other. The tight binary system accretes mass from the wind. Due to the orbital motion of the two stars of the tight binary system around their mutual center of mass, the accreted mass possess specific angular momentum that leads to the formation of an accretion disk around one (or two) of the stars. The plane of the accretion disk is close to being parallel to the orbital plane of the tight binary system, rather than to that of the triple stellar system (Soker, 2004). The accretion disk launches two opposite jets that are almost perpendicular to the orbital plane of the tight binary system.

We do not study the mass accretion process, but assume that two opposite jets are launched from an origin that orbits the AGB star. We do note though, that the accretion rate might be larger by a factor of several relative to the classical Bondi-Hoyle-Lyttleton accretion rate, as the mass transfer might take place via the so called wind Roche-lobe overflow (RLOF; Mohamed & Podsiadlowski 2011). In a full wind RLOF Mohamed & Podsiadlowski (2011) claim that the accretion rates are at least an order of magnitude greater than the classical Bondi-Hoyle-Lyttleton accretion rates. In Fig. 1 we present a schematic drawing of the wind blown by the AGB star and the jets launched by the tight binary system, and the orbit of the jets’ origin around the AGB star.

We simulate two cases differ in the inclination angle $\phi$ of the jets’ axis to the line perpendicular to the orbital plane of the triple system. In the first case $\phi=20^{\circ}$ and in the second case $\phi=30^{\circ}$. The other parameters are identical for the two runs and are as follows. The orbital separation of the jets’ origin (the location of the tight binary system) to the center of the AGB star is $a_{123}=30~{}\rm{AU}$. The orbital period is $P_{123}=67~{}$yrs, corresponding to a total mass of the triple stellar system of $M_{123}=M_{1}+M_{2}+M_{3}=5M_{\odot}$.

At the beginning of each simulation the slow AGB winds fills the grid, and it is continuously replenished along the simulation by injecting mass from the numerical surface of the AGB star at a radius of $r=1145R_{\odot}$ from the center. At this radius the wind almost reaches its terminal velocity. The mass loss rate and velocity into the slow AGB wind are $\dot{M}_{w}=5\times 10^{-5}M_{\odot}~{}\rm{yr}^{-1}$ and $v_{w}=10~{}\rm{km}~{}\rm{s}^{-1}$, respectively.

At $t=0$ we start to launch two jets, along the $+z$ and $-z$ directions. The jets are conical, with a half opening angle of $\alpha_{j}=30^{\circ}$. The opening angle is similar to the values we have been using over the years ($\simeq 30-50^{\circ}$) while studying the shaping of PNe by jets (e.g., Akashi & Soker 2013). The jets opening angle influences the morphology of the nebula. However, we have limited numerical resources and cannot study more cases at this time. The initial velocity of the jets is $v_{j}=500~{}\rm{km}~{}\rm{s}^{-1}$, that is about equal to the escape speed from a main sequence star. The mass loss rate into the two jets together is $\dot{M}_{2j}=2\times 10^{-7}M_{\odot}=0.004\dot{M}_{w}$. This corresponds to about $10\%$ of the Bondi-Hoyle-Lyttleton accretion rate from the wind.

2.2 The numerical code

Our simulations are performed by using version 4.0-beta of the FLASH code (Fryxell et al., 2000). The FLASH code is an adaptive-mesh refinement (AMR) modular code used for solving hydrodynamics or magnetohydrodynamics problems. Here we use the unsplit PPM (piecewise-parabolic method) solver of FLASH. We neither include gravity, as velocities are much above the escape speed in the region we simulate, nor radiative cooling. Radiative cooling is important. At this first study we did not include radiation from the mass-accreting star(s) as well. As we will see below, the gas does suffer a substantial adiabatic cooling. Radiative cooling will further reduce the gas temperature, while radiation from the stars can heat it. As we did not include radiative cooling and radiation, our results cannot reproduce exact details of observed PNe, but non the less give the overall morphological characteristics of ‘messy nebulae’.

We employ a full 3D AMR (8 levels; $2^{11}$ cells in each direction, or 10 levels; $2^{13}$) using a Cartesian grid $(x,y,z)$ with outflow boundary conditions at all boundary surfaces. The center of the grid is at $(0,0,0)$. The AGB star is placed at $(0,0,-a_{123}\sin\phi$). The tight binary system orbits the AGB star along a circle. The projection of the orbit on the $y=0$ plane is shown in Fig. 1. We simulate the whole space with no symmetry assumptions. The jets are continuously launched during the entire simulated time.

The slow wind and the jets start with a temperature of $1000~{}\rm{K}$. The initial jets’ temperature has no influence on the results (as long as it is highly supersonic) because the jets rapidly cool due to adiabatic expansion. For numerical reasons a weak slow wind is injected from the jets’ origin in the sector $\alpha_{j}<\theta<90^{\circ}$. This wind has a small effect on the outcome, as it contains very little energy. Without this wind numerical instabilities might develop on the boundary between the ambient gas and the sphere from where the jets are launched (more numerical details are in Akashi & Soker 2013 ).

3 RESULTS

In Fig. 2 we present density maps in the $y=0$ plane for the $\phi=20^{\circ}$ case (see Fig. 1), at four times. In all panels the AGB star is located just below the center of the grid at $(x,y,z)=(0,0,-10)~{}\rm{AU}$. The source of the jets orbits the AGB star, and its location is $(x,y,z)=(-21.5,-19.4,-18.1)~{}\rm{AU}$ at $t=41~{}\rm{yr}$, $(x,y,z)=(25.1,13.6,-1.1)~{}\rm{AU}$ at $t=72~{}\rm{yr}$, $(x,y,z)=(-19.6,-21.5,-17.4)~{}\rm{AU}$ at $t=109~{}\rm{yr}$, and $(x,y,z)=(26.2,-11,-0.7)~{}\rm{AU}$ at $t=130~{}\rm{yr}$.

The main outcome of the simulation is a nebula that lacks any symmetry. We simulate only the early interaction time, and cannot deduce the exact final nebular shape that is observed at later time. Nonetheless, our results clearly demonstrates the early formation of a ‘messy PN’. The arrows depict the flow direction at each point, but not its value. We can notice a complicated flow structure in the low density regions that are enclosed by the dense shell. These low density regions are filled with the shocked jets’ material. These regions contain little mass, and will not be observed in the visible when the object turns into a PN. The complicated flow structure in these regions, including many vortexes, forms many filaments in the dense gas; these will be observed in the PN phase. We will now turn to describe only the $\phi=30^{\circ}$ case.

In Fig. 3 we present the density maps for the $\phi=30^{\circ}$ simulation. Now the AGB star is located below the center of the grid at $(x,y,z)=(0,0,-15)~{}\rm{AU}$, and the source of the jets that orbits the AGB star is located at $(x,y,z)=(-19.8,-19.4,-26.4)~{}\rm{AU}$ at $t=41~{}\rm{yr}$, $(x,y,z)=(23.2,13.6,-1.6)~{}\rm{AU}$ at $t=72~{}\rm{yr}$, $(x,y,z)=(-18.2,-21.5,-25.5)~{}\rm{AU}$ at $t=109~{}\rm{yr}$, and $(x,y,z)=(24.2,-11,-1)~{}\rm{AU}$ at $t=130~{}\rm{yr}$. The results are not much different from those for the $\phi=20^{\circ}$ case presented in Fig. 2. The two main morphological features are a ’messy’ shape and the formation of a filamentary structure.

We present the temperature maps (in Kelvin) in Fig. 4, corresponding to the same planes and times as in Fig. 3. The jets are shocked to temperatures of few millions degrees Kelvin. The post-shock gas suffers a substantial adiabatic cooling. Nonetheless, the high temperatures regions in the maps correspond to post-shock jets’ material. These are regions of low density, that when surrounded by dense AGB gas are termed bubbles.

In Fig. 5 we present the density maps (again, for the $\phi=30^{\circ}$ case) at $t=130~{}~{}\rm{yr}$ and in four different planes perpendicular to the $x-y$ plane. The angle $\theta$ is the angle between the plane and the $x-z$ plane shown in Fig. 3. In these panels the horizontal axis is not the $x$ axis of the grid that is used in the previous figures, but it is actually $r=(x^{2}+y^{2})^{1/2}$. In Fig. 6 we plot the velocity magnitude, $v=(v_{x}^{2}+v_{y}^{2}+v_{z}^{2})^{1/2}$, given in units of $~{}\rm{cm}~{}\rm{s}^{-1}$ by the color bar, and show the velocity direction with arrows, in the same four planes. The four panels in these two figures emphasize the lack of any symmetry in the outflow, i.e., a ‘messy nebula’, as each panel presents a different density structure. The panels show also the presence of many dense filaments in the nebula.

In Figs. 7 and 8 we present the density structure and temperature maps with velocity direction, respectively, in six planes parallel to the $z=0$ plane, characterized by their distance from the $z=0$ plane, all at $t=130~{}\rm{yr}$. As before, these panels emphasize the messy nebula and the presence of filaments.

In Fig. 9, we present the three-dimensional density structure of the gas at $t=130~{}\rm{yr}$. Each panel emphasizes surfaces with different densities as indicated. These panels fully reveal the messy nature of the nebula and the filaments in the outflow.

A discussion on the limitations of these simulations is in place here. Ideally we should have run the simulations until the AGB envelope is removed, and then blow the fast wind from the central remnant of the AGB star (the hot core to become a WD). This is impossible due to severe resource limitations when we conduct 3D hydrodynamical simulations. Despite this limitation, for two reasons we argue that our simulations, that show that a messy circumstellar matter is formed, can be projected to the later PN phase. (1) The messy nature of the circumstellar matter that appears after one orbit (panel b in Figs. 2 and 3), does not disappear after 2 orbits (panel d in these figures), and the density contrasts between different segments of the circumstellar matter even increase. So we can project and conclude that running the simulations for over hundreds orbits (thousands of year) till the AGB envelope disappear will not erase the messy nature of the nebula. (2) At the end of our simulations the outer shell expands at a velocity of $v\simeq 50~{}\rm{km}~{}\rm{s}^{-1}$, and its sound speed is $C_{s}\simeq 20~{}\rm{km}~{}\rm{s}^{-1}$, i.e., the Mach number is $\mathcal{M}\simeq v/C_{s}\simeq 2.5$. This is high enough to ensure that the messy structure will not be smoothed out. We can summarize this point by stating that although we cannot compare this results directly to images of messy PNe, we can confidently claim that the messy circumstellar we have obtained will develop into a messy PN.

4 SUMMARY

We conducted 3D hydrodynamical simulations of fast jets interacting with a slow wind from an AGB star. The source of the jets orbits the AGB star with an orbital period of $P_{123}=67~{}$yrs. The significant new ingredient in the simulations was that the symmetry axis of the jets is inclined to the angular momentum axis of the triple-stellar system (see Fig. 1). Such a flow might result if the jets are launched by a star that is a member of tight binary system that orbits the AGB star, accretes mass from its wind, and the orbital plane of the tight binary system is inclined to the orbital plane of the triple system (Soker, 2004).

The evolution of the density structure in one plane is presented in Figs. 2 and 3 for two different cases simulated. The density structures in different planes at the end of the simulation for the $\phi=30^{\circ}$ case are presented in Figs. 5 and 7. These figures clearly show that the nebula lacks any symmetry, and has many filaments. This is termed a ‘messy nebula’. The high temperature regions that are seen in Figs. 4 and 8 are regions of post-shock jets’ material. They have low density, and are termed bubbles. The complicated flow structures of the post-shock jets’ material, including vortexes, seen in Figs. 2, 6, and 8, account in part for the filamentary structure of the dense gas. The full glory of the messy nebula is seen in Fig. 9.

One outcome of the motion of the jets in an inclined orbit is that each jet encounters a denser AGB wind than the other jet during half of the orbital period. Over all, the jets do not penetrate in our simulations to large distances, and the bubbles that are formed are small. For example, an ‘hourglass’ nebula is not formed in our simulations. One possible example (Soker, 2004) might be the PN NGC 6210 (images by Balick 1987; Pottasch et al. 2009), that might as well ended in a merger (Soker, 2016).

As stated earlier, we cannot follow the evolution to the PN phase, and hence one to one correspondence with observed messy PNe is impossible. However, we can learn from the recent study of Bear & Soker (2017). They present images of 4 typical messy PNe, that they claim are most likely have been shaped by a triple stellar system. The two PNe H 2-1 (PN G350.9+04.4; original image from Sahai et al. 2011) and Hen 3-1333 (PN G332.9-09.9 or CPD-568032; original image from Chesneau et al. 2006) are relevant to our simulations. These two messy PNe lack large lobes and do not have any pair of two opposite lobes. The PN H 2-1 has a bright broken thin shell (a rim), while the PN Hen 3-1333 has no such rim. We expect that the PNe that result from the flow structure we have simulated will look somewhat like one of these messy PNe. Whether a bright thin shell is formed or not, depends on whether the fast wind during the pre-PN phase is strong or not.

We end by crudely estimating the fraction of PNe shaped by the studied scenario. Bear & Soker (2017) estimate that about $13$ to $21$ per cents of non-spherical PNe have been shaped by triple stellar systems. Soker (2016) estimates that about one in eight non-spherical PNe is shaped by the PN AGB progenitor swallowing a binary system or interacting with a binary system close to its surface. These are subset of all triple stellar scenarios. The scenario studied in this paper, therefore, accounts for some fraction of the 12 per cents of PNe studied by Bear & Soker (2017). We crudely estimate that the scenario studied in this paper accounts for about 3-5 per cents of all non-spherical PNe. There are several scenarios of triple stellar shaping, each accounts for several percents of all non-spherical PNe. The surviving binary system around the central star of the PN is very hard to detect because the binary system is far from the center, and the combined luminosity of the two stars is much smaller than that of the central star. In any case, our results strengthen the call (Bear & Soker, 2017) for the search of a triple stellar systems at the center of some messy PNe.

This research was supported by the Israel Science Foundation. N.S. is supported by the Charles Wolfson Academic Chair.