EDP Sciences
Free Access
Issue
A&A
Volume 585, January 2016
Article Number A69
Number of page(s) 14
Section Interstellar and circumstellar matter
DOI https://doi.org/10.1051/0004-6361/201526653
Published online 21 December 2015

© ESO, 2015

1. Introduction

A planetary nebula (PN) forms when a low- or intermediate-mass star at the end of its life ejects much of its envelope into space. The nebula initially expands with the ejection velocity. As the central star heats up, it ionizes the nebula, and the resulting ionization front and the over-pressure of the ionized region alter the velocity field and density structure. Interaction with the fast wind from the hot central star causes further changes. The velocity field therefore contains information both on the origin and the evolution of PNe.

For spherically symmetric nebula, these interactions have been modelled by the Potsdam group, e.g. Jacob et al. (2013), Schönberner et al. (2014) and references therein. These articles discuss a relation between the true nebular expansion and expansion velocities estimated in different ways. The different kinematical ages are compared with the true post-AGB age. Schönberner et al. (2014) and Gesicki et al. (2014) independently obtained the same result: for more or less regular PNe, the true post-AGB age and kinematical age obtained from the mass-averaged expansion velocity agree with a correction factor of 1.4. In both papers, the “true” post-AGB age is the elapsed time since the star left the asymptotic giant branch (AGB), as obtained from the evolutionary computations of Blöcker (1995).

Observationally, the velocity fields can be derived from the line profiles of different ions where each ion traces a different layer due to the ionization stratification. Kinematical reconstruction of line profiles is done by fitting a photoionization model to the observed nebular parameters and then finding the velocity field that provides the best agreement between the calculated and the observed emission line profiles. The suitably averaged velocity and the outer radius result in kinematical age. Such methods were used successfully for spherical models of PNe in Gesicki et al. (2006), using 1D photoionization models. Using HST imaging combined with VLT spectroscopy, Gesicki et al. (2014) studied the velocity fields of a sample of 36 compact PNe in the Galactic bulge direction. This study also used the 1D Torun photoionization codes, and thus implicitly assumed spherical symmetry. However, most planetary nebulae deviate from spherical symmetry.

The evolution of bipolar PNe is quite successfully described in terms of what is known as the generalized interacting stellar winds model. Balick (1987) schematically outlined the proposed morphological sequences of such an evolution. Recently Huarte-Espinosa et al. (2012), using extensive hydrodynamical simulations, studied the morphological evolution of bipolar nebulae. More complex interactions, such as the wind reflected by a warped disk, allows more complex nebular shapes to be generated (Rijkhorst et al. 2005). However, full 3D radiation-hydrodynamic calculations are very time consuming and completed evolutionary sequences, like those of Potsdam group for 1D, are not available yet.

In this paper for the first time we study the multidirectional structure and velocity field for bipolar PNe. We reconstruct the images and the emission line profiles following a procedure similar to the one elaborated for spherical objects, which is suitable for modelling a larger sample. After deriving Hβ-weighted averaged velocities, we estimate the kinematical ages.

For the kinematical reconstruction we use the concept of pseudo-3D models described in Morisset et al. (2005). The basic idea is to compute a grid of 1D photoionization models with each model corresponding to the given direction (θ, φ) pointing away from the centre of the nebula. A 3D cube of coordinates is built, on which interpolation of any physical parameter is made based on the 1D runs: for each cell of coordinates (r, θ, φ), the program looks for the three closest 1D models that enclose θ and φ, and interpolates on the radii to find the value corresponding to the given cell. In the case of an axially symmetric nebula, there is no need for interpolation on the longitudinal angle φ and only the latitudinal angle θ is used. These models are as accurate as full 3D models as long as the non-radial radiation does not dominate the ionization. The contribution from non-radial directions would be significant in shadowed regions ionized by radiation from the surrounding material or in a nebula ionized by multiple sources, which are not the cases here.

After a short presentation of the observational material we describe the applied pyCloudy python library and the pseudo-3D models derived with it. For each of the selected objects the best-fitting model will be presented in separate figures and subsections. Next we search for correlations between the new parameters and discuss them. Finally, we propose an evolutionary scenario in which the lobes evolve surprisingly fast – the obtained timescales for the bipolar phase are an order of magnitude smaller than the often assumed visibility times (Jacob et al. 2013).

thumbnail Fig. 1

Hα HST images presented at the same angular scale. Each panel is 12 arcsec per side. The intensity scale is linear. The names are indicated in each panel. The first three columns show PNe in the bipolar phase while the rightmost column shows the proposed post-bipolar phase (see text). Six PNe are assumed to belong to the Galactic bulge while the smallest and largest, both shown in the third column, are assumed to be located behind and in front of the bulge, respectively.

Open with DEXTER

2. Observations and analysis

2.1. Selected images and spectra

The observational data are presented in Gesicki et al. (2014) and described there in more detail. The images of compact bulge PNe were obtained with HST/WFPC2 in Hα and [O iii] 5007 Å. For each object a series of 600-second VLT/UVES echelle spectra were obtained with the slit centred on the nebula and aligned in most cases with the minor axis.

The sample of 36 bulge PNe with sizes of 5 arcsec or less was randomly selected from compact PNe listed in the SE CGPN (Acker et al. 1992), and is not biased towards any morphology or brightness. Among these PNe there are 21 objects composed of a central structure (ellipsoid or disk) plus a clear extension. These extensions can be “classic” bipolar lobes (6 PNe), multipolar (quadrupolar or more) lobes (10 PNe), or curved wing-like open structures (5 PNe). Some other objects are ellipsoidal with or without inner structures or are completely irregular. A detailed morphological classification can be found in Sahai et al. (2011).

