Numerical models of planetary nebulae with different episodes of mass ejection: the particular case of HuBi 1.

Planetary nebulae (PNe) are ionized gas envelopes that surround evolved low-intermediate mass stars. This gas was originally part of the stellar envelope. The nebular ejection processes occur during the red giant (RG) and Asymptotic Giant Branch (AGB) epochs and produce nebulae of different types of morphology, density distribution, chemistry, and kinematics depending on the characteristics and conditions of the progenitor star. The physical characteristics of such processes are still debatable, and they depend strongly on the stellar properties, being very different if the star is single, is in a detached binary system, is in a close binary system, and/or if it possesses an accretion disk.

Many PNe show evidence of several ejections at different epochs during the AGB and post-AGB phases; thus, PNe can present several shells of different morphologies, sizes, and other physical characteristics such as density, temperature, ionization degree, and even different chemistry.

On several occasions, it has been found that PN emission might be variable on short time scales or the PN might be ionized by a variable star (e.g., Balick, Guerrero, & Ramos-Larios, 2021; Otsuka et al., 2017; Hadjuk et al., 2015; Peña et al., 2022).

In the HR diagram, single central stars are expected to evolve from the AGB phase towards higher effective temperatures at nearly constant luminosity (see e.g., evolutionary models by Miller Bertolami, 2016 and others). The maximum temperature achieved could be a few hundred thousand K depending on the initial and/or final stellar mass. Afterward, the star enters the cooling track of white dwarfs, causing the evolved PN to slowly recombine while dispersing. However, in some cases, the central star experiences a late thermal pulse (LTP) or a very late thermal pulse (VLTP), in which case the evolution path reverses, the stellar atmosphere expands and cools, and the star makes a loop in the HR diagram, briefly returning to the AGB zone before continuing its evolution towards the PN phase a second time. Central stars that have experienced such an episode are called born-again (see e.g., Iben & Macdonald, 1995; Herwig et al., 1999; Miller Bertolami et al., 2006, Montoro-Molina et al., 2022 and references therein). It has been proposed that in returning to the PN stage, the star becomes a Wolf-Rayet, [WR], star and a H-poor nebula is ejected. However, not all the known [WR] central stars are born-again. It is tempting to propose that born-again stars are produced as a result of binary evolution. However, PNe are produced by low-intermediate mass stars and multiple lines of evidence point to binary systems as the birthplace of most PNe, then it should be expected that several of the eight so far known born-again stars belong to a binary system, although until now, only A 30 has been proposed as a possible binary star (Jacoby et al., 2020). As already said, so far only eight objects have been identified as born-again stars. In all these cases, it is crucial to analyze the evolution of the physical parameters of the nebula and its dynamical behavior, trying to understand the stellar evolution, which is not always directly observable due to the faintness of the central star.

HuBi 1 (PN G012.2+04.9, PM 1-188, IRAS 17514-1555) is one of the PNe whose central star apparently had a VLTP and presented a recent born-again episode (see Guerrero et al., 2018). Hu & Bibo (1990) discovered that the central star of HuBi 1 showed an unusually late [WC 10] stellar type, and later studies determined that it possesses an effective temperature of about 35,000 K (Leuenhagen & Hamann, 1998; Guerrero et al., 2018). The star has been fading with time, losing about 10 mag in the optical in the last 40–50 years. The nebula has two shells, one external, photoionized, emitting H and He recombination lines and weak heavy-element collisionally excited lines, and an inner shocked shell emitting intense collisionally excited lines of N⁺, O⁺ and S⁺. Pollacco & Hill (1994) were the first to mention that HuBi 1 consists of two shells with the inner, very compact one (diameter smaller than 1.5 arcsec), much denser and brighter than the outer one. This structure was confirmed by Peña (2005) and Guerrero et al. (2018) who showed the existence of this inner nebula in HuBi 1 based on images of H$\alpha$ and [N II]$\lambda$6584. They reported that the inner nebula has an inverse ionization structure typical of shocked gas. The inner nebula is responsible for the intense emission of [N II] lines while the outer nebula emits mainly H recombination lines. Rechy-García et al. (2020) and Peña et al. (2021) analyzed the nebular kinematics which shows the presence of an accelerated structure (with velocity larger than 100 km s^-1) in the inner region of the nebula with characteristics of shocked gas.

