Dusty photoionization model of NGC 6781

The Herschel Planetary Nebula Survey (HerPlaNS) - a comprehensive dusty photoionization model of NGC6781


We perform a comprehensive analysis of the planetary nebula (PN) NGC 6781 to investigate the physical conditions of each of its ionized, atomic, and molecular gas and dust components and the object’s evolution, based on panchromatic observational data ranging from UV to radio. Empirical nebular elemental abundances, compared with theoretical predictions via nucleosynthesis models of asymptotic giant branch (AGB) stars, indicate that the progenitor is a solar-metallicity, initial-mass star. We derive the best-fit distance of 0.46 kpc by fitting the stellar luminosity (as a function of the distance and effective temperature of the central star) with the adopted post-AGB evolutionary tracks. Our excitation energy diagram analysis indicate high excitation temperatures in the photodissociation region (PDR) beyond the ionized part of the nebula, suggesting extra heating by shock interactions between the slow AGB wind and the fast PN wind. Through iterative fitting using the Cloudy code with empirically-derived constraints, we find the best-fit dusty photoionization model of the object that would inclusively reproduce all of the adopted panchromatic observational data. The estimated total gas mass () corresponds to the mass ejected during the last AGB thermal pulse event predicted for a initial-mass star. A significant fraction of the total mass (about 70 %) is found to exist in the PDR, demonstrating the critical importance of the PDR in PNe that are generally recognized as the hallmark of ionized/H regions.

ISM: planetary nebulae: individual (NGC 6781) — ISM: abundances — ISM: dust, extinction

M. Otsuka

1 Introduction

Photometry Observations
Obs-Date Telescope Instrument Band Aperture (Nebula+CSPN) Program-ID/PI References
2011-07-25 GALEX GALEX NUV 180″
2008-07-31 ING/INT 2.5-m WFC RGO , Sloan and 320″ I08AN02/P. Groot
2015-05-12 ESO/NTT 3.6-m EFOSC2 Bessel , , 200″ 60.A-9700(D)/Calibration
2009-08-09 ING/INT 2.5-m WFC IPHAS H 320″ C129/J. Casare
1995-07-24 HST WFPC2/PC F555W, F814W (CSPN only) GO6119/H.E. Bond
2010-06-26 UKIRT 3.8-m WFCAM , , 180″
2010-04-13 WISE WISE 3.4, 11.6, 22.1 µm
2004-04-20 Spitzer IRAC 4.5, 5.8, 8.0 µm 240″ 68/G. Fazio
1996-04-28 ISO ISOCAM 14.3 µm 240″ COX 1/P. Cox
2011-10-17 Herschel PACS 70, 100, 160 µm 240″ OT1-tueta-2/T. Ueta (1)
2011-10-11 Herschel SPIRE 250, 350, 500 µm 240″ OT1-tueta-2/T. Ueta (1)
Radio telescopes Various 1.4, 5, 22, 30, 43 GHz (2), (3), (4), (5), (6)
Spectroscopy Observations
Obs-Date Telescope Instrument Wavelength Program-ID/PI References
1997-08-09 ING/WHT 4.2-m ISIS 3600-8010 Å W-97B-41/X.-W. Liu (7), (8)
2005-10-19 Spitzer IRS 5.2-39.9 µm 1425/IRS-Calibration
2011-10-14 Herschel PACS 51-220 µm OT1-tueta-2/T. Ueta (1)
2012-04-01 Herschel SPIRE 194-672 µm OT1-tueta-2/T. Ueta (1)

References. – (1) HerPlaNS1; (2) Condon et al. (1998, 376.5 12 mJy at 1.4 GHz); (3) Stanghellini & Haywood (2010, 323 mJy at 5 GHz); (4) Petrov et al. (2007, 190 mJy at 22 GHz); (5) Pazderska et al. (2009, 264.1 7.1 mJy at 30 GHz); (6) Umana et al. (2008, 710 mJy at 43 GHz); (7) Liu et al. (2004a); (8) Liu et al. (2004b).

Table 1: The log of panchromatic observations of NGC 6781 adopted for the present study.

The life cycle of matter in the Universe is intimately connected with the stellar evolution because stars are the most fundamental building blocks of the Universe. Hence, the chemical evolution of galaxies has always been made possible by stellar nucleosynthesis, convection/dredge-up, and ultimately, stellar mass loss. This stellar mass loss becomes significant when stars evolve into the final stage of stellar evolution, i.e., the asymptotic giant branch (AGB) stage for low-mass stars () and core-collapsed supernovae explosions for high-mass stars ().

Either way, the mass loss process would expel a significant fraction of mass contained in stars as the circumstellar shells, which would eventually become part of the interstellar medium (ISM). Besides gas, molecules and solid state particles (i.e., dust grains) participate in the stellar mass loss and make up a significant part of the circumstellar shells as the photodissociation region (PDR). These cold components of the mass loss ejecta will provide the seed material for the formation of the next generation of stars and planets. Hence, understanding of stellar mass loss is important in characterizing the cosmic mass recycling and chemical evolution in galaxies.

Planetary nebulae (PNe) are low-mass stars that have completed mass-loss during the preceding AGB phase and consist of a hot central star ( K; evolving to become a white dwarf) and an extensive circumstellar shell. While PNe are famous for their spectacular circumstellar structures seen via bright optical emission lines arising from the ionized gas component of the nebula, the ionized part of PNe is surrounded by the neutral gas and dust components (i.e., the PDR). Therefore, being relatively isolated from surrounding objects, PNe provide unique laboratories to further our understanding of the stellar evolution and the chemical evolution of galaxies, from high-temperature fully-ionized plasma to low-temperature dusty molecular gas.

So far, more than 2000 PNe in the Milky Way have been identified (Frew, 2008; Parker et al., 2016). The evolutionary history of the progenitor (the central star of a PN, CSPN) is imprinted in the circumstellar shells. Radiation from the CSPN permeates into the circumstellar shells, controlling the physical conditions and local structures (see e.g., Villaver et al., 2002). Moreover, PNe are in the evolutionary stage in which the circumstellar shells would reach their largest extent before the material at the periphery begins to dissipate into the ISM. Therefore, by investigating spatially-extended emission from each of the ionized, atomic, and molecular gas and dust components, one can infer ionic, elemental, and molecular/dust abundances and the mass-loss and evolutionary histories of the CSPN.

Because PNe are H regions, there is a history of observations that has generated a wealth of archival data in the UV and optical. Similarly, the bright ionized gas in PNe is also bright in the radio continuum. With the advent of new technologies, PN observations in the X-ray, near-IR, and mid-IR follow suit. Recently, a window of opportunity in the far-IR was brought forward by a suite of space telescopes, which filled the remaining hole in the spectral coverage. We seized this opportunity and initiated the Herschel Planetary Nebula Survey (HerPlaNS; Ueta et al., 2014, HerPlaNS1, hereafter) and its follow-up archival study, HerPlaNS, using data collected for a hoard of PNe with the Herschel Space Observatory (Pilbratt et al., 2010, Herschel, hereafter).

In our previous pilot/demonstration study, we focused on the bipolar PN NGC 6781 to empirically characterize its dusty circumstellar nebula based mainly on far-IR data. We confirmed a nearly pole-on barrel structure of the dust shell (of  K, ) rich in amorphous carbon via broadband mapping. We also determined the physical stratification of the nebular gas (of 0.86 ) in terms of the electron density and temperature via spatially-resolved far-IR line diagnostics. Moreover, we yielded a gas-to-dust mass ratio map by a direct comparison between the empirically-derived dust and gas distributions. These analyses were made with the adopted distance of 0.95 kpc. Assuming that all mass-loss ejecta were detected and that the present-day core mass were 0.6 , we concluded that a initial-mass progenitor was about to complete its PN evolution.