From this observed sample, for the 3D analysis we selected objects that appear similar and exhibit both axial-symmetry and point-symmetry. These PNe show both a central structure (which can be interpreted as an equatorial torus) and a pair of rather wide, extended and closed classic lobes. Both lobes are the same and are oriented perpendicular to the equatorial plane. This criterion is fulfilled by six objects only: Hen 2-262, H 2-20, H 2-18, Hen 2-260, Wray 16-286, and H 2-27. It is not unlikely that some other objects in the complete sample should fall into this group. Round PNe could be bipolars seen along the symmetry axis. However, for the clarity of this analysis we did not consider these doubtful cases. We added two more objects to the present work: H 2-15 and H 2-26. Because of their morphology and the old age assigned by Gesicki et al. (2014), we suggest that both are in a post-bipolar phase and thus provide a useful comparison sample; this is discussed in Sects. 3.7 and 4.4.2.

Figure 1 presents the Hα HST images. The intensity is represented as a linear grey scale. The images reveal some details of the bipolar lobes; the lobes are much weaker in intensity than the torus. The angular size of each square is the same (12 arcsec per side) to visualize the differences in angular size between the PNe. The two proposed post-bipolar PNe are shown in Fig. 1 in the rightmost column.

The selected emission lines are shown in the right columns of Figs. 310 together with the fitted models. The spatial resolution of the UVES data is much lower than that of HST: the slit was 0.5 arcsec wide and the seeing was about 1 arcsec (typically between 0.6 and 1.4 arcsec, for H 2-27 even above 2.5 arcsec). For the analysis of the line profiles, the central arcsecond was extracted from the observed slit spectra. Thus, the analysed line profiles are those observed near the centre of each nebula. This is enough to derive the velocity of the inner regions, but the far-away tips of the lobes lie beyond the long-slit so that the extremum velocity is not as well constrained by the observations.

2.2. pyCloudy and 3D photoionization modelling

The pyCloudy code (Morisset 2013) used for the modelling is available on a web page1 together with documentation and instructions. It consists of a photoionization code (for which Cloudy2 (Ferland et al. 2013) is used) and a python library that contains routines for analysing the models. Included are routines for visualizing the models, for computing the differently defined averaged expansion velocities, for integrating the masses, fluxes, etc.

As already described, this method of solving for the structure is pseudo-3D rather than full 3D. The python codes in their current version limit the modelling to axial- and point-symmetric objects. The photoionization models are themselves 1D, and therefore do not consider the effect of the diffuse radiation field correctly. The interpolation introduces further uncertainties. It is, however, computationally fast. Full models, e.g. using Monte Carlo methods (see e.g. Danehkar et al. 2013), are computationally very expensive. In the current work we explore whether the data can be fitted with the pseudo-3D models. Full 3D models are left for later.

We use a black-body spectrum for the central star, as was done in previous 1D models. These models tend to return a lower stellar temperature than using model atmospheres, because they lack the Lyman jump. We ran several tests comparing black-body SEDs with the stellar atmosphere SEDs available at the Cloudy web page, in particular with the PN central star models of Rauch (2003) and found that this had no significant effect on the velocity fields and ages. The adopted black-body temperatures Tbb are listed in Table A.1.

2.3. Density structure

The photoionization part of pyCloudy computes a sequence of 1D regular Cloudy models aimed at representing the structure of a single nebula diagnosed at many different radial directions. Each 1D model actually represents a spherically symmetric PN and consists of an empty inner region extending from the centre up to the specified inner radius, then follows the main nebular body at a constant density, and terminates at the specified outer radius or at the ionization boundary, whichever comes first. Produced are 1D grids of line emissivities, electron temperatures, etc.

Each computed 1D model corresponds to a selected radial direction, inclined to the equatorial plane at a specified (latitudinal) angle with varying inner and outer radii. In our models, the first is at (equator) and the last at 90° (symmetry axis). The angles in between are selected based on the structure of each nebula, and are different for each model. The spacings do not need to be uniform in angle. We used eight separate angles, based in part on the required computational time. A test run with a set of 16 angles showed little improvement, and for our structures sets of 8 angles appear to be sufficient.

pyCloudy re-grids the obtained 1D distributions into a XZ cross-section of a 3D cube. Because of the introduced symmetry centre and symmetry axis, only one quadrant of the plane needs to be filled. The space between the 1D lines is filled in by interpolation while the remaining quadrants of the cross-secting plane are filled in by symmetry conditions. These planes are shown for each analysed object in Figs. 310 in the upper left panels. The symmetry axis is placed horizontally. The gas density is shown in shades of grey in linear scale but is plotted only for the ionized part of the nebula. In one quadrant the magenta line segments indicate the eight directions along which the 1D Cloudy models were computed. The outer end of each segment corresponds to the assumed end of gas distribution i.e. the matter (density) boundary. Wherever the shown density does not reach the end of the line segment, it means that the actual nebular edge is defined by the ionization front in this direction.

The full cube is filled by pyCloudy through rotation around the symmetry axis. To fill the cross-secting plane and the full cube the standard linear interpolation of python was applied. Routines are available to compute the projection of the full cube onto the plane of the sky, defined by the inclination angle between the line of sight and the symmetry axis. In this way model images can be obtained for any of the calculated emission lines. For Hα they are shown in Figs. 310 in the centre left panels. Very different nebular structures can be obtained by varying the shapes of the inner voids and outer lobes.

The computed images are compared visually with the HST images (shown in Figs. 310 in the bottom left panels) and then the inner and outer radii, the inclination of the symmetry axis, the densities, and the latitudinal angles of the individual sectors are optimized by trial and error. Their numerical values are given in Table A.2.

2.4. Velocity field

A velocity field is defined where all velocities are directed radially with respect to the central star, and the velocity only depends on the radial distance from the centre. We note that this is simpler than the density structure: the material has different spatial distribution at different latitudinal angles while the velocity distribution is the same for any inclination. More complex velocity fields (depending for example on the angular distance from the symmetry axis) are possible, but go beyond the scope of this paper or the observational constraints.

