Morpho-kinematic properties of Wolf-Rayet planetary nebulae

The majority of planetary nebulae (PNe) show axisymmetric morphologies, whose causes are not well understood. In this work, we present spatially resolved kinematic observations of 14 Galactic PNe surrounding Wolf-Rayet ([WR]) and weak emission-line stars ($wels$) based on the H$\alpha$ and [N II] emission taken with the Wide Field Spectrograph on the Australian National University 2.3-m telescope. Velocity-resolved channel maps and position--velocity diagrams, together with archival Hubble Space Telescope ($HST$) and ground-based images, are employed to construct three-dimensional morpho-kinematic models of 12 objects using the program SHAPE. Our results indicate that these 12 PNe mostly have elliptical morphologies with either open or closed outer ends. The kinematic maps show the on-sky orientations of the interior shells in NGC6578 and NGC6629, as well as the compact ($\leq 6$ arcsec) PNe Pe1-1, M3-15, M1-25, Hen2-142, and NGC6567, in agreement with the elliptically symmetric morphologies seen in high-resolution $HST$ images. Point-symmetric knots in Hb4 exhibit deceleration with distance from the central star, which could be due to shock collisions with the ambient medium. The velocity dispersion maps of Pe1-1 also disclose the shock interaction between its collimated outflows and the interstellar medium. Collimated bipolar outflows are also visible in the position--velocity diagrams of M3-30, M1-32, and M3-15, which are reconstructed by tenuous prolate ellipsoids extending upward from dense equatorial regions in the kinematic models. The formation of aspherical morphologies and collimated outflows in these PNe could be related to the stellar evolution of hydrogen-deficient [WR] and $wels$ nuclei, which require further investigation.