In the present study reported here, we continue our investigation of NGC 6781 by adopting as much panchromatic data as possible in addition to our own HerPlaNS far-IR data. This time, our focus is to generate a coherent model of NGC 6781 that would satisfy the adopted panchromatic data as comprehensively as possible. To this end, we first derive the empirical characteristics of the central star and its circumstellar nebula with a greater amount of self-consistently based on the adopted panchromatic data set. Then, we use these empirically derived quantities as more constraining input parameters for a dusty photoionization model consisting of ionized, atomic, and molecular gas plus dust grains to construct one of the most comprehensive models of the object ever produced. In doing so, preference is given to adopting a panchromatic data set rather than exploiting the spatially-resolved nature of the data. This is also because, while the existing multi-band images of the nebula certainly help us to empirically establish its 3-D structures, the amount of imaging data (especially emission line maps) still lacks to fit detailed 3-D models of internal stratifications in the nebula.

The organization of the rest of the paper is as follows. We summarize the panchromatic observational data of NGC 6781 adopted in the present study (§ 2) and review each of the ionized, atomic, and molecular gas and dust components of the nebula and the central star to derive empirical parameters that are pertinent to the subsequent dusty photoionization model fitting (§ 3). Then, we present the best-fit dusty photoionization model of NGC 6781 produced with Cloudy (Ferland et al., 2013) while emphasizing improvements made by adopting the panchromatic data comprehensively and self-consistently (§ 4), before describing conclusions drawn from the empirical analyses and modeling (§ 5). This study would demonstrate that the derived best-fit model is robust enough to empirically constrain theoretical stellar evolutionary predictions and that the cold dusty PDR of PNe is at least equally important as the ionized part when characterizing their progenitor’s evolution and mass-loss histories, especially in the context of the cosmic mass recycling and chemical evolution of galaxies.

2 Adopted panchromatic data of NGC 6781

2.1 Photometry data

We collect photometry data – 10 and 27 data points for the CSPN alone and the nebula plus the CSPN, respectively – from previous observations made with various ground- and space-based telescopes as listed in Table 1 and plotted in Fig. 1. We re-reduce the archived data ourselves to perform photometry measurements unless science grade images are already made available. The diameter of the adopted photometry aperture for the entire nebula (including the CSPN) is indicated in Table 1. For photometry of the CSPN alone, we use a circular aperture of 0.4 (HST), 1.2 (EFOSC2), 3.8 (WFC), and 2.2 (WFCAM) centered at the CSPN. In Appendix A, we outline the method of data reduction and photometry for the HST/WFPC2, INT 2.5-m/WFC, ESO NTT 3.6-m/EFOSC2, UKIRT 3.8-m/WFCAM, and INT 2.5-m/IPHAS H broadband images.

Figure 1: The panchromatic photometric and spectroscopic data of NGC 6781 adopted in the present study. Broadband photometry was done over the entire extent of the nebula from the following sources: GALEX (open triangle), ING/INT (open circles), ESO/NTT (pluses), UKIRT (crosses), WISE (asterisks), Spitzer (filled circles), ISO (filled square), Herschel (open squares), and Radio (filled triangles), while photometry of the CSPN (filled stars) was also done using HST/WFPC2 images in addition to the above optical and near-infrared sources. Spectra (grey lines) were sourced from WHT/ISIS, Spitzer/IRS, and Herschel/PACS and SPIRE. The adopted spectra from four instruments are shown in grey lines, with their respective spectral ranges indicated at the bottom. The inset displays the Spitzer/IRS spectra in the mid-IR full of H, polycyclic aromatic hydrocarbons (PAHs) and ionized gas emission features/lines with the dust continuum steadily rising toward longer wavelengths from around 10. See text for how the data were scaled with respect to each other. See also Tables A1 and B1.

2.2 Spectroscopy data

We collected spectroscopy data from previous optical, mid-IR, and far-IR observations made as summarized in Table 1 and plotted in Fig. 1. Detailed accounts of data reduction and spectroscopic measurements are given in Appendix B for each instrument (WHT/ISIS in the optical, Spitzer/IRS in the mid-IR, and Herschel/PACS and SPIRE in the far-IR). Also given in Appendix B is a detailed description as to how the H flux of the entire nebula is estimated using the IPHAS H image. Our choice of the data sets is motivated to ensure that the adopted spectra represent the bulk of the nebula. Fig. 2 shows relative slit positions with respect to the entire nebula.

Figure 2: Relative slit positions of previous spectroscopic observations with respect to the NGC 6781 nebula shown on the NOT/H image (previously presented by Phillips et al. (2011)), in which field stars are subtracted by PSF-fitting. N is up and E is to the left.

Optical WHT/ISIS spectrum

The optical WHT/ISIS spectrum is obtained by scanning the nebula along declination during integrations, with the position angle (P.A., defined to be degrees E of N) of the slit is set at 90: the resulting spectrum, therefore, represents an average of the bulk of the central part of the nebula (X.-W. Liu, private communication). Fig. 2 shows the central scanned region of with a blue box. Flux densities of the WHT/ISIS spectrum are scaled to match the INT/WFC IPHAS H band fluxes (see Appendix B.2).

Mid-IR Spitzer/IRS spectra

The archival Spitzer/IRS (Houck et al., 2004) spectra are obtained with the SL (5.2-14.5 µm; a pair of the vertical light blue slits at P.A. of in Fig. 2) and LL (13.9-39.9 µm; the horizontal slit at P.A. of in Fig. 2) modules. Only the SL spectrum was previously presented (Phillips et al. 2011; Mata et al. 2016), whereas we include the LL spectrum in our analysis. While there is only little flux density offset between the SL and LL spectra, we combine the two spectra by scaling the SL spectrum to match the LL spectrum so that the combined mid-IR spectrum would represent the central part of the nebula (Fig. 1). Flux densities of the combined mid-IR spectrum are then scaled using the results of mid-IR photometry (see Appendix B.3).

Far-IR Herschel/PACS and SPIRE spectra

Far-IR Herschel spectra of the nebula for the present study are adopted from those previously presented (HerPlaNS1). To define a far-IR spectrum representing the bulk of the nebula we combine spectra from all PACS IFU spaxels ( in the apertures), while a single spectrum from the central bolometer of the SPIRE array is included (of and diameter in the short and long wavelength band, respectively; at both the “center” and “rim” positions as depicted as white boxes and gray circles, respectively, in Fig. 2). The combined far-IR spectra are then scaled using the flux density ratios between far-IR lines and H for the entire nebula with the synthesized H map constructed from the H image (see Appendix B.4).

Interstellar reddening correction and Flux measurements

Once we reconstruct spectra in the optical, mid-IR, and far-IR to represent the bulk of the nebula, we measure line fluxes by Gaussian fitting. For the ISIS spectrum, the line fluxes are corrected for the interstellar extinction with the following formula:


where () is the de-reddened line flux, () is the observed line flux, (H) is the reddening coefficient at H, and () is the interstellar extinction function at computed by the reddening law of Cardelli et al. (1989) with .

We measure the reddening correction factor (H) by comparing the observed Balmer line ratios of H and H to H with the theoretical ratios given by Storey & Hummer (1995) for an electron temperature  K and an electron density  cm under the assumption that the nebula is optically thick to Ly (so called “Case B”; e.g. see Baker & Menzel 1938; also see § 3.1.1 for the bases of these and values). The measured (H) turns out to be from the (H)/(H) ratio and from the (H)/(H) ratio. Thus, we adopt (H) of , which is a weighted-mean of the above values. We do not correct for the interstellar extinction at longer wavelengths than -band because extinction would be negligible at those wavelengths. The final de-reddening line fluxes measured in the adopted spectra are listed in Table B1. The quoted fluxes are normalized with respect to (H) = 100.

While we adopt these reprocessed 1-D panchromatic spectra and duly-measured de-reddened line fluxes as representative of the bulk of the nebula, a word of caution appears appropriate at this point. As Fig. 2 shows, the spatial coverage of the nebula by various spectroscopic apertures is not complete and uniform. As would become apparent later from the model fitting (§ 4), there would be some inconsistencies in line emission strengths, especially in neutral and low-excitation lines such as [N i], [O i], and [S ii] (see § 3.1.2). This is primarily because the highest surface brightness regions (the E and W end of the central ring structure; Fig. 2) are missed in the optical data and may be less strongly weighted than they should be in the far-IR data. We would return to these issues when we discuss model fitting in § 4.