A constant expansion velocity is inconsistent with the different widths of lines of different ions, which indicate higher velocities in the outer (low ionization) regions. To limit the number of the model free parameters, the velocity was restricted to be linearly increasing with radial distance, starting from zero at the central star and reaching a maximum at the ends of lobes. This Hubble-type dependency (alternatively called homologous) improved the model. Trial computations showed a need to further modify this simple scheme. We introduced an inner region with a constant velocity. Thus, the adopted velocity fields are composed of an inner plateau (its velocity and radial extent was adopted by trial and error, guided by the [O iii] 5007 Å line profile), connected to the linearly increasing part further out (guided by the [N ii] 6583 Å line). The spherical symmetry of the velocity field combined with axial and point symmetry of the density distribution provides a slow expansion for the dense regions close to the star (dominated by the equatorial torus at low inclination angles) and fast expansion for the far-away extended and thin lobes at high inclination angles.

We stopped at this level of complication. Such a situation appears reasonable for 3D bipolar models and corresponds quite well to many of the 1D models (composed of a dense inner shell and an extended outer halo; e.g. Gesicki et al. 2014). Such velocity is in agreement with the “typical middle-aged and fully ionized nebular structure” of Schönberner et al. (2014) obtained with the 1D hydrodynamical modelling. Even better agreement is with the 2.5D hydrodynamical calculations of Huarte-Espinosa et al. (2012) where some of the models show a Hubble-type expansion combined with closed bipolar lobes.

Our velocity field is also similar to those found observationally for bi-lobed nebulae. For the well-resolved nebula NGC 2346 at a distance of 700 pc, Arias et al. (2001) derived a best model with “outflow velocity directed radially and proportional to the radial distance”. The NGC 2346 expansion velocities of 16 km s-1 at the central ring and of about 60 km s-1 at the tips of the lobes are comparable to many of our values. Similar velocity fields were found for two other nebulae, NGC 6302 (Szyszka et al. 2011) and NGC 7027 (Walsh et al. 1997), although both objects have a rather dramatic appearance. Alcolea et al. (2007) analysed molecular gas in the pre-planetary nebula M 1-92 (Minkowski’s footprint) and found a Hubble-like velocity field not only in the bipolar flow, but also in the equatorial plane. We note that this situation is neither new nor unique and is typical to our whole sample.

Line profiles can be calculated by integrating the particular line of sight through the cube, at any position within the projected image. Summing these over a defined area simulates the line profile obtained through a slit. The spectrograph slit is assumed to be in the shape of a rectangle around 1.5 arcsec per side. This includes the effect of seeing. Positioned centrally on the nebular image, this produces symmetric line profiles. However, when the rectangle is shifted somewhat off centre, calculations show an asymmetric profile. This effect is well known, see e.g. Monteiro et al. (2000). The asymmetry obtained in this way fits quite well to the observed spectra and this was done for some of our objects. In Figs. 310 in the middle left panels we indicate with the (red) rectangle the applied size of the slit superposed on the projected nebula image. The right panels show the computed line profiles (dashed red lines) obtained with this set-up.

However, there is an ambiguity: we cannot tell whether the observed asymmetry comes from a slight offset of the spectrograph slit (which for extended weak objects is very likely) or from real irregularities within the nebulae themselves (which the images indicate is also present in some cases).

2.5. Inclination angle

The derivation of geometrical dimensions and of the expansion velocities for bipolar PNe is much more complicated than for spherical ones. The angle of inclination between the polar axis and the plane of the sky affects several parameters and requires special attention.

thumbnail Fig. 2

Emission line profiles for different inclination angle. Plotted is the [N ii] 6583 Å line computed for the model of the spatially unresolved nebula Hen 2-260. The solid line represents the inclination angle of (symmetry axis within the plane of the sky), long-dashed line: inclination of 30°; short-dashed line: 60°.

Open with DEXTER

The dominant part of the line emission originates in the dense nebular torus. If the symmetry axis lies in the plane of the sky then most of nebular material moves towards and away from the observer and the resulting emission profile is clearly split with a dip at the systemic velocity. At the other extreme, when the torus lies in the plane of the sky the bulk expansion proceeds perpendicularly to the line of sight, the line narrows, and the splitting disappears. The observed splitting of the brightest part of the line profiles can be used to derive the best-fitting inclination of the symmetry axis of the model as illustrated in Fig. 2. This figure was plotted for [N ii] 6583 Å line of Hen 2-260 for inclination angles of , 30°, and 60°. This plot can be compared with Fig. 6 (bottom right panel) where an inclination of 50° was adopted. The nebula is spatially unresolved by the spectrograph; therefore, the behaviour shown in Fig. 2 does not depend on the selection resulting from an applied small slit area. This kind of analysis is possible thanks to the high spectral resolution (about 5 km s-1) and this example shows that the accuracy of the inclination of such a procedure can be roughly estimated as ± 10°. However, the splitting also depends on the assumed velocity field. Furthermore, the derived inclination angle determines the true extent of the model which is used in the velocity field.

Another effect can be deduced from Fig. 2: with increased inclination angle the low-intensity and high-velocity line wings become apparent. This follows from the fact that for the highly inclined object the extended fast lobes expand in the direction closer to the line of sight and the emission becomes Doppler-shifted. This effect causes the lobes to be partially available for the kinematical reconstruction, even despite the limited spectrograph slit area. For the angles of about 40°−50° a significant part of the velocity field can be probed. However, for smaller inclination angles and larger nebulae the lobe’s velocity has to be extrapolated with the assumed homologous expansion law derived based on the central regions. Certainly integral field unit spectroscopy would provide data much better suited to this kind of analysis; nevertheless, even this would not help to measure the gas motions in the plane of the sky.

All these factors introduce a type of degeneracy. The degeneracy is reduced by the restriction to linearly increasing velocity fields and ellipsoidal shapes. Careful examination of details of the images and spectra assures us that the model is a good proxy for the real nebulae, but each solution may not be unique.

3. pyCloudy models for individual objects

thumbnail Fig. 3