The nebula is schematically represented in Fig. 1 where the nebular radii, expansion velocities, and kinematic ages of the different components, as presented by Peña et al. (2021), have been included. Kinematic ages depend on the distance and radii adopted (different authors mention different sizes depending on the depth of their images). For HuBi 1 the heliocentric distance is 6.88$\pm$2.19 kpc, as derived by Frew, Parker, & Bojicic (2016). Ages therefore have large error bars, however the values we are using in Fig. 1 are in agreement with the values reported by Rechy-García et al. (2020) and Peña et al. (2021).

In this paper, we use hydrodynamic simulations to try to replicate the particularly intriguing characteristics of HuBi 1. Numerical simulations have become essential tools in astrophysics. Hydrodynamical models have been computed since long ago with the aim of understanding the physical processes of ejection and evolution of planetary nebulae by considering the stellar and environmental characteristics (Balick & Frank, 2002; Raga et al., 2000; among others). One 3D radiation-hydrodynamic simulation for HuBi 1 was generated by Toalá et al. (2021), by adopting mass-loss and stellar wind terminal velocity values obtained from observations. The model by these authors assumes that the inner shell of HuBi 1 was formed as a result of a very energetic VLTP which ejected H-poor material with velocities of about 300 km s^-1, 200 yr ago. They concluded that the large variations in the stellar wind parameters (mass ejection and the terminal velocity) produce instabilities and these instabilities can form structures that can be seen as clumps and filaments of ionized material between the shells of HuBi 1. This ionized gas produced the [N II]$\lambda$6584 emission. However, Toalá et al. (2021) have not a photoionization code or the atomic/ionic rates, therefore they could not directly conclude about the formation of emission lines produced by shocks, for example, the emission of [N II]$\lambda$6584. Moreover, the results they have obtained are only related to variations in velocity and mass injection rate during the VLTP phase. The temporal distance between the events is not entirely clear. This is because from the observations only the kinematic times are available and the formation of internal structures will be more dependent on the time between the events that allow for the interaction between the shells.

This paper is organized as follows: In §1 an introduction describing previous work on HuBi 1 is included. In §2 the numerical simulations for building the hydrodynamical models are introduced. In §3, we present the numerical models. Our results regarding the density structure and gas properties (velocity, pressure, and temperature structures) can be found in §4. In §5 we present our conclusions.

We developed a set of numerical simulations using a hydrodynamical spherically symmetric code named Coral1D, based on the code presented in Mellema et al. (1998) and Raga et al. (2000), where the photoionization and non-equilibrium radiative cooling functions are considered. The code solves the mass, momentum, and energy conservation equations.

where $\mathbf{U}$ is the vector containing the conservative variables which consider several ions of H, He, C, O, Ne, S and N which are calculated as advection equations. The vector of conserved variables is formed by

where $\rho$ is the density, $u$ is the velocity in the radial direction $R$, $E=\frac{1}{2}\rho u^{2}+\frac{nkT}{\gamma-1}$, is the internal energy, $n$ is the density, $T$ is the temperature giving by state equation $P=(n+n_{e})kT$ ($n_{e}$ is the electron density), and $\gamma$ is the ratio of heat capacities. The vectors are filled with the densities of the different ions, $n_{ion}$. $\mathbf{F}$ is the fluxes vector and is given by

and $\mathbf{S}$ is the sources vector. Therein, we calculate the heating and cooling rates and the ionization and recombination coefficients due to the photoionization and the emission of radiation. The atomic parameters used in Coral1D are described in Raga et al. (1997).

Our main aim in this work is to model PNe that have had more than one episode of matter ejection. In particular, we are interested in studying planetary nebulae that have presented ejection episodes with important differences in the momentum flux injected between episodes. As an example, in this paper we are going to focus on HuBi 1, a well known PN, recently analyzed by Guerrero et al. (2018) and Peña et al. (2021), among others, with at least three distinct episodes of matter ejection, as described below:

1. a)

   An external nebula, with a present radius of $\sim$10 arcsec as derived from the position-velocity diagram (PVD) presented in Fig. 6 by Peña et al. (2021) and Fig. 1 by Guerrero et al. (2018). It shows an expansion velocity of 50 km s^-1, and it was ejected about 6700 yr ago. The gas is photoionized, and in the process of recombination, due to the systematic weakening of the central star.