3 Anatomy of NGC 6781

3.1 The ionized/neutral gas component

Plasma diagnostics

We determine the and pairs for the ionized/neutral gas component2 of NGC 6781 for a few temperature/ionization regions based on various collisionally-excited lines (CELs) and recombination lines (RLs) detected in the adopted panchromatic spectra. In the present plasma diagnostics and the subsequent ionic abundance derivations, we adopt the effective recombination coefficients, transition probabilities, and effective collisional strengths listed in their Tables 7 and 11 of Otsuka et al. (2010), in which all the original references to all the atomic data are found. The diagnostic line ratios used in the present analysis and the resultant and values are summarized in Table 2.

Figure 3: The - diagram based on CEL diagnostic lines. The dashed and solid lines are the and indicators, respectively. The ID numbers indicate the corresponding line ratios listed in Table 2.
ID Ion -diagnostics Ratio Result (cm)
1 [O i] (63 µm)/(146 µm) 590
2 [S ii] (6716 Å)/(6731 Å)
3 [O ii] (3726 Å)/(3729 Å)
4 [N ii] (122 µm)/(205 µm)
5 [S iii] (18.7 µm)/(33.5 µm)
6 [Ne iii] (15.6 µm)/(36.0 µm)
7 [O iii] (4959 Å)/(88.3 µm)
ID Ion -diagnostics Ratio Result (K)
8 [S ii] (6716/31 Å)/(4069 Å)
9 [N ii] (6548/83 Å)/(5755 Å)
10 [N ii] (6548/83 Å)/
(122 µm+205 µm)
11 [O ii] (3726/29 Å)/(7320/30 Å)
12 [Ar iii] (7135 Å+7751 Å)/(9.0 µm)
13 [O iii] (4959 Å)/(4363 Å)
14 [Ne iii] (3868 Å+3967 Å)/(36.0 µm)
He i (7281 Å)/(6678 Å)
Table 2: Summary of plasma diagnostics using nebular lines.

The - plot shown in Fig. 3 summarizes how and relate to each other in the regions of the nebula, from which the particular CELs involved in the diagnostic line ratios would arise: the solid lines are the - curves derived from the -sensitive ratios, while the dashed lines are those from the sensitive-line ratios. Strictly speaking, the diagnostic lines labeled as (1), (7), (8), (10) and (11) in Fig. 3 are sensitive to both and . In the present work, however, we used the lines (1) and (7) as indicators and (8), (10)3, and (11) as indicators, respectively. By doing so, we estimated ([O iii]), ([O ii]), ([N ii]), and ([S ii]) by adopting ([O iii]), [O ii], Ne([N ii]), and ([S ii]). Since we could not de-blend [N i] 5198/5200 Å (its ratio is a density indicator for the neutral region), we used the far-IR [O i] ratio.

Liu et al. (2004b) reported five and four values based on the CELs seen in the ISIS spectra augmented by lines detected in the ISO spectra (see their Table 7). Taking advantage of the fine-structure lines detected at higher sensitivity and better spatial resolution in the Spitzer and Herschel spectra, we calculate seven and eight values. Our values of the CEL and are generally consistent with those determined by Liu et al. (2004b).

The - diagnostic diagram (Fig. 3) suggest that the bulk of the ionized gas appears to have between 6 000 K and 12 000 K. Thus, we adopt a constant  K to derive values. The derived ([Ne iii]) value is more than one order of magnitude larger than the other values. To double-check the above, we analyze the Spitzer/IRS spectra of the nebula nearby the central star obtained with the higher-dispersion SH and LH modules (of the and slit dimensions, respectively; not shown in Fig. 2). From the SH and LH spectra alone, we derive ([Ne iii]) =  cm and ([S iii]) =  cm. Because the spatial coverage of the SH and LH modules is very restrictive around the central star, the higher and values may be influenced heavily by the conditions in the vicinity of the central star. Previously, the [O iii] 52/88 µm ratio in the central part of the cavity yielded 350 cm.

Next, we calculate based on the derived values. The average of = 260 cm among ([S ii], [O ii], [N ii]) is adopted to calculate ([S ii]) and ([N ii]) (ID: 10). To compute ([Ar iii] and [Ne iii]), ([O iii]) of 220 cm is adopted. To calculate ([O iii]), ([O ii]), and ([N ii]) (ID: 9) accurately, we subtract contributions from O, O, and N RLs to the [O iii] 4363 Å, [O ii] 7320/30 Å, and [N ii] 5755 Å lines, respectively, i.e., ([O iii] 4363 Å), ([O ii] 7320/30 Å), and ([N ii] 5755 Å).

We calculate ([O iii] 4363 Å) with


(eqn 3, Liu et al., 2000) for which the O/H ratio (3.02(–5), see § 3.1.2) is computed using the ([O iv] 25.9 µm)/(H) ratio assuming ([Ne iii]) and ([O iii]). In the end, ([O iii] 4363 Å) turns out to be 0.73  of the observed ([O iii] 4363 Å). After we subtract ([O iii] 4363 Å) from the observed ([O iii] 4363 Å), we obtain ([O iii]) by adopting ([O iii]).

([O ii] 7320/30 Å) is calculated with


(eqn 2, Liu et al., 2000) where we adopt the O/H ratio (, see § 3.1.2) derived from the ([O iii] 88.3 µm)/(H) ratio assuming ([O iii]) and ([O iii]). ([O ii] 7320/30 Å) turns out to be 2.19  of the observed ([O ii] 7320/30 Å). After we subtract the recombination contribution from the observed ([O ii] 7320/30 Å), we obtain ([O ii]) by adopting = 260 cm.

Finally, we estimate ([N ii] 5755 Å) using


(eqn 1, Liu et al., 2000) where the N/H ratio (7.01(–5), see § 3.1.2) was calculated using the ([N iii] 57 µm)/(H) ratio assuming ([O iii]) and ([O iii]). ([N ii] 5755 Å) is 0.54  of the observed ([N ii] 5755 Å). After we subtract ([N ii] 5755 Å), we obtain ([N ii]) (ID: 9) by adopting = 260 cm. In § 4 below, we verify the above estimates of the RL contributions based on the best-fit modeling results.

We also determine (He i), which is necessary to estimate He and He abundances, using the He i (7281 Å)/(6678 Å) ratio with the He i recombination coefficients in the case of = 100 cm given by Benjamin et al. (1999). The and pairs derived and adopted from the present plasma diagnostics are summarized in Appendix Table C1.

Ionic abundance derivations

We calculate CEL ionic abundances by solving the equation of population at multiple energy levels with the adopted and (Appendix Table C1, which also lists the adopted and to calculate RL He and C abundances): the resulting ionic abundances are listed in Appendix Table C2. We give the 1- uncertainty for each ionic abundance estimate, which is propagated from 1- uncertainties of line fluxes, (H), , and . Ionic abundances are derived for each of the detected line intensities when more than one lines for a particular target ion is detected. In such cases, we adopt the weighted-average of all of the derived abundances listed at the last line for that particular ion in Italics.

The resulting ionic abundances based on different lines in the optical nebular, auroral, and trans-auroral transitions to IR fine-structure lines turn out to be generally consistent with each other within the associated uncertainties for most of the cases. This indicates that our choice of the - pair for each ionic species is robust and that the adopted scaling of the mid- and far-IR line fluxes to the optical H line flux via the adopted photometry data is reasonable. However, there are a few exceptions, which we briefly discuss below.