Hen 2-262. The top left panel shows the assumed density structure only for the ionized nebula, drawn at the cross-section plane taken along the positioned horizontal symmetry axis. The magenta lines indicate the directions along which the 1D constant density Cloudy models were calculated. The end of each line segment corresponds to the assumed matter (density) boundary. The middle left panel shows the Hα image (i.e. a projection onto a plane of the 3D structure) obtained for the same model, but inclined to the sky plane at 20° to correspond approximately to the HST image shown in the bottom left panel (also in Fig. 1). The intensity grey scale is linear in all three images, the size is scaled to fill the width of the box. The red rectangle in the middle panel indicates the area over which the line profiles are integrated. The nebular emission line profiles are shown on the right; the solid black lines show the observations and the dashed red lines show the computed model for Hα, [O iii] 5007 Å, and [N ii] 6583 Å.

Open with DEXTER

thumbnail Fig. 4

H 2-20. Data presented as in Fig. 3. The axis inclination is 30°.

Open with DEXTER

thumbnail Fig. 5

H 2-18. Data presented as in Fig. 3. The axis inclination is 35°.

Open with DEXTER

thumbnail Fig. 6

Hen 2-260. Data presented as in Fig. 3. The axis inclination is 50°. The spectrograph slit (left middle panel) covers the whole nebula.

Open with DEXTER

thumbnail Fig. 7

Wray 16-286. Data presented as in Fig. 3. The axis inclination is 20°.

Open with DEXTER

thumbnail Fig. 8

H 2-27. Data presented as in Fig. 3. The axis inclination is 35°.

Open with DEXTER

thumbnail Fig. 9

H 2-15. Data presented as in Fig. 3. The axis inclination is 45°.

Open with DEXTER

thumbnail Fig. 10

H 2-26. Data presented as in Fig. 3. The axis inclination is 40°.

Open with DEXTER

Table A.1 lists the parameters of the photoionization models for the program nebulae. The Galactic bulge distance was assumed by default; however, the two objects discussed below were considered to be lying outside the bulge. The stellar black-body temperature Tbb was adopted from previous Torun models (Gesicki et al. 2014). We checked that major modifications were not necessary to improve the fits, although we did allow for minor corrections. The stellar luminosities were adopted to fit the nebular Hβ fluxes. In some cases these luminosities were smaller than predicted by the evolutionary tracks that Gesicki et al. (2003) interpreted in terms of photon leakage through transparent parts of the nebulae. In the axially symmetric objects the low-density and low-opacity lobes are easily identified as being responsible for such leakage therefore the luminosity values should be regarded as lower limits to the true stellar values.

For the line strengths and ratios we list the observed and the pyCloudy values. The observed values are taken from the literature and relate to the entire nebula. The Hβ fluxes are fitted well because the densities and luminosities were varied to obtain that fit. The two most representative strong lines are fitted with an average error of about 25% (in two extreme cases reaching 60%), which is actually not larger than in the old Torun models. The chemical abundances, are taken from Table 7.2 (column “Planetary”) of the Cloudy C13.1 documentation available at the web page. For N and O we changed the values to the available estimations for individual PNe wherever possible; these data are given in Table A.1.

The ionized mass is defined as the mass of the region where hydrogen is ionized, and includes both the ionized and neutral helium in this region3. For each object, two values of the ionized masses are listed: the 1D Torun model and the pyCloudy model. The difference between them is discussed in Sect. 4.1.

The values of the inclination angles between the sky plane and the nebula symmetry axis are presented in Table A.2 and in the figure captions for each of the discussed objects. In Figs. 310 we present the influence of the inclination angle on the observed image: the upper left boxes show the density distribution over the full length of the model, while the middle left boxes show the Hα image obtained for the same model but now with the long axis inclined to the plane of the sky. The effect is more or less visible depending on the assumed angle. These images correspond to the appropriately rotated HST images shown in the lower left boxes. This double presentation illustrates more clearly the intrinsic model structure versus the observed structure of the nebula.

3.1. Hen 2-262 (PN G001.2+02.1)

The Torun model of Hen 2-262 has not been previously published and therefore its parameters are not compared with pyCloudy in Table A.1. It was observed with the same set-up together with the whole Galactic bulge sample; however, it was not included in the analysis of Gesicki et al. (2014) because the complicated velocity field required significant turbulent motions in the 1D model. The pseudo-3D model shown in Fig. 3 is acceptable; moving the slit off-centre improved the fit, but it still appears that adding some turbulence would broaden the line wings and reduce the central splitting in better agreement with observations. However, the presented model is with a simple, linear velocity field.

The density contrast between equator and pole is a factor of 5.6. The model is radiation bounded in the equatorial torus, while the lobes are matter bounded.

3.2. H 2-20 (PN G002.8+01.7)

The best-fitting pyCloudy model (see Fig. 4) has an equator-to-pole density contrast of a factor of 5. The model is radiation bounded in all directions.

The emission line profiles show the presence of an ionization gradient within the PN, which is usually associated with the presence of the radiation boundary, i.e. ionization front. A small offset of the artificial slit assumed at the line profile computations resulted in asymmetric lines which resemble the ones observed. Nevertheless, in Fig. 4 more asymmetry can be seen in the observed profiles which cannot be fitted with cylindrical-symmetry models.

The 1D model of Gesicki et al. (2014) predicted a velocity field that monotonically, almost linearly, increases outwards with a small plateau in the inner region. This is very similar to the velocity field derived in this work.

3.3. H 2-18 (PN G006.3+04.4)

The best-fitting model (see Fig. 5) has a high equator-to-pole density contrast of a factor of 11. The [N ii] 6583Å line is very weak: any model with an ionization front overestimates this line by an order of magnitude. The observed line ratios are best reproduced with a model that is matter bounded and has low densities. However at these low densities there is no ionization gradient within the nebula, which contradicts the observed emission profiles. We built a compromise model which is still matter bounded in all directions but with a density of the equatorial torus that is high enough to show an ionization gradient within the nebula.