2. b)

   An inner nebula with a radius of $\sim$3 arcsec (same references as in item (a)), and an expansion velocity of 40 km s^-1, showing intense [N II]$\lambda$6584 emission that is evidence of being shocked. It was ejected about 1700 yr ago.

3. c)

   An outflow or jet immersed in the central part of the inner nebula with a size of about 2 arcsec and a velocity of more than 150 km s^-1; Rechy-García et al. (2020) claims a maximum velocity of 250 km s^-1 for this internal outflow. It seems to have been ejected about 335 yr ago.

All our models consider that there was a mass injection during the AGB phase of the progenitor star of the planetary nebula. The duration of the AGB phase is sufficiently long to fill the entire simulation box with gas that was injected at a velocity of 10 km s^-1, a mass injection rate of 1$\times$10^-5 M_⊙/yr, and a temperature of 200 K, typical of observed AGB winds (Kwok, Purton, & Fitzgerald, 1978). The ejected gas has solar abundance ¹¹1Kwitter & Henry (2022) have compiled the abundances of a large amount of PNe presented in several disk PN surveys and concluded that mean values of O, Ne and Ar abundances are close to solar, while He, C, N abundances are enhanced only by a few tenths of a dex, thus we have adopted solar abundances for the ejected gas in our models.. The simulation box is 0.325 pc per side, in a one-dimensional grid of 4000 solution cells, $\sim$ 16.72 au/pix. Using this simulation, we have constructed 2D projection maps of 500 cells per side, with a resolution of 131 au/pix.

The planetary nebula ejection has been divided into 3 stages:

1) First, a stage where the star has a luminosity of 10^3.8 L_⊙ and an effective temperature of 35,000 K (Guerrero et al., 2018). In this stage, a stellar wind is produced with a mass injection rate of 2$\times$10^-6 M_⊙/yr and a velocity of 360 km s^-1 (Kwok, Purton, & Fitzgerald, 1978). This stage is injected at a time t = 0 yr.

2) In the second stage, temperature and luminosity do not change with the first stage, as expected, but the wind that will form the inner nebula has a greater mass injection rate: 1$\times$10^-5 M_⊙/yr and a velocity of 50 km s^-1. This stage is injected at a time $t=\Delta t_{in}$ (see Table 1).

3) Finally, there is a last stage with a significant decrease in the injection of ionizing stellar photons, so we consider that the number of photons injected is 4 orders of magnitude lower than for the previous stages. This is consistent with the fact that the central star of HuBi 1 has declined by more than 10 mag in the visual, in the last 40$-$50 years. The wind injected at this stage has a velocity of 150 km s^-1 and a mass injection rate of 2$\times$10^-6 M_⊙/yr (like the first stage). This stage is injected at a time $t=\Delta t_{in}+\Delta t_{out}$ yr (see Table 1).

Table 1 shows the model properties, including the previous wind injected by the star in the phase of AGB, and the three different outflows that are considered in this work.

It is important to note that this study is only focused on quantifying the importance of each injection episode. The velocities and mass injection rates used in each of the episodes were selected to reproduce the observed kinematics in the three regions of the nebula HuBi 1, namely, the outer, inner, and jet nebula, as well as the currently observed velocity. To achieve this, we performed a series of simulations that allowed us to adjust these values within intervals proposed by a standard planetary nebula range from a hundred to a thousand km s^-1 (e.g., Kwok, Purton, & Fitzgerald, 1978; Perinotto et al., 2004; etc.). For the outer PN and for the inner PN nebula, we have explored velocity and mass loss rate values for the VLTP stage, such as those presented in Gomez-González et al. (2020) and Toalá et al. (2021) since tens to a hundred km s^-1. It is clear that there are other possible combinations of velocities and mass loss rates, for each episode, that can give similar results to those presented in this work. However, the physical parameters are within the expected intervals for each of the stages and the dynamic results we obtain with them show a good agreement with the observed values for this object. A more comprehensive study of velocities and mass loss rates for each stage can be conducted in future works, provided that the observational constraints we have for this object are considered, but this is outside the goals of this paper.

As we have said before, from observations it is found that PNe can have many injection stages, however, initially, we have explored the simplest models, models M0 and M1, in which not all the observed injection stages mentioned above are considered.