As pointed out above, the spatial coverage of the nebula in spectroscopic observations is not complete and uniform: especially, the ISIS spatial coverage in the optical missed the brightest E and W “rim” regions, in which low-excitation and neutral lines are particularly strong (Fig. 2). This explains why the O abundances derived from optical lines are much smaller than the abundance based on the [O i] 145 µm line (by a factor of ). Hence, if we were to assume O/H = based solely on the [O i] 145 µm line, we would have N/H = and S/H = by adopting a factor of . Nevertheless, for the O/H abundance we adopt the average of the observed two optical (6300 and 6364 Å) and single far-IR (145 ) lines, because there is no way to ascertain how much line flux is missed by incomplete spatial coverage of the nebula.

To determine the He abundance, we do not include the He i 4712 Å line because the blue wing of this line seems to be contaminated by the [Ar iv] 4711 Å line. Assuming that He is indeed (), (He i 4712 Å) and ([Ar iv] 4711 Å) have to be and , respectively4. The Ar abundance derived from this expected ([Ar iv] 4711 Å) is , which is consistent with the Ar abundance derived from ([Ar iv] 4740 Å).

To derive the RL C abundance, we use the C ii 4267 Å line with its effective recombination coefficient in Case B for = 10 cm defined as a polynomial function of by Davey et al. (2000). This is justified because while the effective recombination coefficient is not available for the case of = 100 cm that is more appropriate here, the RL abundances are in general insensitive to for  10 cm. As for , we adopt ([Ar iii]) because the ionization potential (I.P.) of C is similar to that of Ar.

Overall, we conclude that our derived ionic abundances are improved with new CEL detections in the mid- and far-IR (such as Ne, S, Si, Cl, and Ar) made with Spitzer and Herschel observations, more robust adaptation of and for targeting ions, and the use of a larger number of lines in various ionization stages, compared with those calculated previously by Liu et al. (2004a).

Elemental abundance

By introducing the ionization correction factor (ICF; see, e.g., Delgado-Inglada et al. 2014 for more detail), we infer the nebular abundances of the observed nine elements in the ionized part of the nebula based on their observed ionic abundances. In Appendix Table C2, the ICF(X) value of the element “X” and the resulting elemental abundance, X/H = ICF(X) X/H, are listed in bold at the last line for each element. Here, we exclude C, N, and O from abundance calculations for the respective elements, as these ions are considered to be present mostly in the PDR surrounding the ionized part of the nebula. In Table 3, we compare the derived elemental abundances (X) corresponding to , where (in column 2) and the relative Solar abundances (X/H; in column 3).

We perform an ionization correction using the ICF based on the I.P. of the element in question, except for He, O, Ne, and S (i.e., ICF for these four elements is taken to be unity because unobserved high excitation lines are considered negligible). We will compare these ICFs based on the I.P. and the predicted ICFs by the best-fit modeling in § 4.

In performing ionization correction, the ICF for N, Si, Cl, and Ar is set as follows. We assume that the N abundance is the sum of N, and adopt ICF(N) ICF(O), which is equal to the ) ratio. Similarly, we assume that the Si abundance is the sum of Si, and adopt ICF(Si) ICF(S), which corresponds to the S/S ratio. For Cl and Ar, we assume that the Cl and Ar abundances are the sum of Cl and Ar, respectively. Then, we adopt ICF(Cl) ICF(Ar) ICF(S), which corresponds to the ) ratio.

As for the ICF(C), we originally adopt ICF(C) ICF(N) corresponding to the N/N ratio. With this ICF(C), the derived RL C abundance using the RL C ii 4267 Å line would come out to be . Note that we do not include the CEL C abundance for the elemental C abundance because (1) the [C ii] 157 µm line arises mostly from the PDR as stated above and (2) the nature of these lines is different (C is of RL while C is of CEL).

However, this RL C abundance would be extremely unlikely for NGC 6781. The average abundance between [Cl/H] and [Ar/H] derived for NGC 6781 suggests that the metallicity () of the object is close to the solar metallicity (see also § 3.1.5). Then, such a high RL C abundance is very difficult to explain by current AGB nucleosynthesis models (e.g., Karakas, 2010) for stars with the solar metallicity (, ). Hence, the derived RL C abundance of appears to be overestimated.

It is known that C, N, O, and Ne ionic abundances derived from RLs are sometimes found to be larger than the corresponding abundances obtained from CELs in PNe and H ii regions. This issue is known as the abundance discrepancy problem. (see, e.g., Liu (2006), for more detail). In spite of a number of attempts to explain such ionic/elemental abundance discrepancies, no consensus has been reached yet. Thus, we need other options to estimate the C abundance in light of the abundance discrepancy problem. One option is to compute the expected CEL C abundance by scaling the measured RL C abundance with the average C(RL)/C(CEL) ratio because no UV spectrum is available for NGC 6781.

Previously, Delgado-Inglada & Rodríguez (2014) showed general agreement between measured and scaled CEL abundances, the latter of which was scaled from measured RL abundances with the average C(RL)/C(CEL) ratio of among 37 Galactic PNe (their Table 5). While it is yet unknown whether there is a correlation between the RL and CEL C abundances, the relatively small standard deviation of the measured ratios would indicate that this option has some merit. Because there are no other alternatives, we adopt this option for the present study and use the average C(RL)/C(CEL) ratio of found among 58 PNe in Milky Way and Magellanic Clouds (Otsuka et al., 2011) to obtain the scaled expected CEL C of .

This expected CEL C of () would be more reasonable than the measured RL C abundance of with respect to current AGB nucleosynthesis models for the solar abundance stars (e.g., Karakas, 2010). In addition, Delgado-Inglada & Rodríguez (2014) reported a C(RL)/C(CEL) ratio of 3.63 for NGC 6720, which possesses the central star and nebula properties very similar to those of NGC 6781 (see § 3.4.1).

Further on the C and Cl abundances

Because our present analysis and the previous analysis done by Liu et al. (2004a, listed in Table 3, column 4) are based on the same ISIS optical spectrum, both results should be consistent with each other. However, this is not the case for C and Cl.

The discrepancy in (Cl) arises because we adopt the Cl abundances of and and the corresponding ICF(Cl) value of 1.17, while Liu et al. (2004a) used the Cl abundance of only with the corresponding ICF(Cl) of 3.394. In addition, the adopted could contribute to the discrepancy because the Cl ionic abundances are determined using their CEL lines, whose emissivities are sensitive to . Overall, we would argue again that our (Cl) value is more improved than the previous estimate because we have more robust for the ionic Cl abundances and we derive a Cl abundance that would reduce uncertainties in ICF(Cl).

X (X) [X/H] (X) (X) (X)
(1) (2) (3) (4) (5) (6)
He 11.06 0.17 +0.13 0.17 11.08 11.05 11.06
C(RL) 9.61 0.29 +1.22 0.30   9.17   8.52   9.06
C(CEL) 8.56 – 9.00 +0.17 – 0.61            
N 8.15 0.09 +0.29 0.15   8.38   8.39   8.42
O 8.76 0.04 +0.03 0.08   8.65   8.94   8.94
Ne 8.15 0.05 +0.10 0.11   8.22   8.12   8.27
Si 7.03 0.27 –0.50 0.28       7.57   7.59
S 6.91 0.06 –0.25 0.06   6.97   7.42   7.44
Cl 5.16 0.42 –0.09 0.42   5.43
Ar 6.49 0.10 –0.01 0.14   6.35

Note. – The number density ratio relative to hydrogen is , where . The CEL C abundance, C(CEL), is an expected value.

Table 3: Elemental abundances of NGC 6781 derived in the present analysis, compared with the solar abundances (column 3; , where is taken from Lodders 2010), previous empirical analysis (column 4; by Liu et al. 2004a), and model predictions (columns 5 and 6; by Karakas 2010; see § 3.1.5).

The discrepancy in RL (C) is due to different values of (C ii 4267 Å) (might be caused by different adopted (H)) and adopted ICF(C): our (C) and ICF(C) values are and 2.03 whereas theirs are and 1.624, respectively. In general, C is a very important element as a coolant of the ionized gas component and also a source of C-based molecules in PNe. Thus, we would discuss the C abundance further in this section.