The shape of the three emission lines is reproduced reasonably well, including the unusual quadruple [O iii] 5007 Å line. However, the details of [N ii] 6583 Å cannot be better tuned within the simple model. The 1D model of Gesicki et al. (2014) predicted a linearly increasing velocity field.

3.4. Hen 2-260 (PN G008.2+06.8)

The equator-to-pole density contrast is a factor of 7. The best-fitting model (see Fig. 6) shows a larger than average inclination, which is required to reproduce the lack of splitting in the [N ii] 6583 Å line profile. The nebula is radiation bounded, except for a narrow cone along the symmetry axis where it is matter bounded.

This object was discussed in Hajduk et al. (2014), where it was analysed with spherically symmetric 1D Cloudy and Torun codes. They derived a distance to this PN of 12 kpc. We adopted this value, and also took the other parameters from Hajduk et al. (2014). The object was also in the sample of Gesicki et al. (2014), but there it was assumed to be at the distance of the Galactic bulge and therefore a smaller radius was used. The low radial velocity of Hen 2-260 (Durand et al. 1998) is consistent with a non-bulge location. Figure 1 shows that Hen 2-260 is much smaller than the other PNe, which is understandable if it is located beyond the Galactic bulge (Bensby & Lundström 2001).

The line ratios are reproduced rather well by pyCloudy using a black-body Teff = 26 000 K. Hajduk et al. (2014) used a stellar atmosphere model which requires a higher temperature for the same number of ionizing photons: they find 35 000 K. Hajduk et al. (2014) derived the age of this PN with a direct method (from the heating rate of the central star): their age (980 yr) approximately agrees with that estimated below from the kinematics.

The 1D model of Gesicki et al. (2014) predicted a parabolic shape for the velocity field, with an average expansion velocity somewhat smaller than the 3D estimation presented here. The linear velocity field used here provides a reasonable fit to the line profiles, but there are slight discrepancies in the widths of the various lines. Any further improvement would require a more complex velocity field.

3.5. Wray 16-286 (PN G351.9-01.9)

The equator-to-pole density contrast of the best model (see Fig. 7) is high, at a factor of 9.3. The model is radiation bounded in the equatorial torus, while the lobes are matter bounded.

The old 1D Torun models indicated a parabola-like velocity field with mass-averaged velocity of 13 km s-1. This yields an age that is significantly larger than pyCloudy: we found that the parabola minimum is responsible for the difference. The pyCloudy model reproduces line profiles much better and the narrow [O iii] 5007 Å and wide, flat-topped, and slightly asymmetric [N ii] 6583 Å are explained by the spectrograph pointing off-centre with the slit not covering the midpoint. The pyCloudy model results in a higher average velocity of around 20 km s-1 (depending on the weighting method used in the averaging).

3.6. H 2-27 (PN G356.5-03.6)

The best-fitting model (see Fig. 8) gives a high equator-to-pole density contrast of a factor of 7.9. The model assumes a spectrograph slit position, which is a little offset in order to reproduce the small asymmetry observed in profiles. The model is matter bounded in all directions.

For this PN we assumed a smaller distance than to the bulge. Although reasonable photoionization models can be built for different distances, some aspects point towards a smaller value. In Fig. 1 this PN is shown to have the largest angular dimensions. Durand et al. (1998) report a low radial velocity, which does not require a bulge location. The correlations discussed below suggest that H 2-27 has a smaller kinematical age. The velocity is derived from spectra, i.e. independently from distance; therefore, the age can be reduced by decreasing the physical radius and therefore the distance. Having no more constraints we simply assumed the distance of 4 kpc for the pyCloudy models. For H 2-27 we recomputed the Torun model for the smaller new distance; previous Torun models, which assumed bulge distance, resulted in a greater age for this PN.

The 1D model of Gesicki et al. (2014) predicted a linearly increasing velocity field.

3.7. Two post-bipolar nebulae: H 2-15 (PN G003.8+05.3) and H 2-26 (PN G356.1-03.3)

These two nebulae have a more irregular structure than the others, and it has been suggested that they are in a post-bipolar phase of evolution. They have numerous common characteristics and they are included here for comparison. The models are shown in Figs. 9 and 10.

The models show a relatively low equator-to-density contrast of a factor of 5 and 2.5 for H 2-15 and H 2-26 respectively. They have lower density than the other objects peaking at around 103 cm-3. Both models are composed of a thick disk which is radiation bounded in the directions near the equatorial plane, while the directions towards the poles are matter bounded. For both stars, the lowest stellar luminosities were derived so they can be the most transparent for stellar ionizing radiation (most leakage). Both were fitted with linear velocity fields in Gesicki et al. (2014). The 3D pyCloudy resulted in a very similar averaged velocities. Gesicki et al. (2014) built models for these two PNe with a low density halo attached to the dense shell; this improved the fit of the [N ii] 6583 Å line, but in consequence increased the outer radius and the age. The pyCloudy models assume for both nebulae that the artificial spectrograph slit is a little offset to reproduce the asymmetry observed in profiles. However, the images also show significant asymmetries. The kinematical ages derived here are smaller than those derived in Gesicki et al. (2014), but are still two times larger than the ages of the six bi-lobed PNe.

4. Discussion

4.1. Ionized masses

For all objects the pyCloudy masses (presented in Table A.1) are smaller than those from the 1D Torun models, in most cases (6 out of 8) by a factor of 2–7. In two cases the masses are nearly the same. The dominant effect is the use of a 3D versus a 1D structure. A 3D model provides for a range of density distributions within one object.

The Hβ emission that is proportional to the density squared (actually np × ne) is dominated by regions of high hydrogen density. This gas has different spatial distributions in the spherical and the pseudo-3D models, being a compact sphere or a torus, respectively. Despite these differences the densities and overall masses of these regions are comparable and produce the same Hβ fluxes. The differences in the distributions of low density gas are much larger. In spherical models the low density gas occupies a large volume of the spherical halo, which is a significant fraction of the total nebular mass.