Introduction
Planetary nebulae (PNe) are bright hydrogen-rich envelopes expelled by evolved stars with low-to-intermediate progenitor masses (1-8 M ⊙ ) in the asymptotic giant branch (AGB) phase. These envelopes are fully ionized by ultraviolet (UV) radiation from their hot degenerate nuclei in the post-AGB stage. Photoionization produces visible nebulae and also contributes to some thermal effects. Stellar winds from central stars of planetary nebulae (CSPNe) departing from the AGB phase introduce some hydrodynamic effects to these ionized envelopes that successively change their gaseous structures. The line emission produced by the ionized envelope allows us to observe its morpho-kinematic structure. Kinematic observations of PNe not only reveal their morphologies, but also provide insights into mass-loss processes in the AGB phase and during their transition from the AGB to PN, as well as PN evolution in the post-AGB phase (see e.g. Balick 1987;Corradi & Schwarz 1995;Balick & Frank 2002;Schönberner et al. 2005aSchönberner et al. ,b, 2010Kwok 2010).
Observations of PNe have mostly displayed elliptical and bipolar axisymmetric morphologies (see e.g. Balick 1987;Balick et al. 1987;Schwarz et al. 1992;Sahai et al. 2011;Weidmann et al. 2016), which pose a major challenge to theories describing their formation (see the review by Balick & Frank 2002). Although single stars are expected to produce spherical nebula shells, the interaction with the interstellar medium (ISM) may also result in elliptical shells. However, the theories face considerable difficulties in generating highly axisymmetric morphologies similar to what are seen in observations. According to the interacting stellar wind (ISW) theory (Kwok et al. 1978), a compressed, dense shell is formed when a slow (∼ 10 km s −1 ), dense wind from the AGB phase is pushed away by a fast (∼ 1000 km s −1 ), tenuous stellar wind from the PN phase. The ISW model extended by Kahn & West (1985), known as the generalized interacting stellar wind (GISW) theory, can also produce a highly axisymmetric gas distribution. Nevertheless, the GISW model is unable to fully describe all observed elliptical and bipolar morphologies (Balick & Frank 2002). Moreover, the GISW model needs a density contrast to produce an aspherical morphology. Recently, many complex axisymmetric morphologies have been found by high-resolution imaging observations (see e.g. Sahai et al. 2011), which cannot adequately be described by the GISW theory.
It has been suggested that axisymmetric morphologies are produced by tidal interactions with a binary companion (Soker & Harpaz 1992;Soker & Livio 1994;Soker 2006;Nordhaus & Blackman 2006;Nordhaus et al. 2007). This binary companion is usually expected to be a white dwarf that accretes material from an AGB star, though it could also be an exoplanet (De Marco & Soker 2011;Hegazi et al. 2020; Rapoport et al. 2021). Paczynski (1976) first proposed the role of binarity in shaping PNe through a common-envelope (CE) phase. Alternatively, it has been proposed that a mixture of strong toroidal magnetic fields and rotating stellar winds from a single star can lead to an equatorial density enhancement and jet-like outflows (García-Segura 1997;García-Segura et al. 1999;García-Segura & López 2000;Frank & Blackman 2004). On the other hand, Soker (2006) contended that the magnetic fields measured in PNe cannot play a central role in shaping PNe and that a single star cannot provide enough energy and angular momentum to support complex axisymmetric morphologies. Currently, there is a strong argument that most axisymmetric morphologies of PNe have been shaped through a binary channel (Miszalski et al. 2009a,b;De Marco 2009;Nordhaus et al. 2010). Furthermore, Miszalski et al. (2009b) discovered that nearly 30 percent, possibly as much as 60 percent, of bipolar PNe contain post-CE binaries, implying that the CE phase preferentially shapes axisymmetric PNe. The transformation of the orbital angular momentum of a binary system into the CE unbinds it and shapes an axisymmetric nebula, whose symmetric axis is perpendicular to the orbital plane of the binary system. Highresolution kinematic analyses of some PNe around post-CE binary stars have disclosed alignments between the nebular symmetric axes and binary orbital axes (see e.g. Mitchell et al. 2007;Jones et al. 2010Jones et al. , 2012Tyndall et al. 2012;Huckvale et al. 2013).
Small-scale low-ionization structures (LISs) inside or outside the main shells have been identified in nearly 10% of the Galactic PNe (Corradi et al. 1996;Gonçalves et al. 2001Gonçalves et al. , 2009, which are more visible in low-excited emission like [N II] and [S II]. They were classified as knots, jets, and jetlike systems by Gonçalves et al. (2001). Knots, either in pairs or isolated, are defined as those LISs with an aspect length-towidth ratio close to 1, while those with an aspect ratio much larger than 1 are classified as filaments. Jets are defined as highly collimated filaments that symmetrically appear as a pair on both sides of the central star and move with velocities much larger than the expansion velocity of the main shell. Those filaments with no evidence of velocities higher than the expansion velocity of the main body are called jetlike. However, projection effects make it difficult to distinguish between jets and jetlike systems. Highly collimated highvelocity jets or high-velocity pairs of knots constitute around half of the LISs, which are named fast, low-ionization emission regions (FLIERs; Balick et al. 1993Balick et al. , 1994. FLIERs have radial velocities of 25-200 km s −1 with respect to the main bodies (Balick et al. 1994). The formation of FLIERs, as well as their associated density and velocity contrasts with the main shells, are not yet well understood. It has been hypothesized that FLIERs are produced by axisymmetric superwind massloss through a CE channel, tidal interaction with a low-mass companion, and angular momentum deposition of a binary system (Soker 1990;Soker & Harpaz 1992).
Most CSPNe have hydrogen-rich stellar atmospheres, but a substantial fraction ( 25%) of them have been found to possess hydrogen-deficient fast expanding atmospheres characterized by high mass-loss rates (Tylenda et al. 1993;Leuenhagen et al. 1996;Leuenhagen & Hamann 1998;Acker & Neiner 2003). Their surface abundances exhibit helium, carbon, oxygen, and neon, which are products of the helium-burning phase and a post-helium flash (Werner & Herwig 2006 (Leuenhagen et al. 1996;Leuenhagen & Hamann 1998). A few CSPNe have narrower and weaker emission lines (C IV 5805Å and C III 5695Å) that are not identical to those of [WR] stars, the socalled weak emission-line stars (wels; Tylenda et al. 1993). These wels CSPNe are surprisingly higher towards the Galactic bulge and closer to the Galaxy's center, and they appear to have evolved from different stellar populations than [WR] stars (Górny et al. 2004(Górny et al. , 2009). However, medium-resolution spectroscopic studies of 19 stars, that had previously been classified as wels using low-resolution spectra, pointed to 12 normal H-rich, 2 probably H-deficient, and 5 unclassified stars (Weidmann et al. 2015).
This work is aimed at the morpho-kinematic properties of PNe surrounding [WR] and wels nuclei by means of optical integral field unit (IFU) spectroscopy. This paper is structured as follows. Section 2 describes our observational and data reduction processes. In Section 3, we present our IFU kinematic results. Three-dimensional (3D) morpho-kinematic models and their results are presented in Section 4, followed by conclusions and discussion in Section 5.

Observations
The IFU observations of the PNe studied in this work were collected with the Wide Field Spectrograph (WiFeS; Dopita et al. 2007Dopita et al. , 2010 at the Siding Spring Observatory in April 2010 under program number 1100147 (PI: Q. A. Parker). The WiFeS is an image-slicing IFU designed for the Australian National University (ANU) 2.3-m telescope that feeds a double-beam spectrograph. It samples 0.5 arcsec along each of twenty five 38 arcsec × 1 arcsec slitlets, resulting in a fieldof-view of 25 × 38 arcsec 2 with a spatial resolution element of 1 × 0.5 arcsec 2 . The spectroscopic output is projected onto a CCD detector of 4096 × 4096 pixels. Each slitlet is mapped onto 2 pixels on the CCD detector, which leads to a reconstructed point spread function (PSF) with a full width at half-maximum (FWHM) of ∼ 2 arcsec. The spectrograph uses volume phase holographic gratings to provide a spectral resolution of R ∼ 7000 or 3000, corresponding to a velocity resolution of either ∼ 21 or ∼ 60 km s −1 , respectively.
Our targets were observed with a spectral resolution of R ∼ 7000 covering the wavelength range of 4415-7070Å. As each spectrum is 4096 pixels long in the CCD detector, it provides a linear wavelength dispersion per pixel of 0.36Å in the blue arm and 0.45Å in the red arm with the R ∼ 7000 grating. The WiFeS spatial resolution with the reconstructed PSF FWHM of 2 arcsec enables us to identify morphologies of PNe with angular diameters larger than 6 arcsec. Although this spatial resolution is not ideal for resolving compact (≤ 6 arcsec) PNe, it can determine the on-sky orientations of morphologies from The mass-loss rates calculated using the formula given by Nugis & Lamers (2000).
b The evolutionary tracks of the helium-burning model by Blöcker (1995).
c Calculated from the nebular excitation class using the correlation by Dopita & Meatheringham (1990. their extended faint lobes or exterior halos (see Danehkar & Parker 2015). All targets were acquired in the classical data accumulation mode. We also collected bias frames, flat-field frames, twilight sky flats, arc lamp exposures, wire frames, and standard stars for bias reduction, flat-fielding, wavelength and spatial calibrations, and flux calibration. Table 1 presents a journal of the WiFeS observations. The usual PN names, PN in Galactic coordinate (PNG) numbers (Acker et al. 1992), coordinates (RA and Dec.), stellar types of the CSPN, as classified by Tylenda et al. (1993), Acker & Neiner (2003) and Depew et al. (2011), are given in Columns 1-4, respectively. Columns 5-6 provide the effective temperature (T eff ), stellar luminosity (L), and stellar wind terminal velocity (V ∞ ), along with their references, respectively. The mass-loss rates (Ṁ ) in Column 7 were calculated using the formula given by Nugis & Lamers (2000), with stellar luminosity and typical [WR] chemical composition of Y = 0.43 and Z = 0.56. Column 8 presents the stellar mass (M ⋆ ) estimated from the helium-burning evolutionary models of Blöcker (1995). The exposure times used for the observations and the observing date are given in Columns 9 and 10, respectively.
Our sample includes 10 well-known PNe with [WR] central stars chosen from the literature (Crowther et al. 1998;Acker & Neiner 2003). These PNe allow us to identify morphologies at various evolutionary stages with different stellar classes from early-to late-type [WR]. The morpho-kinematic structures of two PNe (Hen 3-1333 and Hen 2-113) with [WC 10] studied by Danehkar & Parker (2015) are very compact, since they are too young in comparison to the [WCE] PNe. The CSPN M 1-32 was classified by Acker & Neiner (2003) as the peculiar [WO 4] based on the wide FWHM of the C IV-5801/12 doublet, which corresponds to a terminal stellar wind velocity of V ∞ ≃ 4900 km s −1 . This stellar velocity is higher than V ∞ ≃ 2000-3000 km s −1 as typically seen in the usual [WR] CSPNe. The PN Th 2-A, whose morphological features were disentangled by Danehkar (2015), is another object with [WO 3] pec classified by Weidmann et al. (2008).
The data reduction was accomplished using a variety of techniques with the IRAF pipeline wifes. 1 We performed the flat-fielding using the frames taken with exposures from a quartz iodine (QI) lamp. A medium-averaged bias was subtracted from the calibration frames. We carried out the wavelength calibration with Cu-Ar arc exposures taken at the beginning of the night. The spatial calibration was conducted by employing wire frames collected from the diffuse illumination of the coronagraphic aperture with a QI lamp. Cosmic rays and bad pixels were removed from the raw data before the sky subtraction. We employed the IRAF task lacos im (LA-Cosmic package; van Dokkum 2001) to eliminate cosmic rays from the raw data, while bad pixels and any remaining cosmic rays were manually cleaned off afterward with the IRAF/STSDAS task imedit. For the purpose of sky subtraction, we chose a suitable sky window from the science data. We calibrated the science data to absolute flux units using observations of the spectrophotometric standard stars EG 274 and LTT 3864. The flux-calibrated data were corrected for the atmospheric extinction. We also rebinned the spaxel map by a factor of 2 and then performed a three-point median smoothing to reduce the sampling artifacts created by the image-slicing IFU. The high-resolution Hubble Space Telescope (HST) images of the PNe were retrieved from the Mikulski Archive for Space Telescopes (doi:10.17909/t9-pdkg-tg57), and the 3.5-m ESO New Technology Telescope (NTT) narrow-band images from Schwarz et al. (1992), which were utilized to constrain our morpho-kinematic models. Table 2 lists the HST images, together with their program IDs, instruments, filters/gratings, wavelength bands, spatial resolutions, and exposure times. The IRAF task lacos im was utilized to eliminate cosmic rays from the HST images. Table 3 lists the kinematic features measured from the integrated emission-line profiles in the spectroscopic observations of our sample. Column 3 presents the systemic velocity (v sys ) in the frame of the local standard of rest (LSR) based on the Hα emission that approximately represents the systemic velocity of the whole nebula, along with the previous value from Durand et al. (1998) in Column 4. We transferred the radial velocity to the LSR frame using the IRAF/ASTUTIL task rvcorrect. The mass-weighted mean expansion velocities (v HWHM ) associated with the half width at half maximum (HWHM) of Hα λ6563, [N II] λλ6548,6584 and [S II] λλ6716,6731, and their average velocity are presented in columns 5-8, respectively. The HWHM expansion velocities, V HWHM = 8 ln(2)σ true /2, were corrected for the instrumental width and the thermal broadening as follows: σ true = σ 2 obs − σ 2 ins − σ 2 th − σ 2 fs , where the instrumental width σ ins is measured from the [O I] λ5577 and λ6300 night sky lines having a typical value of ≈ 18 km s −1 for the WiFeS (R ∼ 7000), the thermal broadening σ th is estimated from the Boltzmann's equation

Mass-weighted Mean Expansion and Systemic Velocities
(Z is the atomic weight of the atom or ion), and the fine structure broadening σ fs is typically ≈ 3 km s −1 for Hα (Clegg et al. 1999). For comparison, Column 9 lists the expansion velocity, together with its reference, from the literature.
The radiation-hydrodynamic simulations conducted by Schönberner et al. (2010) demonstrated that the HWHM velocities of volume-integrated line profiles always underestimate the true expansion velocity, so the HWHM method is suitable for slowly expanding objects, but cannot provide true expansion velocities of extended PNe. For the Magellanic Cloud PNe, Dopita et al. (1985Dopita et al. ( , 1988 assumed that v exp = 1.82 V HWHM , nearly twice the widespread use of the HWHM velocity. This is associated with the maximum gas velocity behind the outer shock, and it may not be related to the mean expansion velocity of the nebular shell. Here, we assume v exp = V HWHM , which may correspond to the mean expansion of a spherical gaseous shell. However, the velocity-resolved channel maps in Section 4 imply that the HWHM velocities of integrated emission-line profiles do not correctly describe the mean expansion velocity of the nebular shell due to aspherical morphologies and collimated bipolar outflows extending upward from the main shell.

Spatially-resolved Kinematic Maps
Figure 1 presents spatially resolved maps of the flux intensity, continuum, radial velocity, and velocity dispersion of the Hα and [N II] line emission derived from the WiFeS observations of all the PNe in our sample (the figures for all the objects can be found in the online version of the journal). This information is obtained by fitting Gaussian functions to spaxels of the IFU datacube using the IDL library MPFIT (Markwardt 2009) developed for the nonlinear least-squares minimization fitting. The emission line profile is resolved if its emission line width is wider than the instrumental width (σ ins , see § 3.1). The contour lines in each figure depict the boundaries at ∼ 10 percent of the mean surface brightness of each object in the Table 3. Systemic velocities in the LSR frame and mass-weighted mean expansion (HWHM) velocities measured from the integrated emission-line profiles.   Figure 2 shows the velocity-resolved channels of the Hα and [N II] flux intensity, on a logarithmic scale, observed in a sequence of 12 or 18 velocity channels with a velocity resolution of ∼ 21 km s −1 relative to the LSR systemic velocity (the corresponding velocity channel maps for all the objects are available in the online journal). The stellar continuum map derived for each line profile was subtracted from the associated flux intensity maps. Figure 3 also presents position-velocity (P-V) diagrams of the Hα and [N II] line emission of each object extracted from the IFU datacube for two slits passing through the central star that are oriented with the position angles (PA) along and vertical to the symmetric axis of the best-fitting morpho-kinematic model constructed in § 4 (the corresponding P-V diagrams for all the PNe are available in the online version of the journal). The velocity axes on the P-V arrays are relative to the systemic velocity in the LSR frame. The stellar continuum of each line profile was also subtracted from the associated P-V diagrams. In § 4, we use the velocity channels, along with the P-V arrays, to constrain our morphokinematic models.
Our IFU velocity maps of PB 6 taken with a long exposure time (1200 sec) reveal that two groups of ionized gas are moving in opposite directions on both sides of the central star. The HST image of PB 6, taken with STIS/MIRVIS and a short exposure time of 22 sec on April 26, 2012 (Program ID: 12600, PI: R. Dufour), also depicts a complex morphology with many multi-scale filamentary structures and several knots, which are analogous to the morphology of NGC 5189, also around a [WO 1] star (Crowther et al. 1998). This object could contain low-ionization envelopes similar to those identified in NGC 5189 ). However, the short exposure time and the MIRVIS long-pass filter (4950-9450Å) of the HST image were inadequate for a detailed morphological analysis of PB 6.
The kinematics maps of M 3-30 likely suggest a pair of bipolar inner knots embedded in the elliptical shell. Previously, Stanghellini et al. (1993) also classified this object as an elliptical morphology with inner knots. However, our morphokinematic model (see § 4) constrained by the velocity channel maps and P-V diagrams indicates that they are associated with the outer ends of tenuous collimated bipolar outflows extending from the dense toroidal shell. This toroidal morphology is also visible in the narrow-band Hα+[N II] image taken with the 3.5 m ESO NTT (see Figure 4).
The IFU kinematic maps of Hb 4 illustrate a pair of pointsymmetric knots or collimated jets detached from the main shell, with a torus morphology, based on the HST imaging observations in Figure 4. These point-symmetric structures have extremely low brightness, but high velocity dispersion. The Corradi et al. (1996) also depicted the point-symmetric knots. Taking the inclination angle of i = −40 • found by the kinematic model (in § 4), the maximum values of the point-symmetric knots reach 160 ± 10 km s −1 with respect to the central star, which agrees with the long-slit kinematic analysis (López et al. 1997).
The IFU observations of IC 1297 likely imply an elliptical morphology, which is also supported by the Hα+[N II] image (Schwarz et al. 1992). However, Zuckerman & Aller (1986) described it as a broken ring, and Stanghellini et al. (1993) classified it under the irregular nebulae. Furthermore, Corradi et al. (1996) proposed the presence of a single isolated knot, which was later dismissed by Gonçalves et al. (2001), and turned out to be a field star. The Hα+[N II] and [O III] images taken by Corradi et al. (1996) again pointed to a broken ringshaped shell.  The velocity dispersion maps of Pe 1-1 depict two regions with high values of σ true = 35 ± 10 km s −1 on both sides, which are analogous to those in Hb 4. Its HST image in Figure 4 shows a barrel morphology having two collimated outflows in the regions with high velocity dispersion in the IFU maps. The on-sky orientation of the bipolar outflows determined from the IFU velocity maps is consistent with what is seen in the high-resolution HST image.
It can be seen that the IFU velocity and dispersion maps of M 1-32 are akin to what is found for Th 2-A, which is also viewed almost pole-on (Danehkar 2015). Similarly, the P-V arrays ( Figure 3) imply that this object also has collimated bipolar outflows with very low inclination (i = 4 • ; in § 4) relative to the line of sight. Peña et al. (2001)  The velocity IFU maps of K 2-16 also point to an elliptical nebula with upward bipolar extension. The narrow Hα+[N II] image of K 2-16 (Schwarz et al. 1992) exhibits a faint, circular nebula with a diameter of 23 ′′ . Kohoutek (1977) also suggested that it could be a large and old PN.
The velocity maps of Sa 3-107 likely hint at an axisymmetric morphology, though the detailed structures of this compact object cannot be resolved by WiFeS. Similarly, WiFeS observations of the compact PN PB 8 around a [WN/WC]hybrid star pointed to an aspherical morphology but without any precise details (Danehkar 2018). In particular, tenuous , and (n) Sa 3-107, followed by the associated synthetic velocity-resolved channel maps obtained for all the PNe, except for PB 6 and Sa 3-107, produced by the best-fitting morpho-kinematic models with the parameters given in Table 4. Each observed slice has a ∼ 21 km s −1 width, whose central velocity is given in km s −1 unit at the top of the panel. The LSR systemic velocity (vsys) of each object is given in km s −1 unit in the right bottom corner of each observed velocity channel map. The flux color in each observed slice is in logarithm of 10 −15 erg s −1 cm −2 spaxel −1 unit. The gray contour in each panel corresponds to ∼ 10 percent of the mean surface brightness of each object in the Hα emission (or R-band) retrieved from the SHS (or SSS). North is up and east is toward the left-hand side. The complete figure set (52 images) is available in the online journal. lobes extending from a compact nebula can spatially be mapped through IFU spectroscopy, which may help us determine the on-sky orientation of its axisymmetric morphology (e.g. Danehkar & Parker 2015). Thus, our WiFeS observations of Sa 3-107 could expose the on-sky orientation of its elliptical or bipolar morphology.

Morpho-kinematic modeling
To construct 3D morphological models based on velocityresolved channel maps and P-V diagrams, we have used the interactive kinematic modeling program SHAPE v5.0 (Steffen & López 2006;Steffen et al. 2011), which applies molded polygon meshes to 3D geometrical models of gaseous nebulae. It assembles a morphological model on a grid of volumetric cells and utilizes a ray-casting process to transfer radiation fields to these cells. The program produces outputs that can be directly compared with kinematic observations, namely the P-V diagram and velocity channel maps. P-V diagrams have been used for long-slit observations of many PNe (e.g. García-Díaz et al. 2009;Jones et al. 2010;López et al. 2012;Clark et al. 2013), whereas velocity channels have recently been employed to interpret kinematic IFU maps of some objects (Danehkar 2015;Danehkar et al. 2016). It also produces 2D rendered images of ionized nebulae that can be matched with high-resolution images. In particular, SHAPE can also export its 3D mesh model to the standard tessellation language (STL) file format, which can be converted to the extensible 3D graphics (X3D) format by mesh processing programs such as MeshLab (Cignoni et al. 2008) for inclusion in publications. This program has been used to create 3D morpho-kinematic models of many PNe, including NGC 6337 (García-Díaz et al. For this work, the velocity field is defined by a homologous expansion ) that increases linearly with distance from the nebular center. However, to reproduce the point-symmetric knots in the P-V diagrams of Hb 4, we use a decelerating velocity law that linearly decreases with distance from the nebular center. To model each object, we carefully choose a 3D geometry model and adjust various geometric parameters, along with the inclination and origination angles, in order to reproduce those kinematic features seen in the P-V arrays and velocity channels of each object, while the 2D rendered images are also expected to resemble the associated archival image. The inclination (i) relative to the line of sight, the position angle (PA) of the on-sky orientation, and geometric parameters are modified in an iterative process until the model velocity maps and P-V diagram closely resemble the kinematic observations. We were able to construct 3D morpho-kinematic models of the 12 PNe. For the compact PN Sa 3-107, we do not have any available high-resolution images. The HST/MIRVIS image of PB 6 is not suitable for morpho-kinematic modeling, but it may suggest a complex morphology with multi-scale features similar to NGC 5189 . To reproduce the WiFeS observations, we configured the seeing parameters in SHAPE according to the WiFeS spatial resolution that resulted in the convoluted P-V and channel outputs, while 2D rendered images were made  without the convolution processing for comparison with the high-resolution archival images. Figure 4 depicts the final SHAPE mesh models of the 12 objects viewed from the top (i = 0 • ) and front (i = 90 • ), as well as the best-fitting inclinations (i) and orientations (PA; see Table 4) before rendering and the rendered images, respectively (for all the objects available in the online journal). Additionally, an interactive X3D viewer of each SHAPE mesh model is provided in Figure 5 for the online journal. 2 In Figure 2, we also present the synthetic velocity-resolved channel maps created by the best-fitting kinematic models aimed at reproducing the velocity channel maps of the Hα and   Table 4 lists the key parameters of the best-fitting morphokinematic models. Columns 3-5 provide the position angle (PA) of the nebular symmetric axis projected onto the sky plane in the equatorial coordinate system, the Galactic position in the Galactic coordinate system, and the inclination (i) of the symmetric axis with respect to the line of sight, respectively. Column 6 presents the linear velocity factor (k shell ; km s −1 arcsec −1 ) describing the homologous expansion law v(km s −1 ) = k shell · r(arcsec) of the primary shell. Columns 7-8 give the mass-weighted average expansion velocity of the primary shell (v exp ) adopted according to the integrated flux HWHM measurements, and the polar expanding velocities of collimated outflows (v out ) with respect to the central star estimated from the P-V diagrams and velocity channels in Hα and [N II], respectively. For comparison, we also include the kinematic parameters of the PNe Th 2-A and M 2-42, which were determined by Danehkar (2015) and Danehkar et al. (2016), respectively. The SHAPE models of the 12 PNe studied here, besides Th 2-A and M 2-42 previously modeled, before rendering together with the rendered results at their best-fitting inclinations and orientations, are summarized in Figure 5, with fully interactive 3D mesh models in the X3DOM framework (Behr et al. 2009) available in the online version of this article.
The primary and secondary morphological features determined from our morpho-kinematic models are also summarized in Columns 10 and 11 of Table 4, respectively. Following Sahai et al. (2011), the classification codes include four primary classes: bipolar (B), elliptical (E), multipolar (M), and irregular (I). The B class is defined as objects with two primary, diametrically opposed lobes, centered on their expected location. The E class represents objects that are elliptical along a specific axis, and it also comprises ringshaped or torus morphologies. The M class defines objects showing more than one primary lobe pair whose axes are not aligned, such as NGC 5189 (see Danehkar et al. 2018). The I class is defined by objects do not display any geometrical symmetry. The structure may be open or closed at its outer ends, denoted by 'o' or 'c', respectively. Other morphological features described by secondary descriptors include point symmetry (ps), two or more pairs of diametrically opposed lobes (ps(m)), diametrically opposed ansae (ps(an)), the overall shape of the lobes (ps(s)), rings projected on lobes (rg), and the central region with a toroidal structure (t).
The SHAPE model of M 3-30 was constructed using a thick toroidal shell together with a thin prolate ellipsoid describing collimated outflows. For the thin prolate ellipsoid, we employed an inhomogeneous density model defined by a sinusoidal function in the spherical coordinate system, where the density profile decreases as it moves from the polar angle φ = 90 • (equator) toward 0 • and 180 • (polar) with a slightly higher density reduction in the outflow moving toward us. For the PN BD +30 • 3639 around a [WC9] star, which has P-V diagrams similar to M 3-30, Akras & Steffen (2012) similarly adopted a prolate ellipsoid with thinner and faster polar regions. To reproduce the [N II] observations, a different density law is applied to the toroidal shell, which has a powerlaw radial distribution with a power-law index of 3, similar to the [N II] model built for Abell 48 (Danehkar 2022). Moreover, the prolate ellipsoid is deactivated to generate the [N II] synthetic results. To improve the velocity channels, the thick toroidal shell was also modified by the shear operator in SHAPE. A sheared toroidal model was built in a similar way by Jones et al. (2010) to reproduce the P-V arrays of the PN Abell 41. A sheared elliptical shell could be a result of the PN-ISM interaction that occurs as the PN passes though the ISM. The rendered SHAPE representation of M 3-30 is akin to the narrow band Hα+[N II] image from Schwarz et al. (1992). The toroidal shell seen in the image is similar to the PN Abell 48 ), while M 3-30 also contains collimated outflows. The maximum velocity of the collimated outflows in the P-V diagram of the Hα emission corresponds to the outflow velocity of v out ≈ 80 km s −1 at the inclination of 25 • found by our model. The waist expansion velocity of the equatorial shell is about v exp = 36 ± 10 km s −1 similar to v HWHM = 36 ± 5 km s −1 estimated from the integrated emission-line profiles.
The HST observation of Hb 4 on a linear scale in Figure 4 shows that the main shell has a deformed ring-like structure, probably owing to the interaction with the ISM. Moreover, the HST observation shows that the point-symmetric knots or filaments are not vertically aligned, which again hints at the ISM interaction. Our IFU observations on a logarithmic scale indicate that the ring-shaped dense shell is also surrounded by a tenuous exterior halo. Accordingly, we used a SHAPE model that consisted of an equatorial dense torus, a tenuous oblate spheroid, and two sets of four point-symmetric knots   Figure 5 for an interactive viewer in the online journal), followed by the rendered image (fourth column) and the associated archival imaging observation (fifth column), respectively. The parameters of the best-fitting models are listed in Table 4. The archival data include the HST images of Hb 4, Pe 1-1, M 3-15, M 1-25, Hen 2-142, NGC 6578, NGC 6567 and NGC 6629, with the instruments and filters/gratings listed in Table 2; and the narrow-band Hα+[N II] filter images of M 3-30, IC 1297, M 1-32 and K 2-16 taken with the 3.5-m ESO NTT from Schwarz et al. (1992), with the image scales shown by solid lines. North is up and east is toward the left-hand side in each archival image. The complete figure set (12 images) is available in the online journal. using a string of four knots for each of them. Our P-V diagrams imply that the expansion velocity of the dense shell reaches 80 ± 20 km s −1 at the inclination of 40 • found by the morpho-kinematic model, while the mass-weighted average HWHM expansion velocity of v HWHM = 23 ± 1 km s −1 is estimated from the line profiles integrated over an aperture covering the main shell. The deprojected velocities of the bipolar collimated outflows are found to be in the range of v out = 30-160 km s −1 at i = 40 • ± 3 • . To match better the Hα P-V arrays, we also include a thin oblate ellipsoid around the dense toroidal shell, but it is deactivated in the [N II] model.
The P-V diagrams reveal the deceleration of the pointsymmetric knots in Hb 4, which were also identified by Derlopa et al. (2019). This deceleration is not the homologous expansion of typical PNe. We should note that hydrodynamic simulations of bipolar jets and ring-shaped nebulae (see e.g. García-Segura et al. 2001;Lee & Sahai 2004) and nebular shells (Schönberner et al. 2005a) indicate that the radial expansion velocity increases with the distance from the central star. Decelerating bow shocks have been discovered in the Herbig-Haro (HH) object HH34 (Devine et al. 1997) and post-AGB water-fountain star IRAS 18113-2503 (Orosz et al. 2019), which could be due to shock collisions and interactions with ambient media. Point-symmetric knots and LISs in other PNe also demonstrate some signs of shock excitation (Akras & Gonçalves 2016) and molecular H 2 emission (Akras et al. 2017(Akras et al. , 2020 because of the interaction with the ISM or their exterior nebula halos. To model the decelerating features of the point-symmetric structures, we employed two sets of four knots whose velocities linearly decrease with distance from the nebular center. This allowed us to reproduce the knot features in the P-V arrays and velocity channels. Although the knots adopted in our SHAPE model can reproduce our IFU observations, having a reconstructed PSF FWHM of ∼ 2 arcsec, we caution that the knots could be smaller, according to the high-resolution HST observations. The velocity dispersion maps of Hb 4 also depict higher values in the locations of these knots, which again suggest shock collisions with the surrounding medium. Based on the narrow-band Hα+[N II] image of IC 1297 (Schwarz et al. 1992;Corradi et al. 1996), a broken thick elliptic torus and a thin prolate spheroid are used to assemble a 3D model for this object. Additionally, another shifted, rotated tenuous prolate spheroid is included in the Hα kinematic model to recreate the high Hα expansion toward us projected mostly in the southern part. The velocity dispersion maps exhibit higher values in the southern part of the nebula, where a part of the torus is missing. High dispersion could be a sign of the ISM interaction or mass-accretion to a companion. Although the radial velocity map of IC 1297 is similar to that of M 3-30, the velocity channel slices and P-V arrays do not suggest the same kinematic structure. The Hα+[N II] image depicts a broken ring-shaped nebula with the same on-sky orientations seen in the WiFeS velocity maps. A maximum polar expansion (a) 3D mesh models (b) Rendered mesh models velocity of v out = 100 ± 10 km s −1 is determined from the [N II] P-V diagrams, while a mass-weighted expansion velocity of v exp = 31 ± 1 km s −1 is pointed out by the HWHM measurements.
The HST image of Pe 1-1 (Figure 4) allowed us to mold a 3D morphological model for this object. Moreover, the IFU velocity dispersion maps in Figure 1 manifest two regions with possible shock collisions on both sides of the nebula, similar to the point-symmetric knots seen in the dispersion maps of Hb 4. The high values of the velocity dispersion in Pe 1-1 could be due to the interaction of collimated outflows with the ambient medium. We therefore included a tenuous prolate ellipsoid associated with bipolar outflows protruding from both sides of a dense torus.
To model M 1-32, we exploited an equatorial thick torus resembling its Hα+[N II] image (see Figure 4), and a faint bi-conical shape reproducing collimated bipolar outflows reaching v out = 170 ± 20 km s −1 relative to the central star in the P-V arrays (see Figure 3). Additionally, a tenuous exterior halo was included to match the velocity channels, resulting in a model similar to Akras & López (2012). The velocity dispersion maps depict high values in the central region analogous to that of Th 2-A, which is similarly viewed almost pole-on (Danehkar 2015). Those high values in the dispersion maps are produced by high-velocity collimated bipolar outflows with a very low inclination (i = 4 • ). Akras & López (2012) also estimated v out = 200 km s −1 relative to the central star. The HWHM measurements of the integrated emission-line profiles also yield a mean expansion velocity of 31 ± 3 km s −1 .
The SHAPE model of M 3-15 was fabricated using a stretched, sheared torus based on the HST observation (Figure 4), and a thin prolate ellipsoid associated with collimated bipolar outflows with v out = 110 and 70 km s −1 according to the P-V diagrams (Figure 3) for Hα and [N II], respectively. Two dissimilar linear velocity factors are employed in the velocity law of the prolate ellipsoid to generate the synthetic Hα and [N II] results. The main shell in this PN could be deformed similar to M 3-30 and Abell 41 (Jones et al. 2010) owing to the ISM interaction. However, we caution that our WiFeS spatial resolution is insufficient for properly constraining the ring-shaped shell of this compact (≤ 6 arcsec) PN, so our 3D model is one of several possible solutions that can be achieved by modifying the main shell and its orientation. Akras & López (2012) also found outflow velocities of 100 km s −1 in a long-slit observation of M 3-15. Our inclination of 4 • with respect to the line of sight is also in agreement with Akras & López (2012). The HST image of  (2015) and Danehkar et al. (2016), respectively.
b The deprojected outflow velocities are not available for PB 6 and Sa 3-107.
Note. The morphological classification codes in Columns 10 and 11 are based on Sahai et al. (2011). The symbol '?' means that further observation is required to confirm the structure. The deprojected outflow velocities (vout) of the collimated outflows or polar expansions are with respect to the central stars. The mass-weighted average expansion velocities (vexp) are from the HWHM measurements of the integrated emission-line profiles (see Table 3). The linear factor k shell (km s −1 arcsec −1 ) of the homologous velocity law, v(km s −1 ) = k shell · r(arcsec), corresponds to the primary shell in each model. this PN looks very similar to the PN Vy 1-2 with [WR] or wels, which has a similar morphology seen almost pole-on with an inclination angle of 10 • (Akras et al. 2015). K 2-16's SHAPE model is made up of an elliptical shell and a thinner interior prolate spheroid with nonuniform density laws (defined with sin φ) decreasing from the equator (φ = 90 • ) toward the polar direction (φ = 0 • and 180 • ). A small interior sphere was also included to better match the P-V arrays. This model reproduces the P-V diagrams in the Hα emission (see Figure 3), while the lower part of the P-V diagram, extracted from the slit over the symmetric axis, was not covered by the WiFeS footprint. Taking the inclination of −30 • used by the 3D model, we derived a maximum polar expansion velocity of v out = 70 ± 10 km s −1 . The HWHM measurements also point to a mean expansion velocity of v exp = 31 ± 3 km s −1 .
The high-resolution HST images of M 1-25, Hen 2-142, NGC 6578, NGC 6567, and NGC 6629 disclose their elliptically symmetric morphologies. Although the spatial resolution of our IFU observations did not provide much detail of their elliptical morphologies, the IFU velocity maps display their on-sky orientations similar to what are seen in the HST images. Their 3D morpho-kinematic models (Figure 4) were constructed with modified prolate ellipsoids in a way where the rendered images resemble their inherent shapes, as seen in the HST images. We also included tenuous exterior halos in Hen 2-142, NGC 6578, and NGC 6629, which are visible in the HST observations. An exterior halo was utilized only for the Hα kinematic model of NGC 6567. Similarly, second larger exterior halos were added to the Hα models of NGC 6578 and NGC 6629, which are able to reproduce the Hα kinematic observations. The orientations and inclinations of these models were constrained by the P-V arrays ( Figure 3) and velocityresolved channels (Figure 2).

Nebular Evolutionary Ages
The expansion velocity v and radius r can be used to estimate the dynamical age t = r/v of each nebular component. However, a semi-empirical analysis by Dopita et al. (1996, D96) suggested that the actual evolutionary age should be 1.5 times greater than the dynamical age, i.e. τ = 1.5(r/v). Moreover, Gesicki & Zijlstra (2000, G00) assumed an AGB velocity of v init = 0.5v as the initial value, resulting in a mean explanation velocity of (v + v init )/2 and a true dynamical age of τ = 1.33(r/v).
The nebular evolutionary ages are presented in Table 5, including the optical angular dimensions a × b (Column 3), the distance (Column 4), the optical radius r opt = √ ab/2 (Column 5) associated with the optical angular dimensions at the distance D, and the radial distances r out of the collimated outflow or polar expansion from the central star based on the Hα and [N II] kinematic models (Columns 6 and 7). The mean actual dynamical ages τ mean of the main structure given in Columns 8 and 9 were calculated using the optical radius r opt and the mass-weighted average HWHM expansion v exp from Table 4 under the different assumptions, namely τ mean = 1.5 r opt /v exp (D96) and τ mean = 1.33 r opt /v exp (G00). Similarly, the actual dynamical ages τ out of the bipolar outflow or polar expansions (Columns 10-13) were determined from the radial distance r out and the deprojected outflow velocities v out (listed in Table 4) for the Hα and [N II] models. Note. The actual evolutionary ages estimated from vexp and vout as listed in Table 4 using the 1.5× dynamical age derived by Dopita et al. (1996, D96) and the initial expansion velocity of v init = 0.5v assumed by Gesicki & Zijlstra (2000, G00), discussed in text. A distance to M 2-42 was derived by Danehkar et al. (2016) using the Hα surface brightness-radius relation of Frew et al. (2016). We note that the optical radius r opt based on the optical angular dimensions (typically measured at 10-20% nebular brightness in the literature) may not be the same as the shell radius used in the SHAPE model, which was constructed in a way to replicate the observed P-V arrays and velocity channels. For example, the outer radius of 6 arcsec is employed for the primary shell in the 3D model of Hb 4, corresponding to r out = 0.148 pc, whereas r opt = 0.113 pc according to the optical angular dimensions. This is consistent with a radius of 0.8 times the outer radius assumed by Gesicki & Zijlstra (2000) in order to account for the velocity gradient in the age calculation using a mass-weighted average velocity. The adopted optical radius and mass-weighted HWHM expansion velocity estimated from the integrated Hα emission might yield the average age of the main nebula, while the collimated outflows and point-symmetric structure could be younger than the primary shell.
The true dynamical ages of the point-symmetric structures in Hb 4 were calculated using the maximum deprojected outflow velocity of v out = 160 km s −1 at the radial distance r out = 0.415 pc (16.8 arcsec) being projected onto 10.8 arcsec from the central star where the nearest knots are located. We obtain the mean evolutionary age of τ mean = 6420 yr (G00) at D = 5096 pc (Stanghellini & Haywood 2010), which is lower than the kinematic age of 3190 yr estimated by Derlopa et al. (2019) at D = 2.88 kpc. Moreover, Derlopa et al. (2019) derived a kinematic age of around 890 yr for the fastest knots (v out = 160 km s −1 ), whereas we get a true dynamical age of τ out = 3380 yr (G00). These discrepancies can be explained by dissimilar distances and the AGB velocity v init = 0.5v out assumed in our dynamical age calculation. Although the pointsymmetric knots in Hb 4 may not be as young as what was found by Derlopa et al. (2019), they could still be around 3000 yr younger than the primary shell.
It can be seen that the mean ages τ mean estimated from the mass-weighted average expansion velocities are not always the same as the outflow ages τ out estimated from the outflow or polar expansion velocities. Although these values are close in the PN Th 2-A, the outflow ages τ out are 1500 to 7200 years younger than the mean ages τ mean in other objects. In particular, the dynamical ages of the collimated outflows are about 500 yr in Pe 1-1 and Hen 2-142, which are lower than the mean ages of 2700 yr. Furthermore, the dynamical ages of the polar Hα expansion were found to be slightly below 1000 yr in M 1-25 and NGC 6567. This could be an indication of some recent outbursts from the post-AGB stars. However, the mean age of a nebula may be overestimated by using a mass-weighted average velocity, so the actual age could be lower.

Conclusions and Discussion
We have studied the morpho-kinematic features of a sample of 14 Galactic PNe surrounding [WR]-type and wels nuclei, using spatially resolved kinematic maps and position-velocity profiles attained from WiFeS integral field observations. The kinematic IFU data, together with archival high-resolution images, have been deployed to carry out 3D morpho-kinematic modeling of 12 PNe. Our results imply that these 12 objects mostly have elliptical axisymmetric morphologies, with either open or closed outer ends. Collimated outflows are found to be present in the PNe Hb 4, Pe 1-1, M 3-30, M 1-32, and M 3-15, while point-symmetric knots or jets are detached in Hb 4. Axisymmetric morphologies and collimated outflows have been recognized in other PNe around [WR] stars (Akras & López 2012). Some elliptical PNe with FLIERs were also observed to have "diffuse X-ray" sources (Kastner et al. 2012). Other aspherical PNe around [WR] stars, such as BD +30 • 3639 (Arnaud et al. 1996), NGC 40 (Montez et al. 2005, and NGC 5315 (Kastner et al. 2008), have been reported to contain diffuse X-ray sources and hot bubbles. The soft X-ray emission could be associated with collimated jets and wind-wind shocks (Kastner et al. 2012), whereas the hard X-ray emission may be an indication of the mass accretion into a binary companion.
The PNe in our sample have elliptically symmetric morphologies, apart from PB 6, IC 1297, and Sa 3-107, which could imply a possible link between the nebulae and their WR-type nuclei. PB 6 around a [WO 1] star likely has a complex morphology containing multi-scale structures similar to NGC 5189 with a [WO 1] star . IC 1297 also has an elliptical broken ring, while the missing ring part exhibits high velocity dispersion, which could be due to either the ISM interaction or massaccretion into a companion. Although we do not have any imaging observations of Sa 3-107, its IFU velocity maps suggest an axisymmetric morphology. The theoretical evolutionary models by Schönberner et al. (2005a,b) predict that the expansion velocity increases as the nebula evolves and the central star becomes hotter. However, we cannot distinguish any link between the nebular kinematics and stellar characteristics listed in Table 1. Some PNe with hot stars such as Pe 1-1 and Hb 4 have low mass-weighted average expansion velocities, while some PNe around cool stars such as K 2-16 demonstrate higher values. This tendency does not seem to be explained by single star evolution, since the mean nebular expansion velocities are not linked to their stellar parameters. Moreover, collimated bipolar outflows were found in some objects, Hb 4, M 3-30, M 1-32, and M 3-15, reaching outflow velocities in the rang of 60-200 km s −1 with respect to the central stars, which cannot be expected from the GISW model.
Recent hydrodynamic simulations demonstrate that binary stellar evolution can shape axisymmetric morphologies (e.g. Akashi et al. 2015;Chen et al. 2016;Akashi & Soker 2016;García-Segura et al. 2018, 2021Zou et al. 2020;Castellanos-Ramírez et al. 2021). It is worthwhile mentioning that magnetohydrodynamic simulations of the envelope of a single AGB star could not support bipolar PNe at physically possible rotating speeds (García-Segura et al. 2014). However, 3D hydrodynamic simulations of a binary system reproduced the wide-bipolar lobes of an AGB Mira-like variable that may be in transition to a PN (Chen et al. 2016). Recently, 2D hydrodynamic simulations of a system that has undergone a CE phase by García-Segura et al. (2018) were able to generate elliptical shapes similar to M 1-25, Hen 2-142, and NGC 6578, as well as ring-shaped elliptical morphologies with a pair of external bipolar point-symmetric structures, similar to Hb 4 and M 2-42. Additionally, 3D hydrodynamic CE simulations illustrate how the interaction between CE ejecta and spherical wind from the post-AGB star results in the formation of highly collimated outflows (Zou et al. 2020). The implication of triple systems was also investigated by Akashi & Soker (2017) through 3D hydrodynamic simulations that generate highly asymmetrical filamentary structures similar to those seen in PB 6 and NGC 5189.
As confirmed by HST imaging observations, our WiFeS integral field observations reveal the on-sky orientations of the elliptical morphologies of the interior shells in the PNe NGC 6578 and NGC 6629, as well as the compact (≤ 6 arcsec) PNe Pe 1-1, M 3-15, M 1-25, Hen 2-142, and NGC 6567. Although the WiFeS spatial resolution is inadequate for unveiling morphological details of compact objects with angular diameters of less than 6 arcsec, it can help us unmask the on-sky orientations of overall morphologies, as previously shown by Danehkar & Parker (2015). This study demonstrates that the WiFeS has the capability to conduct an extensive IFU survey of numerous objects where the precise morphological details are unimportant.
Future high-resolution imaging and kinematic observations of other PNe around [WR] and wels nuclei will definitely shed light on the mechanism that shapes the axisymmetric morphologies in these objects. In-depth spectroscopy and lightcurve monitoring of their central stars will help us gain a better understanding of them. The presence of binary companions should be investigated as they can be responsible for axisymmetric morphologies. Nevertheless, an aspherical morphology can also be formed by an exoplanet whose existence is extremely difficult to verify. Further studies of [WR] CSPNe will unravel the puzzle of elliptical and bipolar morphologies in PNe.
This work is partially based on a dissertation by the author (Danehkar 2014), supported by an international Macquarie University Research Excellence Scholarship (iMQRES) and a Sigma Xi Grant-in-Aid of Research (GIAR). The author wishes to thank the referee for useful suggestions and comments, Quentin Parker for supporting the 2010 ANU observations, Milorad Stupar for his guidance on IRAF data reduction, Wolfgang Steffen for helpful discussions and guidance on SHAPE, Nico Koning for upgrading SHAPE and support, David Frew for helping with the observing proposal, Kyle DePew for carrying out the ANU observations in 2010, and the ANU telescope time allocation committee and the staff at the ANU Siding Spring Observatory for their support. Based on observations made with the NASA/ESA Hubble Space Telescope, obtained from the Data Archive at the Space Telescope Science Institute, which is operated by the Association of Universities for Research in Astronomy, Inc., under NASA contract NAS 5-26555.
Software: IDL Astronomy User's Library (Landsman 1993), IDL Coyote Library (Fanning 2011).    and (n) Sa 3-107, followed by the associated synthetic velocity-resolved channel maps obtained for all the PNe, except for PB 6 and Sa 3-107, produced by the best-fitting morpho-kinematic models with the parameters given in Table 4. Each observed slice has a ∼ 21 km s −1 width, whose central velocity is given in km s −1 unit at the top of the panel. The LSR systemic velocity (v sys ) of each object is given in km s −1 unit in the right bottom corner of each observed velocity channel map. The flux color in each observed slice is in logarithm of 10 −15 erg s −1 cm −2 spaxel −1 unit. The gray contour in each panel corresponds to ∼ 10 percent of the mean surface brightness of each object in the Hα emission (or R-band) retrieved from the SHS (or SSS). North is up and east is toward the left-hand side.  Table 4. Slits oriented with the position angles (PA) along and vertical to the symmetric axis of the morpho-kinematic model passing through the central star, respectively. The velocity in each observed P-V arrays is with respect to the systemic velocity of the object, given in km s −1 unit. The angular offset at 0 arcsec is the nebular center. , and the best-fitting inclination and orientation (third column; see Figure 5 for an interactive viewer in the online journal), followed by the rendered image (fourth column) and the associated archival imaging observation (fifth column), respectively. The parameters of the best-fitting models are listed in Table 4. The archival data include the HST images of Hb 4, Pe 1-1, M 3-15, M 1-25, Hen 2-142, NGC 6578, NGC 6567 and NGC 6629, with the instruments and filters/gratings listed in Table 2; and the narrow-band Hα+[N II] filter images of M 3-30, IC 1297, M 1-32 and K 2-16 taken with the 3.5-m ESO NTT from Schwarz et al. (1992), with the image scales shown by solid lines. North is up and east is toward the left-hand side in each archival image. , and the best-fitting inclination and orientation (third column; see Figure 5 for an interactive viewer in the online journal), followed by the rendered image (fourth column) and the associated archival imaging observation (fifth column), respectively. The parameters of the best-fitting models are listed in Table 4. The archival data include the HST images of Hb 4, Pe 1-1, M 3-15, M 1-25, Hen 2-142, NGC 6578, NGC 6567 and NGC 6629, with the instruments and filters/gratings listed in Table 2; and the narrow-band Hα+[N II] filter images of M 3-30, IC 1297, M 1-32 and K 2-16 taken with the 3.5-m ESO NTT from Schwarz et al. (1992), with the image scales shown by solid lines. North is up and east is toward the left-hand side in each archival image.