Our expected C(CEL) of ((C)) adopted in the previous section, in comparison with the observed O(CEL) of ((O)), would suggest a slightly C-rich nature for NGC 6781 (C/O number density ratio of ). Indeed, the Spitzer/IRS mid-IR spectrum (Fig. 1, inset) shows polycyclic aromatic hydrocarbon (PAH) emission at  µm (mostly from ionized PAH) and at 11.3 µm (from neutral PAH) and dust continuum due to amorphous carbon, while the spectrum does not clearly show any O-rich dust features such as amorphous silicates at 9 µm and 18 µm and crystalline silicates around 30 µm.

Guzman-Ramirez et al. (2014) reported detection of PAH emission in O-rich PNe in the Galactic bulge and suggested that PAHs could be formed in the compact/dense torus (i.e., the “waist” region of bipolar PNe) using C atoms liberated from CO molecules by photodissociation. At this point, there is no clear evidence to suggest this possibility for NGC 6781 based on the spatially-resolved spectroscopic data.

If we adopt RL C of and ICF(C) of 1.634 as previously used by Liu et al. (2004a) and convert the RL C abundance to the CEL C abundance by the average C(RL)/C(CEL) ratio of 4.10 (Otsuka et al., 2011), we would obtain the expected CEL C abundance of 3.61(–4), which would correspond to (C) of 8.56. This would result in a C/O ratio of 0.76, indicating that NGC6781 is slightly O-rich. Hence, the possibility of NGC 6781 being O-rich is not completely ruled out.

As seen above, the C abundance depends on many factors, from the (C ii 4267 Å) measurements to the ICF(C) and C(RL)/C(CEL) values adopted. Therefore, in the present work, we opt to allow a range of the expected CEL abundance for NGC 6781 as (correspondingly, ) based on the arguments presented above.

Comparison with the previous model predictions

We compare the derived with the values predicted by AGB nucleosynthesis models. As for the metallicity of the progenitor of NGC 6781, it is best to reference elements that can never be synthesized within AGB stars. Thus, we adopt Cl and Ar as good indicators. The average between the observed [Cl/H] and [Ar/H] values of corresponds to .

However, the S abundance () suggests a much lower . So far, this S abundance anomaly has been found in many Milky Way and M31 PNe (Henry et al., 2012, see their Fig. 1). Henry et al. (2012) concluded that the sulfur deficit in PNe is generally reduced by increasing the S abundance and selecting a proper ICF(S). Such an S depletion may indicate that a significant part of the atomic S mass is locked up as sulfide grains in the nebula (e.g., MgS and FeS in C- and O-rich environments, respectively). However, the Spitzer/IRS spectrum displays neither the broad 30 µm feature often attributed to MgS nor narrower emission features around 30 µm attributed to FeS. The discrepancy between the observed and the AGB model S abundances may thus be related to the adopted reaction rates; Shingles & Karakas (2013) demonstrated a possibility that the S depletion could be explained by introducing a large Ne(,)Mg reaction rate. Here, we propose that the apparently low [S/H] abundance is attributed to missing fluxes of low-excitation [S ii] lines as discussed above (by adopting the revised S/H in § 3.1.2, we would obtain , which is consistent with (S)).

Now, we compare our empirically-derived elemental abundances with those predicted with AGB nucleosynthesis models o f stars (Karakas, 2010) in Table 3: the values in columns (5) and (6) are the predicted values for initially 2.25  and 3.0  stars, respectively. To assess the goodness of fit of the model prediction, we evaluate chi-square values () between our derived abundances and the model-predicted abundances for stars in the initial mass range from 1.5 to 4.0 . Adopting the lower CEL abundance limit of , a good fit to the observed is achieved with the 2.25  model (reduced = 15.5).

Meanwhile, adopting the upper CEL abundance limit of , the values suggest that the observed is most consistent with the 2.5  model (reduced = 16.15). The reduced value = 17.5 of the 3.0  model is equally good. Therefore, based on these results we conclude that the initial mass of the CSPN is between 2.25 and 3.0 .

3.2 The molecular gas component

Given the number of molecular lines seen in the spectra, especially with the rare OH detection (Aleman et al., 2014), NGC 6781 has to be treated as a PN rich in neutral gas. Then, it is critical to include the PDR of the nebula for a complete understanding of all of its components (ions, atoms, molecules, and dust). In this section, therefore, we investigate the physical conditions of the most abundant species in the PDR, H, to articulate our understanding of the PDR in NGC 6781.

Physical conditions: spatial distribution

Figure 4: (a) A narrow band image of H S(1) at 2.122 µm taken with the 3.6-m CFHT/WIRCAM, overlaid with yellow contours of [N ii] 6583 Å emission taken with the 2.5-m NOT/ALFOSC (Phillips et al. 2011; provided to us by M. A. Guerrero). (b) A close-up of the central part of the nebula, showing the filamentary structure of the H emission in the central region, from which the adopted spectra arose. The location of the central star is also indicated.

We obtain the H image taken with the Wide-field Infrared Camera (WIRCAM, Puget et al., 2004) on the 3.6-m Canada France Hawaii Telescope (CFHT) from the Canadian Astronomy Data Centre (CADC). The observations were done on 2006 April 14 (PI: S. Kwok, Prop. ID: 06AT03) through Taiwan CFHT time. The basic calibrated data set retrieved from the CADC archive is reduced into a single image after bad pixel masking and geometric distortion correction using IRAF. Fig. 4 shows the H S(1) image at 2.122 µm overlaid with contours of [N ii] 6583 Å emission and the close-up of the central region from which emission of the spectra adopted in the present study arose (cf. Fig. 2). Fig. 4a shows that the spatial distribution of the molecular gas component in NGC 6781 seen via H emission is very similar to that of the cool low I.P. gas component seen via [N ii] emission (and also via H emission; Fig. 2). The same similarities in the spatial distributions are seen between the dust and ionized gas components delineating the nearly pole-on cylindrical barrel structure (Fig. 3 of HerPlaNS1). Highly localized distributions of the molecular gas component are apparent from the filamentary appearance of the H emission (Fig. 4b). These H filaments (and maybe clumps, too) are patches of H survived in the ionized region.

Physical conditions: shocks vs. UV radiation

Table 4 summarizes near- and mid-IR H lines detected in NGC 6781. As reported by Phillips et al. (2011) and Mata et al. (2016), pure rotational H lines are detected in the Spitzer/IRS spectra (Fig. 1, inset). Observations made by Arias & Rosado (2002) show that the intensity of H S(1) at 2.248 µm is much fainter than that of H S(1) at 2.122 µm, which indicates collisional excitation. The kinematic studies of Hiriart (2005) pointed to a post-shock origin for the H emission. If the observed H lines are radiatively excited through the absorption of far-UV photons (11 – 13 eV) in PDRs, the upper vibrational level would have to have a larger population, resulting in a relatively high H (2.248 µm)/(2.122 µm) via UV fluorescence (e.g., Kwok, 2007). Collisional excitation, on the other hand, can occur in both shocks and PDRs. Excitation mechanisms of H in PNe are examined by evaluating H (2.248 µm)/(2.122 µm) ratio (e.g., Otsuka et al., 2013), even though it is not easy to do with -band data alone.

Interestingly, the expansion velocity of H ( km s, Arias & Rosado, 2002; Hiriart, 2005) is found to be greater than the expansion velocity measured from the [O iii] line (10 km s, Weinberger, 1989) and [N ii] line (12 km s, Arias & Rosado, 2002). Hiriart (2005) concluded that the average H S(1) surface brightness could be explained by shocks at  km s heading into the pre-shock region of the H density at  cm.