Another contributing factor is the irradiation of the outer regions. In 1D models, the radiation reaching the outer regions first passes through dense inner layers and the attenuation gives a much reduced intensity. To reach a significant line strength of the low-ionization lines (which arise in the outer regions), a larger amount of gas is required. In 3D models, there are directions in which the outer regions are illuminated by radiation which has passed through little intervening gas. Thus, less gas is needed to reproduce the line strengths.

This important result shows that ionized masses of planetary nebulae have very large uncertainties attached. One-dimensional models are best viewed as giving upper limits to the mass. Pseudo-3D models provide a better approximation to true nebular ionized masses. The total mass lost by the central star might be somewhat larger than the model’s ionized mass because we do not know how much neutral material exists in the dense equatorial disk beyond the ionization boundary.

4.2. Averaged expansion velocities

The true velocity distribution within a planetary nebula can in principle be quite complicated. The current analysis considers only the large-scale flows which in our paper are only approximated, as described in Sect. 2.4. Further simplification is related to limiting the velocity as a function of radius only, and not of latitudinal angle. Defining a single-valued expansion velocity for a non-constant velocity field requires a weighting. The velocity at the outer radius (or ionization front) is typically not useful as it is not directly measured; the available emission lines come from a region a bit further in (e.g. the “post-shock” velocity discussed by Jacob et al. 2013). For 1D models, the mass-averaged expansion velocity (Gesicki & Zijlstra 2000) has been shown to be a robust parameter. It represents the bulk motion, and is easily and reproducibly measured from the models, even where very few emission lines are available.

For the 3D models, a similar methodology has not been developed yet. We investigate two options: weighting by mass (similar to the 1D models) and weighting by the Hβ emissivity. Weighting by mass has the drawback that the pyCloudy models use constant density for each direction. This gives a high weight to the outermost radius where velocities are highest but neither the densities nor the velocities are well constrained. The Hβ emissivity decreases outwards tracing the ionized mass distribution. Weighting by Hβ emissivity gives smaller weight to the outer radius – a behaviour closer to the 1D Torun models – where the density decreases towards the outer radius.

Table A.1 lists the averaged expansion velocities Vexp and nebular radii together with the corresponding 1D Torun model values. Both averaging methods are compared in Fig. 11, plotted versus old Torun model values. The pyCloudy mass-averaged velocities are higher than the Hβ-averaged velocities – exactly as expected. However, we note that in either case the result is dominated by the torus where most of the gas is located.

thumbnail Fig. 11

Averaged expansion velocity from Torun models compared with pyCloudy. The abscissa shows the Torun models mass-averaged velocities, while at ordinate the pyCloudy Hβ-averaged values are shown as filled symbols, and the pyCloudy mass-averaged as open symbols. The (red) circles show the bipolar nebulae and the (blue) squares the post-bipolars. Data from Table A.1.

Open with DEXTER

Comparing the values of Vexp from pyCloudy averaged by Hβ with the values of Vexp from Torun models averaged by mass indicates a good correspondence (see Fig. 11). We conclude that the averaged velocity remains a robust parameter for describing a planetary nebula. Because of the rather scarce observational and hydrodynamical data for comparison, weighting by the Hβ line emissivity is preferred over weighting by mass for the simplified models used here.

thumbnail Fig. 12

Averaged expansion velocity vs. averaged outer radius for the pyCloudy modelled PNe. Velocity are Hβ-averaged and the outer radius is mass averaged (see text). Data from Table A.1. The filled circles show the bipolar nebulae and the filled squares the post-bipolars.

Open with DEXTER

In Fig. 12 we plot the velocities versus the outer radius radii. The velocities shown are Hβ averaged while the outer radii are mass-weighted. The circles represent the bipolar nebulae. They show almost linear relation. The two PNe tentatively classified as post-bipolar are represented by squares: they have the largest radii but also show low expansion velocities.

4.3. Kinematical ages

The kinematical age can be derived from the (averaged) expansion velocity and the outer radius. However, this only provides an approximation. The ionization front travels faster than the gas velocity: the gas density declines as the nebula expands, allowing more gas to become ionized. The velocities can increase very quickly at the ionization front owing to the over-pressure. Quantifying the effects requires hydrodynamical models. Based on the hydrodynamical models of Perinotto et al. (2004), in the 1D case the kinematical age should be reduced by a factor of 1.4 (Gesicki et al. 2014).

For 3D models the derivation is more complicated. First, there is no unique outer radius. Second, the torus and the lobes may also have different ages. Third, the correction factors from 1D hydrodynamical models may not be directly applicable to a non-spherical 3D structures where pressure gradients may cause non-radial velocities. Because of these uncertainties, we focus on deriving limits to the ages.

An upper limit can be obtained from the outer radius and the expansion velocity of the torus. We calculate the mass-averaged outer radius of the ionized nebula. This is dominated by the dense torus; the value of rout is listed in Table A.1. We divide this radius by the Hβ weighted expansion velocity to obtain an age. The Hβ-averaged value is used because it is smaller than the mass-averaged value, and therefore yields an older age. This age is listed as “nebula” age, and should be considered an upper limit. No correction factor other than the averaging has been applied. For comparison, the 1D Torun model ages are also listed; for these the correction factor of 1.4 has been used (Gesicki et al. 2014; Schönberner et al. 2005).

A lower limit can be obtained from the lobes, which are expected to be younger than the main nebula (Huggins 2007). The linear velocity field indicates a Hubble-type flow in the lobes. The inherent feature of homologous expansion is the same age derived at any distance, therefore the exact estimation of the extent of the lobes is not an issue. We assume that no acceleration or deceleration has taken place and obtain the age for the tip of the lobes. These ages are listed as “lobe” age in Table A.1.

Our “nebular” age is different from the Hubble-type age of the lobes. The homologous expansion does not start from zero at the centre: there is assumed a dense and constant-velocity inner region and the applied averaging is biased towards this low-velocity region. Encouragingly, the upper limits are always higher than the lower limits. The upper limits are in most cases close to the 1D Torun ages.

