Shaping of Planetary Nebulae by Exoplanets

: (1) Background: We investigate the hypothesis that exoplanet engulfment can help explain the observed non-spherical planetary nebula population, as a complementary shaping mechanism to the binary hypothesis. The aim is to investigate the extent to which massive planets can explain the population of non-spherical planetary nebulae; (2) Methods: This research utilises a new tool to calculate the planet-fraction of planetary nebulae progenitor stars called SIMSPLASH ; (3) Results: we conclude that ∼ 15–30% of non-spherical planetary nebulae around single stars will have a history in which they engulfed a massive planet on the AGB; and (4) Conclusions: Engulfment of massive exoplanets may contribute signiﬁcantly to the formation of non-spherical planetary nebulae around single stars, yet appears to be insufﬁcient to explain them all.


Introduction
According to the generalised interacting stellar winds (ISW) theory, to form a non-spherical PN, there needs to be a mechanism which can produce a density enhancement in the equatorial region of a mass-losing AGB star, or a mechanism by which the mass-loss is collimated as the fast wind blows. Fortunately, the formation of non-spherical PNe can be explained by a combination of stellar rotation, magnetic fields, and binaries. It is clear that binarity plays the most major role in facilitating the shaping mechanisms. However, there are discrepancies between the stellar binary progenitor population and the population of non-spherical PNe [1]. A new hypothesis has emerged which says that a binary companion does not necessarily need to be stellar in nature and that the definition of a binary interaction in the shaping of PNe can be extended to include brown dwarfs and, as is the subject of this work, massive planets ( [2][3][4]).
Planets orbiting within a few AU of their parent star are vulnerable to engulfment when the star evolves to the Red Giant Branch (RGB) or Asymptotic Giant Branch (AGB) ( [4][5][6]). Throughout this work, 'engulfment' is taken to be the point at which the planet makes contact with the star. The planets may then transfer significant gravitational potential energy and angular momentum as they are destroyed within the star, and thus break the spherical symmetry of the resulting nebula. Changing tidal forces over time mean that a simple comparison of maximum stellar radius versus planetary orbit is not adequate to determine whether or not a planet will be engulfed.
We model two sets of exoplanet systems (i) the currently known massive planets, and (ii) a larger synthetic planet population (found in the works of Ida and Lin [7]) introducing a new research tool. By investigating both observational and numerical samples, we aim to get a robust general conclusion. In addition, relevant to the planet shaping hypothesis, is the idea that planet engulfment can effectively enhance mass-loss on the RGB, causing it to lose most of its envelope [8]. The result of such mass-loss on the RGB is that the star cannot evolve to the AGB which in turn prevents the formation of the planetary nebula.

Materials and Methods
SIMSPLASH (SIMulationS for the PLAnet Shaping Hypothesis) [9] is a new tool developed by Dr L. Boyle (author) to combine the theory, data, and techniques to explore the process of planet engulfment and ultimately to calculate the planet-fraction of planetary nebulae. All of the results presented here were obtained using SIMSPLASH which contains three user-selected modes of operation: Orbital mode, Single-epoch population mode, and Multi-epoch population mode.
When run in Orbital mode, SIMSPLASH calculates the orbital evolution of a planet as the host star evolves. It employs the tidal dissipation theory [10][11][12] and MIST (MESA Isochrones and Stellar Tracks) stellar models of Choi et al. [13], and determines whether a planet will be engulfed or not, for a given set of initial conditions. The Single-epoch population mode of SIMSPLASH is an efficient tool for analysing the known exoplanets. It can process samples from the Exoplanet Orbit Database, as well as synthetic planet populations derived from planet formation and evolution theory, with regard to planet engulfment in the progenitor population of planetary nebulae. While there is certainly a bias in the statistics for known exoplanets, population synthesis models such as these will further inform our understanding through observational developments.