We investigate the conditions in the H emitting regions by comparing the flux ratios of mid-IR H lines to the S(3) line at 9.67 µm with the theoretical continuous shock (C-shock) models by Flower & Pineau Des Forêts (2010). The observed (17.04 µm)/(9.67 µm) ratio suggests a match for a model with the shock velocity of  km s and pre-shock hydrogen density of  cm, while the observed (12.29, 8.02, 6.91, 6.11, 5.51 µm) to (9.67 µm) ratios point to a model with  km s and  cm. Here, the possible line flux contamination from the H i 12.3 µm line to the H 12.29 µm line, estimated to be (H i 12.3 µm) = 0.971 when (H)  = 100 in the case of = 10 K and = 200 cm, is removed.

Bachiller et al. (1993) reported a CO expansion velocity of 22 km s. Recently, Bergstedt (2015) reported a velocity of 16 km s via 3-D structure modeling using CO velocity maps. A model by Flower & Pineau Des Forêts (2010) with a shock velocity of  km s and pre-shock hydrogen density of  cm would explain the observed far-IR CO line flux ratios with respect to the CO line at 371.6 µm obtained from our Herschel PACS and SPIRE spectra (HerPlaNS1).

Based on the arguments above, excitation of H and CO lines in NGC 6781 appears to be caused by thermal shocks at a velocity in the range of  km s impinging onto the pre-shock region at  cm. These shocks may be be the consequence of interactions between the slow AGB wind and fast PN wind emanating from the CSPN in the context of the PN evolution. The slow-fast wind interactions could cause diffuse X-ray emission in the interaction regions. No X-ray detection in NGC 6781 may thus be because of extinction (see § 4.2.8). Together with the filamentary/clumpy appearance of the H emission regions (Fig. 4), we would conclude that these structures represent high-density regions delineating the locations of thermal collisional excitation embedded in an lower density ionized gas. Such high H clumps (so called “cometary knots”, O’Dell & Handron, 1996) within the ionized gas are detected in nearby PNe (see e.g., O’Dell et al., 2002). Recently, Manchado et al. (2015) detected cometary H knots within the ionized gas region in the bipolar PN NGC 2346.

One might think that the H distribution in NGC 6781 is similar to that in NGC 7293 (Helix nebula), in which the H emission is considered to arise from H clumps. For NGC 7293, there is no evidence to suggest that the H emission from its cometary knots is due to shocks (Aleman et al., 2011, reference therein). Another possible H excitation mechanism is due to the structure and steady state dynamics of advective ionization front/dissociation front (Henney et al., 2007). However, our Cloudy models with turbulence velocity of  km s in the nebula by following Henney et al. (2007) failed to reproduce the observed H line fluxes. While these are definitely issues that needs to be resolved in future, we tentatively conclude that the observed H emission in NGC 6781 has a shock origin based on the arguments presented above.

Physical conditions: H excitation diagram

Figure 5: The excitation diagram of pure rotational transitions of H lines. We fit the observed data (Table 4) with a single excitation temperature (top frame; a; with all but the 17.04 µm line;  K) and with two excitation temperatures (bottom frame; b; with all lines;  K and  K).
\arraybackslashTransition Average intensity
(µm) \arraybackslash (erg s cm sr)
\arraybackslash0-0 S(1)
\arraybackslash0-0 S(2)
\arraybackslash0-0 S(3)
\arraybackslash0-0 S(4)
\arraybackslash0-0 S(5)
\arraybackslash0-0 S(6)
\arraybackslash0-0 S(7)
\arraybackslash1-0 S(1) 2.70(–4)
Table 4: The average H intensities of NGC 6781 measured with Spitzer/IRS. See, also, Fig. 1.

Assuming that H lines are thermally excited and are in local thermodynamic equilibrium (LTE), the H excitation temperature and column density can be estimated via an excitation diagram. The H column density in the upper state is written as


where (H) is the H line intensity in erg s cm sr, is the transition probability taken from Turner et al. (1977), is the Planck constant, and is the speed of light. In LTE, the Boltzmann equation relates to the excitation temperature (H) via


where is the vibrational degeneracy, is the energy of the excited level taken from Dabrowski (1984), is the Boltzmann constant, and is the rotational constant (60.81 cm).

In Fig. 5, we plot the vs.  for each of the H lines detected in NGC 6781 (Table 4). The of the H S(1) (magenta circle) is calculated using the average line intensity of  erg cm s sr (Hiriart, 2005). The rotational diagram suggests that the bulk of the H 17.04 µm line emission is produced in a region with different physical conditions from the other H line emitting regions.

First, we determine the conditions of the H emitting region by fitting the line fluxes at 12.29, 9.67, 8.02, 6.91, 6.11, 5.51, and 2.12 µm (i.e., all but 17.04 µm) with Equation 6 using a single excitation temperature (Fig. 5a):  K and  cm. The derived ) is comparable to  K and  K, previously derived by Phillips et al. (2011) and Mata et al. (2016), respectively (with a single temperature model using all but the 2.12 µm and 17.04 µm lines).

Next, we fit all H lines (including 17.04 µm) using two excitation temperatures (Fig. 5b). The warm component is found to have  K and  cm, whereas the cold component is found to have  K and  cm (while lack of the H 0-0 S(0) line at 28.2 µm makes the fitting results relatively less certain). Nonetheless, the 17.04 µm line is expected to arise from such colder and denser regions.

Empirically determined molecular gas mass

To conclude this subsection, we estimate the mass of the molecular gas component in the nebula by adopting the distance of 0.46 kpc (§ 3.4.1). Based on the H and CO emission maps (Hiriart 2005 and Bachiller et al. 1993, respectively), we see that molecular emission increases at 54-55″ away from the CSPN with the thickness of 12″. Using H densities of the warm and cold components as derived above ( cm and  cm, respectively), we estimate the H gas mass of and for the warm and cold components, respectively.

Previously, we derived  cm (excitation temperature at 60 K) based solely on our Herschel spectra (HerPlaNS1). Bachiller et al. (1997) measured  cm (excitation temperature at 25 K) based on sub-millimeter data. Assuming that each of the above estimates based on data in the different wavelength/temperature realms would represent the warm and cold component, respectively, the warm and cold CO gas masses are estimated to be and , respectively. These estimates are combined to yield the total molecular gas mass (of H and CO) of .

The empirical (CO)/(H) ratio turns out to be and for the warm and cold temperature regions, respectively. Assuming that the (CO)/(H) ratio translates roughly to , we can estimate (C) of for the molecular gas component. Compared with the adopted CEL expected (C) of for the ionized gas component, 11-60 % of the C-atoms were estimated to be locked in as molecules.

3.3 The dust component: summary of HerPlaNS I

The surface brightness distribution of thermal dust continuum emission from NGC 6781 is spatially resolved in far-IR Herschel broadband images (?)see Fig. 3 of][]Ueta:2014aa. The bright ring structure with 60″ outer radii represents the bulk of the nearly pole-on cylindrical barrel structure (originally proposed by Schwarz & Monteiro, 2006), and the elongated nebula of 200″ in the total north-south extent indicates the distribution of dust along the polar axis of the nebula. The spatial extent of thermal dust continuum emission in far-IR wavelengths is nearly identical with that of atomic gas and molecular emission lines in optical and near-IR wavelengths.

Previously, we performed spectral energy distribution (SED) fitting of the Herschel 70/160/250/350/500 µm images using a modified blackbody function, and found that dust grains are composed mostly of amorphous-carbon based material (i.e., the power-law dust emissivity index is 1 across the nebula) having the dust temperature in the range between 26 and 40 K (HerPlaNS1). Moreover, after removing the contribution to the continuum flux in the far-IR by fine-structure lines and molecular emission lines (amounts to 8-20 % of the total flux), spectral fitting of the integrated far-IR fluxes yielded =  K and = . Indeed, the Spitzer/IRS spectrum (Fig. 1, inset) shows PAH bands and featureless dust continuum, This is consistent with the dusty nebula of NGC 6781 containing more amorphous carbon dust and PAHs than amorphous silicate dust.

3.4 The central star

Distance, luminosity, and effective temperature