4.4. Evolution

thumbnail Fig. 13

Averaged expansion velocity vs. kinematical age for the pyCloudy modelled PNe. The filled symbols show the nebular ages (derived from radii and velocities shown in Fig. 12), and the open symbols the lobe ages (see Sect. 4.3). Circles show the bipolars and squares the post-bipolars.

Open with DEXTER

The aim of this work was to test for possible evolutionary sequences among our objects. Figure 13 shows the expansion velocity versus the kinematical age. The open symbols show the ages of the lobes, and the filled symbols the nebular ages. The latter is defined as above and may be considered as upper limits. All bipolar PNe occupy a narrow region in age in the range 1–2 kyr. Some of the lobes are a bit younger than 1 kyr. The nebular ages are 0.5–1 kyr higher than the lobes. For the two post-bipolars, the lobes are 2.5 kyr old; their nebular ages are around 4 kyr, a much larger difference than for the bipolars.

4.4.1. Six bipolar PNe

The interesting result is that the selected bipolar PNe have nearly the same age (see Fig. 13). However, their central star temperatures are spread over a wide range. This may indicate a range of stellar masses (higher mass stars evolve more quickly to higher temperatures). This result suggests that the lobes may form at a well-defined point during the nebular expansion, independent of stellar mass.

None of the selected bipolars are older than 1.52 kyr. A sample of six objects is too small to draw any general conclusions, but this may point to a fading time for the lobes. A short window of visibility for bipolar lobes would have consequences for their rate of occurrence. The visibility time of the nebulae within the size range selected for HST (5 arcsec or less) is 34 kyr (Gesicki et al. 2014, their Fig. 11) (assuming the accelerated Schönberner tracks). This is a general visibility time that does not depend on the nebular structure. Our sample of six represents 17% of the sample of compact PNe observed with HST. The lobes can be seen for about half of the total visibility period. Correcting for this, we find that 30%–35% of the compact bulge PNe may experience a bipolar phase. This number will be higher if we can identify somehow the bi-lobed PNe seen along the symmetry axis.

Our conclusion concerning ages is not unique among related objects. Huggins (2007) has shown that jets and tori in post-AGB stars have very similar ages; jets are younger by a few hundred years. The expansion times range from a few hundred to a few thousand years and the torus pre-dates the jet. This was also found for the torus and lobes of NGC 6302 (Szyszka et al. 2011) where the lobes are younger by ~600 yr from the 2900 yr old torus. The pre-PN M 1-92 (Alcolea et al. 2007) kinematical age of 1200 yr (the same in polar and equatorial directions of a bipolar molecular outflow) is also very similar to our data.

4.4.2. Two post-bipolar PNe

In the sample of 36 PNe we tried to identify possible descendants of the bipolar PNe. We searched amongst the oldest from the data of Gesicki et al. (2014). Two objects seem to be the likely candidates, namely H 2-15 and H 2-26. Both are older than 3.5 kyr (among the four oldest), both have very hot central stars (supporting their evolutionary advanced status), and both show morphologies resembling a torus, already well expanded with evidence of fragmentation. They also reveal interactions with the interstellar medium. Some faint structures might be the remnants of bipolar lobes. The derived velocity fields contain a linear part, consistent with residual lobe emission. However, more detailed observations would be needed to confirm this.

For both PNe the pyCloudy models confirmed their old kinematical age although the new ages are somewhat smaller. Both are larger than the classic bipolar PNe (Fig. 12) and approximately two times older (Fig. 13).

Their expansion velocities are low. The velocities may be biased due to the weakness of the lobes. However, it is not unlikely that these two are descendants of the subgroup of slow expanding PNe, while the objects that expand faster after a couple of thousand years evolve beyond the detection limit, but this could also be a result of the limited sample and we could just be seeing two of these objects.

Some hints concerning the evolved nature and the timescales of similarly looking PNe can be found in Manchado et al. (2015) where a sequence of hydrodynamical models has been presented, intended to explain the clumpy structure observed in PN NGC 2346. The published cross-sections through the nebular emission measures (unfortunately not the rendered images) show how the equatorial toroid already breaks up into clumps about 1000 yr after the decline of the fast stellar wind which earlier produced the bipolar lobes, which is, interestingly, similar to our post-bipolar objects.

It cannot be proven that the two PNe form late stages of bipolar evolution, but in view of the short visibility time of bipolar lobes, the presence of such objects in the bulge sample would be expected. As noted in the previous subsection the phase of symmetric lobes lasts about half of the total visibility time, so the other half can appear like these.

5. Conclusions

The steps we took in our research were the following:

  • We selected targets from the HST survey of 36 compact bulge PNe, not constrained on brightness or morphology. Six PNe in the sample have clear bi-lobed, axi-symmetric structure suitable for pseudo-3D modelling. Two other targets are possibly in a post-bipolar phase. Most of the targets are believed to be at the distance of the Bulge, at 8 kpc, and two of them are modelled at distances of 12 kpc and 4 kpc.

  • The pyCloudy models are able to fit the available constraints. The constant density per segment and the simple velocity field limit the accuracy, but also keep the number of free parameters reasonable.

  • The 3D models yield much lower ionized masses than the 1D models, typically reduced by a factor of 2–7. We find ionized masses in the range 0.02–0.1 M. This is an important result with implications for the evolution of planetary nebulae.

  • The spectra are successfully reproduced with a linear radius-velocity relation in the outer regions. This simple modelling recipe appears to be in agreement with the hydrodynamical calculations of Huarte-Espinosa et al. (2012) and with a very detailed analysis of the bipolar PN NGC 6302 by Szyszka et al. (2011). We caution that the current models did not explore a wide range of velocity fields, and the spectra were extracted at the central position of the nebula.

  • Averaged expansion velocities derived from the pyCloudy models are similar but not identical to the 1D models. Weighting the velocity field by the Hβ line emissivity is preferred to weighting by mass owing to the limitations on the density distribution in this work. The average expansion velocity ranges from 15 to 30 km s-1, with one object (H 2-18) expanding at 40 km s-1.

  • We derive ages for the lobes based on the model Hubble flow ranging from 500 to 1500 yr. The two post-bipolars show older lobe ages of 2500 yr. The lobe ages provide lower limits to the age of the nebula. We derive upper limits to this age from the average outer radius and expansion velocity: these range from 1300 to 2000 yr for the bipolars and 4000–4500 for the two post-bipolars.

  • The age range of the bipolars is remarkably narrow. The suggestion is that the bright lobes last a relatively short time and fade after ~1000 yr.

  • We cannot be confident that the two post-bipolar PNe actually are descendants of bi-lobed objects. However, the derived parameters do not contradict such a relation. A well-known hypothesis recently summarized by Frew & Parker (2012) states that “several stellar evolutionary pathways form objects classified as PNe”. Our morphological selection appears to have picked one such pathway, within which the clear bi-lobed structure can be seen at an early phase at ages of 1–2 kyr, while the post-bipolar phase can be identified at ages of about 3 kyr, all within the typical total visibility time of about 3–4 kyr. We still do not know whether this pathway concerns single or binary star evolution.