Fate of Known Exoplanets
The Exoplanet Orbit Database was accessed with SIMSPLASH on 25 June 2017 and contained a total of 5287 planets orbiting 4345 stars. After the appropriate cuts were made to the data, the sample consisted of 300 star-planet systems and a further 56 planets in multiple planetary systems. Each planet flagged as being in a multiple planet system was first evolved with the Orbital mode. The MIST stellar models of [13] are adopted for this analysis and the known exoplanet population is evolved firstly using the measured metallicity corresponding to each host star. Secondly, in order to determine the effect of evolving the entire population using a single metallicity for each host star (for example the mean metallicity of the sample), the known exoplanets are also evolved with stellar models of the five different metallicities available in SIMSPLASH. Each star-planet system in the exoplanet orbit database was first evolved with the nearest of these initial metallicities to the measured metallicity of the host star. The remaining planet sample was evolved with the Orbital mode of SIMSPLASH. After taking the multiple planet systems into account according to their engulfment, the final sample consisted of 327 star-planet systems (300 data points with a single planet, and a further 27 which were in multiple planet systems), where 262 planets were found to be in the RGB engulfment zone of their host, 37 in the AGB engulfment zone, and 28 were at a safe semi-major axis. Figure 1 shows the number of planets in the respective engulfment zones as a function of host star mass when evolved with an initial metallicity according to the (rounded) measured host star metallicity. Note that there is a large inconsistency of sample sizes for the various mass bins, which explains the gaps in the data. It appears that higher mass stars are more likely to host a massive planet within the AGB engulfment zone than lower mass host stars. This is unsurprising as it is shown in Boyle [9] that the RGB engulfment zone is much larger than the AGB engulfment zone for a 1 M star, in contrast to the significantly larger AGB engulfment zone of a 1.5 M star. However, there is of course an observational bias inherent in the radial velocity method-close-in massive planets around solar mass stars are the easiest to find, and so the sample is dominated by planets which will be engulfed during the RGB. However, it is still worthwhile to investigate the probability of a star engulfing a known massive planet, since these results can be readily updated in future as the sample of exoplanets becomes more complete in terms of larger semi-major axes and host star masses. The planet-fraction (the fraction of PNe whose morphologies were influenced by the engulfment of a planet during the AGB phase), the binary-fraction (the fraction of PNe whose morphologies were influenced by a binary companion), and the fraction of stars which will evolve to a spherical PN were then obtained, and are summarised in Table 1. Table 1. The resulting PN progenitor population proportions when adopting the probabilities of engulfment according to the known exoplanets evolved with stellar models according to the measured metallicity of each host star. The planet and binary fractions of PN progenitors are denoted as f planet and f binary , respectively, while the fraction of PN progenitors which will evolve as single stars to form spherical PNe (i.e., single stars without a planet engulfment, and binary systems which do not interact) is denoted as f other . The predicted fraction of non-spherical PNe, f NS is 22%, of which 14% of the non-spherical population are expected to be the result of a planet engulfment ( f NS(Pl) ). According to the distribution of currently known exoplanets then, the fraction of PN progenitor stars which will engulf a planet on the AGB is only around 3%. However, if planet engulfment and binary interactions are responsible for all non-spherical planetary nebulae, this means that 14% of all non-spherical PNe evolving from a population similar to the current PN progenitor population of the Milky Way would be shaped by planets. This is a non-negligible proportion, particularly considering that 80% of the observed population of PNe are non-spherical. Accordingly, this gives 14% of this 80% population of PNe which we interpret as 11% which would have been shaped by planet engulfment. This leaves the remaining 69% being shaped by binary interactions. This still does not account for the discrepancy between the observed non-spherical and spherical fractions and those calculated to evolve via single-star evolution. However, as stated previously, the current sample of the known exoplanets is incomplete for intermediate mass stars and biased towards close-in planets, so this analysis will need to be repeated in the future as the sample of the known exoplanets is updated. Fortunately, as the sample of the known exoplanets becomes more complete, repeating the above experiment in the future will be a relatively efficient process now due to the development of SIMSPLASH.