PNe (He) (C) (C) (N) (O) (O) (Ne) (S) (Cl) (Ar) (kK) (cm s) References
NGC 6781 11.06 9.61 8.56-9.00 8.15 8.76 8.15 6.91 5.16 6.49 80 – 123 6.0 – 7.0 (1), (2), (3), (4)
NGC 6720 11.05 9.10 8.59 8.22 9.18 8.80 8.23 6.86 5.19 6.54 80 – 135 6.9 – 7.0 (5), (6), (7)

References. – (1) This work for abundances (see § 3.1.3); (2) Schwarz & Monteiro (2006) for and via photoionization modeling; (3) Rauch et al. (2004) for and via stellar absorption fitting; (4) Liu et al. (2004a) for abundances; (5) McCarthy et al. (1997) for and via stellar absorption fitting; (6) Napiwotzki (1999) for and via stellar absorption fitting; (7) van Hoof et al. (2010) for via Cloudy photoionization modeling.

Table 5: Similarities between NGC 6781 and NGC 6720.

A vast variety of distance estimates are proposed for NGC 6781, including 0.3 kpc (Tajitsu & Tamura, 1998; Phillips, 2002), 0.7 kpc (Stanghellini et al., 2008; Frew et al., 2016), 0.9 kpc (Maciel, 1984), 0.95 kpc (Schwarz & Monteiro, 2006), and 1.27 kpc (Ali et al., 2013), to name a few. For the present study, rather than adopting any of the previous investigations, we elect to determine our own value by comparing the observed photometry of the CSPN (Fig. 1, Table A1) with the post-AGB evolutionary tracks produced by Vassiliadis & Wood (1994) augmented with a grid of synthesised spectra by Rauch (2003). Although several new evolutionary tracks have been produced since then, there has been no AGB nucleosynthesis models constructed based on such new tracks. In comparing observed data with theoretical models, we would regard internal consistencies between models more important. Especially when we aim at determining the state of evolution of the CSPN of NGC 6781, the most critical is adopting AGB nucleosynthesis models that are consistent with evolutionary tracks. Therefore, in the following discussion, we adopt the AGB nucleosynthesis models by Karakas (2010) based on Vassiliadis & Wood (1994).

We start by estimating the CSPN luminosity using a grid of non-LTE line-blanketed plane-parallel hydrostatic atmospheric models generated by Rauch (2003) as templates. We adopt the solar abundance () models for the CSPN based on the results of our own nebular abundance analysis presented in § 3.1.5.

To characterize the stellar atmosphere fully, we also need the effective temperature and surface gravity of the CSPN. Previously, Rauch et al. (2004) suggested  K and  cm s based on the stellar absorption line fitting. If this were true, the CSPN would have been still burning hydrogen in a thin surface layer while increasing . However, detection of strong He ii 4686 Å and [O iv] 25.88 µm lines in the ISIS and Spitzer/IRS spectra, respectively (Fig. 1 and Table B1) requires  K, refuting the previous suggestion. The noisy spectrum due to the faintness of the CSPN might have compromised the previous absorption line fitting analysis.

Thus, we decide to look for the appropriate and values in a PN similar to NGC 6781 in terms of nebula and CSPN properties. Amongst Galactic PNe, NGC 6720 is very similar to NGC 6781 in many respects, especially in their abundance pattern as shown in Table 5. Spectroscopically, both PNe show PAH features and pure rotational H lines in their Spitzer/IRS spectra (Phillips et al., 2011; Cox et al., 2016) as well as rotational-vibrational H emission (e.g., Hiriart, 2005; van Hoof et al., 2010). Both PNe possess a structure due to a heavy equatorial concentration (i.e., a generic bipolar/barrel shape) viewed nearly pole-on (Schwarz & Monteiro, 2006; Sahai et al., 2012; Ueta et al., 2014).

The CSPN of NGC 6720 has a  kK based on the absorption line analysis done by McCarthy et al. (1997) and Napiwotzki (1999). Thus, based on the similarities listed above we adopt  kK and  cm s for the CSPN of NGC 6781 as well. Consistent results were previously obtained from detailed SED fitting with Cloudy photoionization models of NGC 6720 (van Hoof et al., 2010, also see § 4).

Then, we scale the synthesized Rauch model spectra of the adopted CSPN characteristics of  kK with a constant 10 kK step with  cm s fixed so that the observed photometry from the WFC -band to WFCAM -band (see Table A1) matches with the model spectra (Fig. 6, showing the  kK case). The scaled spectra are integrated to yield , which is then parameterized with and the distance in the form of


where is in kpc and is in K. Note that is not very sensitive to . For instance, increases only by 0.8 % when is reduced from the adopted 6.9 cm s to 6.6 cm s. Thus, our choice of single value is warranted.

Figure 6: The synthesised spectrum of a star with = 0.02, = 120 kK, and = 6.9 cm s (Rauch, 2003, the black line) fit with the observed photometry points of the CSPN (the blue circles; Table A1).
Figure 7: The distance-fitting comparison among the post-AGB evolutionary model tracks of 1.5, 2.25, 2.5, and 3.0  initial-mass stars (Vassiliadis & Wood, 1994) and the CSPN luminosity, , computed for = 0.95 and 0.46 kpc (blue squares and red circles, respectively) and , and 140 kK (from right to left, respectively), when  cm s. Also shown is the pair adopted by Schwarz & Monteiro (2006), with which they concluded = 0.95 kpc (black triangle). The light-blue box indicates the parameter range based on our Cloudy model calculations (§ 4) at  kpc. See, also, Fig. 12.

Finally, we compute at =  kK and for a range of , and plot the resulting (, ) pairs over the post-AGB evolutionary tracks of the 1.5, 2.25, 2.5, and 3.0  initial mass stars produced by Vassiliadis & Wood (1994), as shown in Fig. 7. Our choice of the initial mass of the adopted post-AGB evolutionary tracks is dictated by the results of our abundance analysis that indicated the CSPN initial mass being between 2.25 and 3.0  (§ 3.1.5). Also, the previous analysis by Schwarz & Monteiro (2006) suggested the CSPN initial mass of 1.5 .

We find that  kpc fits the initially post-AGB evolutionary tracks the best for the adopted range (light-blue box in Fig. 7). Therefore, we adopt  kpc, which is the the intermediate value between 0.34 and 0.52 kpc (red circles in Fig.7). Accordingly, we find . This evolutionary track fitting suggests that the CSPN of NGC 6781 is in the cooling phase. The results of the fitting are not significantly altered even when we adopt more recent post-AGB evolution tracks such as the ones computed by Miller Bertolami (2016) ( kpc, using the post-AGB evolutionary tracks for 2.0  and 3.0  stars with ; see also Fig. 12).

Previously, Schwarz & Monteiro (2006) concluded that the progenitor of NGC 6781 was a initial mass star based on their derived and values, provided  kpc suggested from their photoionization model fitting (black triangle in Fig. 7; also suggesting that NGC 6781 was in cooling phase). At  kpc, our estimates would be consistent with the 1.5  evolutionary track (blue squares in Fig. 7). However, the progenitor CSPN mass of NGC 6781 would most likely exceed 1.5  because of its empirically determined elemental abundances (§ 3.1.3) and H detection in this object (§ 3.2.2).

With a survey of H S(1) emission in Galactic PNe, Kastner et al. (1996) suggested that H-rich PNe evolved from relatively massive progenitors because H was exclusively detected in bipolar PNe (see also e.g., Guerrero et al., 2000). Bipolar PNe are known to be associated with massive () progenitors based on the distribution of bipolar PNe in the Milky Way with respect to that of elliptical PNe (Corradi & Schwarz, 1995). Hence, the detection of H supports our adaptation of the initial mass for the CSPN of NGC 6781 and the distance of 0.46 kpc based on the fitting.

The filamentary appearance of the nebula (Fig. 4) and low even in the central ionized regions (§ 3.1.1) are also suggestive that NGC 6781 is a highly evolved PN. Referring back to the similarity to NGC 6720, comparisons between and , where is based on Cloudy model fitting of the SED by van Hoof et al. (2010) with the evolutionary tracks by Vassiliadis & Wood (1994) for initially 3.0  stars of , also suggest that NGC 6720 is in the cooling phase.