Appendix A: The photoionization model parameters and results

Table A.1

Photoionization model parameters and modelling results for the program nebulae.

Table A.2

Geometry defining parameters for the program nebulae.


3

pyCloudy returns the mass of ionized hydrogen. It should be multiplied by a factor of 1.4 to obtain the mass of the ionized region.

Acknowledgments

This publication was supported by the Polish National Science Centre (NCN) through the grant 2011/03/B/ST9/02552. We acknowledge the Mexican project that made pyCloudy possible: CONACyT grant CB-2010/153985.

References

All Tables

Table A.1

Photoionization model parameters and modelling results for the program nebulae.

Table A.2

Geometry defining parameters for the program nebulae.

All Figures

thumbnail Fig. 1

Hα HST images presented at the same angular scale. Each panel is 12 arcsec per side. The intensity scale is linear. The names are indicated in each panel. The first three columns show PNe in the bipolar phase while the rightmost column shows the proposed post-bipolar phase (see text). Six PNe are assumed to belong to the Galactic bulge while the smallest and largest, both shown in the third column, are assumed to be located behind and in front of the bulge, respectively.

Open with DEXTER
In the text
thumbnail Fig. 2

Emission line profiles for different inclination angle. Plotted is the [N ii] 6583 Å line computed for the model of the spatially unresolved nebula Hen 2-260. The solid line represents the inclination angle of (symmetry axis within the plane of the sky), long-dashed line: inclination of 30°; short-dashed line: 60°.

Open with DEXTER
In the text
thumbnail Fig. 3

Hen 2-262. The top left panel shows the assumed density structure only for the ionized nebula, drawn at the cross-section plane taken along the positioned horizontal symmetry axis. The magenta lines indicate the directions along which the 1D constant density Cloudy models were calculated. The end of each line segment corresponds to the assumed matter (density) boundary. The middle left panel shows the Hα image (i.e. a projection onto a plane of the 3D structure) obtained for the same model, but inclined to the sky plane at 20° to correspond approximately to the HST image shown in the bottom left panel (also in Fig. 1). The intensity grey scale is linear in all three images, the size is scaled to fill the width of the box. The red rectangle in the middle panel indicates the area over which the line profiles are integrated. The nebular emission line profiles are shown on the right; the solid black lines show the observations and the dashed red lines show the computed model for Hα, [O iii] 5007 Å, and [N ii] 6583 Å.

Open with DEXTER
In the text
thumbnail Fig. 4

H 2-20. Data presented as in Fig. 3. The axis inclination is 30°.

Open with DEXTER
In the text
thumbnail Fig. 5

H 2-18. Data presented as in Fig. 3. The axis inclination is 35°.

Open with DEXTER
In the text
thumbnail Fig. 6

Hen 2-260. Data presented as in Fig. 3. The axis inclination is 50°. The spectrograph slit (left middle panel) covers the whole nebula.

Open with DEXTER
In the text
thumbnail Fig. 7

Wray 16-286. Data presented as in Fig. 3. The axis inclination is 20°.

Open with DEXTER
In the text
thumbnail Fig. 8

H 2-27. Data presented as in Fig. 3. The axis inclination is 35°.

Open with DEXTER
In the text
thumbnail Fig. 9

H 2-15. Data presented as in Fig. 3. The axis inclination is 45°.

Open with DEXTER
In the text
thumbnail Fig. 10

H 2-26. Data presented as in Fig. 3. The axis inclination is 40°.

Open with DEXTER
In the text
thumbnail Fig. 11

Averaged expansion velocity from Torun models compared with pyCloudy. The abscissa shows the Torun models mass-averaged velocities, while at ordinate the pyCloudy Hβ-averaged values are shown as filled symbols, and the pyCloudy mass-averaged as open symbols. The (red) circles show the bipolar nebulae and the (blue) squares the post-bipolars. Data from Table A.1.

Open with DEXTER
In the text
thumbnail Fig. 12

Averaged expansion velocity vs. averaged outer radius for the pyCloudy modelled PNe. Velocity are Hβ-averaged and the outer radius is mass averaged (see text). Data from Table A.1. The filled circles show the bipolar nebulae and the filled squares the post-bipolars.

Open with DEXTER
In the text
thumbnail Fig. 13

Averaged expansion velocity vs. kinematical age for the pyCloudy modelled PNe. The filled symbols show the nebular ages (derived from radii and velocities shown in Fig. 12), and the open symbols the lobe ages (see Sect. 4.3). Circles show the bipolars and squares the post-bipolars.

Open with DEXTER
In the text

Current usage metrics show cumulative count of Article Views (full-text article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.

Data correspond to usage on the plateform after 2015. The current usage metrics is available 48-96 hours after online publication and is updated daily on week days.

Initial download of the metrics may take a while.