Fate of Simulated Planet Population
The synthetic planet data described here consisted initially of 3000 disk systems for each of the star mass populations. After the appropriate cuts were made to the data, the number of planet star systems per mass was reduced significantly as was the total number of systems. In order to demonstrate the effect of the adopted host star metallicity for each star-planet system over an entire population, the number of engulfments were calculated again, but this time by evolving each star-planet system in the population with each of the five different metallicities available in SIMSPLASH. Thus, five population models were calculated.
To determine the probability of a star having and engulfing a planet, the occurrence rates of massive planets as a function of host star mass and metallicity must be determined. The massive planet occurrence rate of Gaidos and Mann [14] is adopted and is given by the term f GP . The number of engulfments in the population resulting from the orbital evolution models, for each of the five different stellar metallicity models, were fit with a logistic regression to determine the probabilities of stars engulfing a planet on either the RGB or AGB as a function of the host star mass. The probabilities from the five regression models were then interpolated with respect to metallicity to compute a finer grid of engulfment probabilities as a function of both mass and metallicity. The resulting grids were then combined with the occurrence rates of massive planets of [14] to determine the probability of having and engulfing planets, for all metallicities between 0.005 and 0.04. The probabilities of having and engulfing planets on the RGB and AGB, as a function of mass and metallicity, are plotted in Figure 2. The Kroupa [15] IMF is employed in this chapter to determine the number of stars in each mass bin, which are calculated in SIMSPLASH.
The synthetic models used in this section were kindly provided by Ida, and were generated with the Ida & Lin planet formation and evolution models. Figure 2 reveals the overall range of probabilities for RGB and AGB engulfments, respectively. This can be explained by again considering the distribution of engulfments in Figure 1 and also the [14] occurrence rates of massive planets. The majority of RGB engulfments in the sample lay within the 0.8 and 1.2 M host star mass bins in the histogram. The occurrence rates of massive planet of stars in this mass range for all metallicities are very low, but particularly at low metallicities where the occurrence rate is less than 10%. Since the occurrence rates of massive planets for low mass host stars and metallicities never exceeds 50-60%, the probability of a star having and engulfing a planet on the RGB remains below 60%. On the other hand, since the majority of AGB engulfments are by stars with masses 1.2 M , and the occurrence rates of massive planets in this host star mass range is much higher, the range of probabilities of a star having and engulfing a planet on the AGB reflects this.
The probability of having and engulfing a planet for a given mass varies largely depending on the mean metallicity of stars in a population, suggesting that high metallicity environments are more conducive for planet shaping of planetary nebulae. It has long since been known that there is a metallicity correlation for stars hosting massive planets [16,17]. Table 2 shows the planet-fraction ranges between 1-13% for the low and high metallicity models, respectively. Furthermore, because the numbers of RGB engulfments increase dramatically for low-mass stars with increasing metallicity, the relative numbers of PN progenitors in the population decreases. While this may seem counter intuitive, this can be explained by again considering the distribution of engulfments in Figure 2 and also the [14] occurrence rates of massive planets. The majority of RGB engulfments in the sample lay within the 0.8 and 1.2 M host star mass bins in the histogram. The occurrence rates of massive planets of stars in this mass range for all metallicities are very low, but particularly at low metallicities where the occurrence rate is less than 10%. Since the occurrence rates of massive planets for low mass host stars and metallicities never exceeds 50-60%, the probability of a star having and engulfing a planet on the RGB remains below 60%. On the other hand, since the majority of AGB engulfments are by stars with masses 1.2 M , and the occurrence rates of massive planets in this host star mass range is much higher, the range of probabilities of a star having and engulfing a planet on the AGB reflects this. Figure 2. Probability of a star both having and engulfing a massive planet on the RGB or AGB as a function of host star mass and metallicity, according to synthetic planet population models created for varying metallicity. The trend is for engulfment probability to rise with metallicity and for engulfment to be more likely to occur on the AGB for more massive host star. Finally, the effect of adopting Z * = 0.014 vs. Z * = 0.025 for the mean metallicity of planet-hosting stars, is that the planet-fraction increases from 3-9%, resulting in the fraction of non-spherical PNe doubled from 14-30%, respectively. Based on the mean metallicity of all planet hosting stars [18], the typical metallicity in orbital evolution models is Z * = 0.02 [6,19]. A statistical analysis of the final fraction of non-spherical planetary nebulae is very complex [20], but, based on many planet and stellar population synthesis models, we estimate that ∼15-30% of non-spherical planetary nebulae around single stars have a history in which they engulfed a massive planet on the AGB [9].

Discussion
Engulfment of massive exoplanets may contribute significantly to the formation of non-spherical planetary nebulae around single stars, yet appears insufficient to explain them all [9]. If additional physical mechanisms are not at work, it could be that the fraction of non-spherical PN is currently significantly overestimated.
It was demonstrated that the planet-fraction of all PN progenitors, depending on the adopted mean metallicity (between Z * = 0.014 and Z * = 0.025) of planet-hosting stars for the Milky Way, can range from 3-9% when adopting the distribution of the known exoplanets. The fractions translate to 15-30% of the population which would evolve to form non-spherical PNe. However, the current population of PN progenitors is not representative of the current population of PNe so estimating the planet-fraction from these results is not an adequate approach.