If the CSPN of NGC 6781 were still in the final H-burning phase, the distance estimate would have to be  kpc. According to Vassiliadis & Wood (1994), is nearly constant at 6300  along the horizontal part of the post-AGB track for a 2.5  initial mass star with . In this case, the number of the ionizing photons is 4.25(+47) s for  K and  cm s. The Strömgren radius for this radiation field in a constant hydrogen density of 300 cm (see Table 2, Fig.8) with a filling factor () of unity would be 0.41 pc. This corresponds to the apparent radius of 237 at  kpc, which disagrees with the observed ionization radius of 55. Because the Strömgren radius is proportional to , it would be consistent with the observed ionization radius at  kpc if were 0.12. However, according to the empirical method introduced by Mallik & Peimbert (1988), the value of NGC 6781 is estimated to be 0.4 at  kpc and almost unity at 0.46 kpc. Therefore, we conclude that the CSPN of NGC 6781 already evolved off to the cooling track presently with and  kK at  kpc.

Possibility of the presence of a binary companion

At present, binary evolution would appear to be one of the most viable explanations for the formation of bipolar nebulae via the inevitable equatorial density enhancement (e.g., Jones & Boffin 2017). Our motivation to collect photometry measurements of the CSPN exhaustively in the UV to near-IR is also intended to establish the presence or absence of a near-IR excess, which would suggest the presence of a cooler binary companion.

From a comparison between the observed colors ( and ) and the grid of theoretical color indices as a function of , Douchin et al. (2015) argued that CSPN of NGC 6781 shows near-IR excess owing to an M1-M7 type companion star. However, we do not observe any IR excess in the SED of the CSPN (Figs. 1 and 6).

It is true that the IR excess detection can be influenced by the way the interstellar extinction is corrected for. With our adopted H, the extinction corrected and colors of the CSPN were and , respectively. If we used (corresponding to H) as adopted by Douchin et al. (2015), the respective and colors would become redder, and , which would be in perfect agreement with Douchin et al. (2015). This would negate the necessity for a M1-M7 type companion star.

Thus, whether NGC 6781 possesses a binary central system is still an open question because the evolutionary effects from the secondary, even if it existed, would still be negligible at this point, based on the observed spectra and photometry. Therefore, we would simply keep the adopted  kpc and other quantities for which there is distance dependency, in our analyses as outlined in the previous sections and in the subsequent modeling section.

4 Cloudy dusty photoionization models

4.1 Modeling approach

In the previous sections, we outlined how we mustered the most comprehensive observational data set yet assembled for NGC 6781 (§ 2) and performed various analyses to determine empirically the CSPN and nebula characteristics for this object (§ 3). In this section, we outline how we construct a realistic input numerical model of NGC 6781 for Cloudy (version C13.03, Ferland et al., 2013), comprising the CSPN and the nebula, the latter of which consists of the ionized/neutral/molecular gas and dust components, based on the collected data.

Our aim is to converge on self-consistent physical conditions of the entire NGC 6781 system from the highly-ionized region to the PDR through iterative model fitting that comprehensively reproduces all of the observational data that we collected: the spatially-integrated fluxes and flux densities from to UV to radio (37 broadband photometry fluxes, 19 flux densities, and 78 emission line fluxes) plus 8 elemental abundances. The empirically derived quantities of the CSPN and nebula provide the input parameters, while the observational data from the UV to radio provide the vital constraints in iterative fittings of the model parameters. For the sake of consistency, we substituted the same transition probabilities and effective collision strengths of CELs used in our plasma diagnostics and nebular abundance analyses in the Cloudy code.

4.2 The input model

SED of the CSPN

As the incident SED from the CSPN, we adopt the theoretical atmospheric model grid by Rauch (2003) for a star with and  cm s (see Fig. 6 for the case of ,  cm s, and  kK). We keep the distance of 0.46 kpc to NGC 6781, and vary and within the possible ranges, and  kK, as determined in § 3.4, during the iterative model fitting to search for the best-fit model parameters that would reproduce the observational data.

Nebular elemental abundances

For the elemental abundances of the nebula, we adopt the empirically-determined abundances (Table 3; § 3.1) as the input values. The nebular abundances are then refined via model iterations within 3- of the input values so that the best-fit abundances would reproduce the observed emission line intensities.

It should be pointed out here that the metal abundances would affect cooling of the nebula, and hence, would alter the nebula’s temperature and ionization structures. As we saw in § 3.1.3, the derivation of the C abundance is definitely a source of uncertainties. The only option of the empirical derivation available to us suggests the expected CEL C abundance of (Table 3). Hence, for the purpose of the present modeling, we set to be at the lower limit of 8.56 and keep it fixed during the model iteration. This will ensure that the best-fit model always satisfies at least the lower limit of the progenitor mass of 2.25  (see § 3.1.5).

The expected CEL of 8.56 is also consistent with the AGB nucleosynthesis model for the stars (Karakas, 2010). As we demonstrated in §3.4.1, NGC 6781 is very similar to NGC 6720 in terms of the elemental abundance pattern of the nebula and evolutionary state of the CSPN (Table 5). The adopted CEL of 8.56 for NGC 6781 is indeed very much consistent with that of 8.59 for NGC 6720. In addition, we adopt a constant C/C ratio of 20 determined by Bachiller et al. (1997).

As for the unobserved elements including heavy metals, we adopt the abundance values predicted with the AGB nucleosynthesis model of the 2.5  initial mass star with (Karakas, 2010). However, the Fe abundance is another exception, because we overpredict the Fe lines when setting as determined by Karakas (2010). The model ([Fe ii] 17.9 µm) and ([Fe iii] 4880 Å) line fluxes turn out to be 31.2 and 2.6 (with respect to H), respectively.

Nevertheless, such strong Fe lines are seen neither in the WHT/ISIS spectrum nor in the Spitzer/IRS spectrum. Therefore, we must adopt a lower Fe abundance. Previously, Liu et al. (2004a) measured in NGC 6720. Thus, we adopt , following the same similarity argument between NGC 6781 and NGC 6720 as in § 3.4.1. For other Fe-peak elements such as Cr, Mn, Co, and Ni, we adopt their solar values simply because their elemental abundances are unknown in NGC 6781.

Geometry of the nebula

Figure 8: The adopted geometry and hydrogen density ((H)) profile of NGC 6781 in Cloudy model.

Many authors suggested that NGC 6781 possessed a nearly pole-on cylindrical barrel structure, which surrounds the central cavity filled with tenuous highly ionized gas (e.g., Bachiller et al. 1993; Hiriart 2005; Schwarz & Monteiro 2006; Bergstedt 2015, as well as HerPlaNS1). Hence, with the 1-D code Cloudy, we represent the barrel wall structure by thin, concentric layers of ionized gas and dusty PDR. Such an “onion skin” configuration naturally explains the observed co-spatial distributions of various components at different temperature by the projection effect (Fig. 4a). While clumps/filaments of H surviving in the ionized region would be plausible (Fig. 4b), we simply adopt this “onion skin” configuration for the sake of 1-D model calculations, assuming that such molecular clumps/filaments would not significantly alter the nebular energetics.

However, we do take into account the barrel geometry of NGC 6781 by invoking the “cylinder” option of Cloudy, which approximates the cylindrical structure by removing polar caps from a hollow sphere (which is the default 1-D spherically symmetric configuration). We set the polar height of the cylinder to 90″, which is the average value between 72″ (suggested from the velocity channel maps in H; Hiriart, 2005) and 117″ (suggested from the velocity channel maps taken in CO at 345.796 GHz (866.96 µm); Bergstedt, 2015), assuming that the H and CO emission arose from the same regions because of the similarities between H and CO maps (Bachiller et al., 1993; Bergstedt, 2015). Fig. 8 shows a schematic of the adopted geometry.

Hydrogen density radial profile of the nebula

The input radial hydrogen density profile,