Our first simplest developed model, M0, considers a constant mass injection rate of 2$\times$10^-6 M_⊙/yr, an ejected gas velocity of 360 km s^-1, and constant effective temperature and luminosity during all lifetime of 6500 yr.

In model M1 we study the effect of the second stage of the material ejection, considering an increase in the mass injection rate and a decrease in the material injection velocity from the central source, as mentioned above.

Finally, to study the effect of the third ejection stage on the evolution of the planetary nebula, we have included a variation in velocity, mass injection rate and the number of photons injected by the central source (in the form of luminosity and effective temperature). For this stage, we consider that the star emits no ionizing photons and therefore, during the last 300 yr of evolution of this object, we have decreased the luminosity and effective temperature of the ionizing source by 10 orders of magnitude. We present the results of 3 models, M2, M3 and M4, for which we have only changed the duration of the second stage (900, 1000, and 1200 yr respectively). It is true that the observational articles propose a (kinematic) age of 1700 yr, which corresponds to the maximum age of the inner nebula. With models M2, M3 and M4, we explore the effect of this last ejection stage on the emission maps obtained from our simulations and compare them with the observations.

We have presented a set of models for the evolution of HuBi 1-like planetary nebulae, considering several stages of mass injection. We have carried out 1D numerical simulations using a hydrodynamical spherically symmetric code named Coral1D, where the photoionization and non-equilibrium cooling and heating functions are considered. In our models, the chemical abundances used for the plasma are solar. There was a previous work to model this object by Toalá et al. (2021). However, they have not included the non-equilibrium cooling and heating functions. Here we present a model that solves the hydrodynamics and the photoionization evolution of the gas, capable of studying not only the hydrogen but other 25 species (ions and neutral). In this paper, we present the [N II]$\lambda$6584 and H$\alpha$ lines to compare with the observations.

Our numerical models consider a previous AGB stage, where the medium has been structured by the supersonic gas blown during this stage. In this structure, there is a first stage of the planetary nebula in which wind-type material injected with a velocity of a few hundred km s^-1 and ionizing stellar photons emitted for approximately 5000 yr, forms the outer nebula that appears in HuBi 1. Then we added a new stage of a fast wind where the gas that is observed in the inner nebula will be injected for approximately 1000 yr. Finally, we injected the last stage where the central “jet-like” object is injected with an injection velocity much higher than in stage 2, but we have turned off the ionizing photon source (as seen in HuBi 1).

In this work, we have shown that the drastic fading of the central star determines the observed structure of this type of object. This is because the decrease in photon flux translates into a reduction of temperature and pressure in the internal structures, creating holes in pressure and forming many shock waves that interact with the gas being ejected by the central object.

As we mentioned above, we can not directly compare our PV synthetic diagrams obtained from the spherically symmetric models with the PV observed diagrams presented in Peña et al. (2021), but we can see that the gas from the latter ejection is shocked and is shown in the internal part, in the intense [N II]$\lambda$6584 line, as predicted in Peña et al. (2021), therefore this result is in a good agreement with the observations. The main reason for not being able to compare the PVs is that the last stage of the gas evolution is complex and since our model is spherical, we lose the geometry. The observed object has a high speed ejecta in a cylindrical or jet shape. We can infer this from Figure 6 in Peña et al. (2021), where they show two slits at different positions (0^∘ and 90^∘) where we can compare their right upper panel (position at 90^∘) with the high speed ejecta ($\sim$150 km s^-1) and their right lower panel (position at 0^∘) with no high speed ejecta, this could give us a hint of the central object morphology as a “belted outflow”. Moreover, Rechy-García et al. (2020) performed two “pseudo-slits” (as they named them) and constructed the PV diagrams at two different angles 50^∘ and 140^∘, showing high velocity outflows ($\sim$200 km s^-1). If we join all these observed angles, we can picture a ”belted outflow” shape.

Finally, we present the model for HuBi 1 as an example to show the capability of our model to explain the hydrodynamical and photoionization evolution in ionization nebulae. Although Coral1D is a code presented by Raga et al. (2008) with an excellent result, this is our first step with a 1D code to study these two physical phenomena at the same time. We are working on a 2D code to analyze different morphologies and different kinds of central ionization sources.