The Herschel Planetary Nebula Survey (HerPlaNS): A Comprehensive Dusty Photoionization Model of NGC6781^*

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 796 × 1'' slit 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). Figure 2 shows the central scanned region of 796 × 84'' 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).

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 36 × 57'' slits at P.A. of $-10^\circ$ in Figure 2) and LL (13.9–39.9 μm; the horizontal 168'' × 105 slit at P.A. of 86^◦ in Figure 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 (Figure 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 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 (5 × 5 in the $50^{\prime\prime} \times 50^{\prime\prime}$ apertures), while a single spectrum from the central bolometer of the SPIRE array is included (of $21^{\prime\prime}$ and $42^{\prime\prime}$ diameter in the short- and long-wavelength band, respectively; at both the "center" and "rim" positions as depicted with white boxes and gray circles, respectively, in Figure 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).

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:

$$\begin{eqnarray}&&I(\lambda )=F(\lambda )\cdot \,{10}^{c({\rm{H}}\beta )(1+f(\lambda ))},\end{eqnarray} \tag{ 1 }$$

where I(λ) is the dereddened line flux, F(λ) is the observed line flux, c(Hβ) is the reddening coefficient at Hβ, and f(λ) is the interstellar extinction function at λ computed by the reddening law of Cardelli et al. (1989) with ${R}_{V}=3.1$.

We measure the reddening correction factor c(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 ${T}_{{\rm{e}}}$ = 10,000 K and an electron density ${n}_{{\rm{e}}}$ = 200 cm⁻³ under the assumption that the nebula is optically thick to Lyα (so-called "Case B"; e.g., see Baker & Menzel 1938; also see Section 3.1.1 for the bases of these ${n}_{{\rm{e}}}$ and ${T}_{{\rm{e}}}$ values). The measured c(Hβ) turns out to be 0.951 ± 0.091 from the F(Hγ)/F(Hβ) ratio and 1.014 ± 0.033 from the F(Hα)/F(Hβ) ratio. Thus, we adopt c(Hβ) of 1.007 ± 0.031, which is a weighted mean of the above values. We do not correct for the interstellar extinction at longer wavelengths than K band because extinction would be negligible at those wavelengths. The final dereddening line fluxes measured in the adopted spectra are listed in Table 12. The quoted fluxes are normalized with respect to I(Hβ) = 100.

While we adopt these reprocessed 1D panchromatic spectra and duly measured dereddened line fluxes as representative of the bulk of the nebula, a word of caution appears appropriate at this point. As Figure 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 (Section 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 Section 3.1.2). This is primarily because the highest surface brightness regions (the E and W end of the central ring structure; Figure 2) are missed in the optical data and may be less strongly weighted than they should be in the far-IR data. We will return to these issues when we discuss model fitting in Section 4.

We determine the ${n}_{{\rm{e}}}$ and ${T}_{{\rm{e}}}$ pairs for the ionized/neutral gas component¹³ 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 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 ${n}_{{\rm{e}}}$ and ${T}_{{\rm{e}}}$ values are summarized in Table 2.

Table 2.  Summary of Plasma Diagnostics Using Nebular Lines

ID	Ion	${n}_{{\rm{e}}}$-diagnostics	Ratio	Result (cm−3)
1	[O i]	I(63 μm)/I(146 μm)	11.423 ± 2.039	590+1190
2	[S ii]	I(6716 Å)/I(6731 Å)	1.201 ± 0.048	230 ± 60
3	[O ii]	I(3726 Å)/I(3729 Å)	0.848 ± 0.035	270 ± 50
4	[N ii]	I(122 μm)/I(205 μm)	4.902 ± 0.991	280 ± 120
5	[S iii]	I(18.7 μm)/I(33.5 μm)	0.939 ± 0.092	1020 ± 300
6	[Ne iii]	I(15.6 μm)/I(36.0 μm)	13.789 ± 1.471	12 600 ± 7590
7	[O iii]	I(4959 Å)/I(88.3 μm)	1.438 ± 0.178	220 ± 50
ID	Ion	${T}_{{\rm{e}}}$-diagnostics	Ratio	Result (K)
8	[S ii]	I(6716/31 Å)/I(4069 Å)	14.891 ± 3.270	10 520 ± 1820
9	[N ii]	I(6548/83 Å)/I(5755 Å)	81.931 ± 2.956	10 800 ± 170
10	[N ii]	I(6548/83 Å)/	57.325 ± 6.201	12 360 ± 980
		I(122 μm+205 μm)
11	[O ii]	I(3726/29 Å)/I(7320/30 Å)	50.262 ± 1.949	9650 ± 200
12	[Ar iii]	I(7135 Å+7751 Å)/I(9.0 μm)	1.211 ± 0.098	9350 ± 400
13	[O iii]	I(4959 Å)/I(4363 Å)	52.943 ± 3.584	10 050 ± 210
14	[Ne iii]	I(3868 Å+3967 Å)/I(36.0 μm)	9.578 ± 0.806	10 340 ± 250
	He i	I(7281 Å)/I(6678 Å)	0.156 ± 0.021	7070 ± 1880

Download table as:  ASCIITypeset image

The ${n}_{{\rm{e}}}$–${T}_{{\rm{e}}}$ plot shown in Figure 3 summarizes how ${n}_{{\rm{e}}}$ and ${T}_{{\rm{e}}}$ 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 ${n}_{{\rm{e}}}$–${T}_{{\rm{e}}}$ curves derived from the ${T}_{{\rm{e}}}$-sensitive ratios, while the dashed lines are those from the ${n}_{{\rm{e}}}$-sensitive line ratios. Strictly speaking, the diagnostic lines labeled as (1), (7), (8), (10), and (11) in Figure 3 are sensitive to both ${n}_{{\rm{e}}}$ and ${T}_{{\rm{e}}}$. In the present work, however, we used the lines (1) and (7) as ${n}_{{\rm{e}}}$ indicators and (8), (10),¹⁴ and (11) as ${T}_{{\rm{e}}}$ indicators, respectively. By doing so, we estimated ${T}_{{\rm{e}}}$([O iii]), ${T}_{{\rm{e}}}$([O ii]), ${T}_{{\rm{e}}}$([N ii]), and ${T}_{{\rm{e}}}$([S ii]) by adopting ${n}_{{\rm{e}}}$([O iii]), ${n}_{{\rm{e}}}$[O ii], Ne([N ii]), and ${n}_{{\rm{e}}}$([S ii]). Since we could not deblend [N i] λ5198/λ5200 (its ratio is a density indicator for the neutral region), we used the far-IR [O i] ratio.

Zoom In Zoom Out Reset image size

Liu et al. (2004b) reported five ${n}_{{\rm{e}}}$ and four ${T}_{{\rm{e}}}$ 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 ${n}_{{\rm{e}}}$ and eight ${T}_{{\rm{e}}}$ values. Our values of the CEL ${n}_{{\rm{e}}}$ and ${T}_{{\rm{e}}}$ are generally consistent with those determined by Liu et al. (2004b).

The ${n}_{{\rm{e}}}$–${T}_{{\rm{e}}}$ diagnostic diagram (Figure 3) suggests that the bulk of the ionized gas appears to have ${T}_{{\rm{e}}}$ between ∼6000 and ∼12,000 K. Thus, we adopt a constant ${T}_{{\rm{e}}}$ = 10,000 K to derive ${n}_{{\rm{e}}}$ values. The derived ${n}_{{\rm{e}}}$([Ne iii]) value is more than one order of magnitude larger than the other ${n}_{{\rm{e}}}$ 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 47 × 113 and 111 × 223 slit dimensions, respectively; not shown in Figure 2). From the SH and LH spectra alone, we derive ${n}_{{\rm{e}}}$([Ne iii]) = 4930 ± 2780 cm⁻³ and ${n}_{{\rm{e}}}$([S iii]) = 1240 ± 60 cm⁻³. Because the spatial coverage of the SH and LH modules is very restrictive around the central star, the higher ${n}_{{\rm{e}}}$ and ${T}_{{\rm{e}}}$ 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 ${T}_{{\rm{e}}}$ based on the derived ${n}_{{\rm{e}}}$ values. The average of ${n}_{{\rm{e}}}$ = 260 cm⁻³ among ${n}_{{\rm{e}}}$([S ii], [O ii], [N ii]) is adopted to calculate ${T}_{{\rm{e}}}$([S ii]) and ${T}_{{\rm{e}}}$([N ii]) (ID: 10). To compute ${T}_{{\rm{e}}}$([Ar iii] and [Ne iii]), ${n}_{{\rm{e}}}$([O iii]) of 220 cm⁻³ is adopted. To calculate ${T}_{{\rm{e}}}$([O iii]), ${T}_{{\rm{e}}}$([O ii]), and ${T}_{{\rm{e}}}$([N ii]) (ID: 9) accurately, we subtract contributions from ${{\rm{O}}}^{3+}$, ${{\rm{O}}}^{2+}$, and ${{\rm{N}}}^{2+}$ RLs to the [O iii] λ4363, [O ii] λ7320/30, and [N ii] λ5755 lines, respectively, i.e., ${I}_{{\rm{R}}}$([O iii] λ4363), ${I}_{{\rm{R}}}$([O ii] λ7320/30), and ${I}_{{\rm{R}}}$([N ii] λ5755).

We calculate ${I}_{{\rm{R}}}$([O iii] λ4363) with

$$\begin{eqnarray}&&\displaystyle \frac{{I}_{{\rm{R}}}([{\rm{O}}\,{\rm{III}}]\,4363\,\mathring{\rm A} )}{I({\rm{H}}\beta )}=12.4{\left(\displaystyle \frac{{T}_{{\rm{e}}}}{{10}^{4}}\right)}^{0.59}\displaystyle \frac{{{\rm{O}}}^{3+}}{{{\rm{H}}}^{+}}\end{eqnarray} \tag{ 2 }$$

(Equation (3) in Liu et al. 2000), for which the O³⁺/H⁺ ratio (3.02(−5); see Section 3.1.2) is computed using the I([O iv] 25.9 μm)/I(Hβ) ratio assuming ${T}_{{\rm{e}}}$([Ne iii]) and ${n}_{{\rm{e}}}$([O iii]). In the end, ${I}_{{\rm{R}}}$([O iii] λ4363) turns out to be 0.73% of the observed I([O iii] λ4363). After we subtract ${I}_{{\rm{R}}}$([O iii] λ4363) from the observed I([O iii] λ4363), we obtain ${T}_{{\rm{e}}}$([O iii]) by adopting ${n}_{{\rm{e}}}$([O iii]).

${I}_{{\rm{R}}}$([O ii] λ7320/30) is calculated with

$$\begin{eqnarray}&&\displaystyle \frac{{I}_{{\rm{R}}}([{\rm{O}}\,{\rm{II}}]\,7320/30\,\mathring{\rm A} )}{I({\rm{H}}\beta )}=9.36{\left(\displaystyle \frac{{T}_{{\rm{e}}}}{{10}^{4}}\right)}^{0.44}\displaystyle \frac{{{\rm{O}}}^{2+}}{{{\rm{H}}}^{+}}\end{eqnarray} \tag{ 3 }$$

(Equation (2) in Liu et al. 2000), where we adopt the O²⁺/H⁺ ratio ($2.78(-4)$; see Section 3.1.2) derived from the I([O iii] 88.3 μm)/I(Hβ) ratio assuming ${T}_{{\rm{e}}}$([O iii]) and ${n}_{{\rm{e}}}$([O iii]). ${I}_{{\rm{R}}}$([O ii] λ7320/30) turns out to be 2.19% of the observed I([O ii] λ7320/30). After we subtract the recombination contribution from the observed I([O ii] λ7320/30), we obtain ${T}_{{\rm{e}}}$([O ii]) by adopting ${n}_{{\rm{e}}}$ = 260 cm⁻³.

Finally, we estimate ${I}_{{\rm{R}}}$([N ii] λ5755) using

$$\begin{eqnarray}&&\displaystyle \frac{{I}_{{\rm{R}}}([{\rm{N}}\,{\rm{II}}]\,5755\,\mathring{\rm A} )}{I({\rm{H}}\beta )}=3.19{\left(\displaystyle \frac{{T}_{{\rm{e}}}}{{10}^{4}}\right)}^{0.33}\displaystyle \frac{{{\rm{N}}}^{2+}}{{{\rm{H}}}^{+}}\end{eqnarray} \tag{ 4 }$$

(Equation (1) in Liu et al. 2000), where the N²⁺/H⁺ ratio (7.01(−5); see Section 3.1.2) was calculated using the I([N iii] 57 μm)/I(Hβ) ratio assuming ${T}_{{\rm{e}}}$([O iii]) and ${n}_{{\rm{e}}}$([O iii]). ${I}_{{\rm{R}}}$([N ii] λ5755) is 0.54% of the observed I([N ii] λ5755). After we subtract ${I}_{{\rm{R}}}$([N ii] λ5755), we obtain ${T}_{{\rm{e}}}$([N ii]) (ID: 9) by adopting ${n}_{{\rm{e}}}$ = 260 cm⁻³. In Section 4 below, we verify the above estimates of the RL contributions based on the best-fit modeling results.

We also determine ${T}_{{\rm{e}}}$(He i), which is necessary to estimate He⁺ and He²⁺ abundances, using the He i I(7281 Å)/I(6678 Å) ratio with the He i recombination coefficients in the case of ${n}_{{\rm{e}}}$ = 100 cm⁻³ given by Benjamin et al. (1999). The ${n}_{{\rm{e}}}$ and ${T}_{{\rm{e}}}$ pairs derived and adopted from the present plasma diagnostics are summarized in Table 13.

We calculate CEL ionic abundances by solving the equation of population at multiple energy levels with the adopted ${n}_{{\rm{e}}}$ and ${T}_{{\rm{e}}}$ (Table 13, which also lists the adopted ${n}_{{\rm{e}}}$ and ${T}_{{\rm{e}}}$ to calculate RL ${\mathrm{He}}^{+,2+}$ and ${{\rm{C}}}^{2+}$ abundances); the resulting ionic abundances are listed in Table 14. We give the 1σ uncertainty for each ionic abundance estimate, which is propagated from 1σ uncertainties of line fluxes, c(Hβ), ${n}_{{\rm{e}}}$, and ${T}_{{\rm{e}}}$. Ionic abundances are derived for each of the detected line intensities when more than one line 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 ${n}_{{\rm{e}}}$–${T}_{{\rm{e}}}$ 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 (Figure 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 7.5 ± 4.8). Hence, if we were to assume O⁰/H⁺ = (5.38 ± 1.05)(−4) based solely on the [O i] 145 μm line, we would have N⁰/H⁺ = (3.69 ± 2.34)(−4) and ${{\rm{S}}}^{+}$/H⁺ = (8.97 ± 6.09)(−6) by adopting a factor of 7.5 ±4.8. Nevertheless, for the O⁰/H⁺ abundance we adopt the average of the observed two optical (6300 and 6364 Å) and single far-IR (145 μm) 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 1.08(−1), I(He i λ4712) and I([Ar iv] λ4711) have to be 0.47 ± 0.18 and 0.87 ± 0.27, respectively.¹⁵ The Ar³⁺ abundance derived from this expected I([Ar iv] λ4711) is $1.99(-7)\pm 6.41(-8)$, which is consistent with the Ar³⁺ abundance derived from I([Ar iv] λ4740).

To derive the RL ${{\rm{C}}}^{2+}$ abundance, we use the C ii λ4267 line with its effective recombination coefficient in Case B for ${n}_{{\rm{e}}}$ = 10⁴ cm⁻³ defined as a polynomial function of ${T}_{{\rm{e}}}$ by Davey et al. (2000). This is justified because while the effective recombination coefficient is not available for the case of ${n}_{{\rm{e}}}$ = 100 cm⁻² that is more appropriate here, the RL abundances are in general insensitive to ${n}_{{\rm{e}}}$ for ≲10⁸ cm⁻³. As for ${T}_{{\rm{e}}}$, we adopt ${T}_{{\rm{e}}}$([Ar iii]) because the ionization potential (I.P.) of ${{\rm{C}}}^{2+}$ 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 ${\mathrm{Ne}}^{+,2+}$, ${{\rm{S}}}^{2+}$, Si⁺, Cl³⁺, and Ar²⁺) made with Spitzer and Herschel observations, more robust adaptation of ${n}_{{\rm{e}}}$ and ${T}_{{\rm{e}}}$ 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).

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 Table 14, the ICF(X) value of the element "X" and the resulting elemental abundance, ${\rm{X}}/{\rm{H}}=\mathrm{ICF}({\rm{X}})\,\times \,({{\rm{\Sigma }}}_{{\rm{m}}=1}$ ${{\rm{X}}}^{{\rm{m}}+}$/${{\rm{H}}}^{+})$, are listed in bold at the last line for each element. Here, we exclude ${{\rm{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 ${\mathrm{log}}_{10}({\rm{X}}/{\rm{H}})+12$, where ${\mathrm{log}}_{10}({\rm{H}})=12$ (in column (2)), and the relative solar abundances (X/H; in column (3)).

Table 3.  Elemental Abundances $\epsilon ({\rm{X}})$ of NGC 6781 Derived in the Present Analysis, Compared with the Solar Abundances (Column (3); $[{\rm{X}}/{\rm{H}}]=\epsilon ({\rm{X}})-{\epsilon }_{\odot }({\rm{X}})$, Where ${\epsilon }_{\odot }({\rm{X}})$ 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 Section 3.1.5)

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 $\epsilon ({\rm{X}})={\mathrm{log}}_{10}({\rm{X}}/{\rm{H}})+12$, where ${\mathrm{log}}_{10}({\rm{H}})=12$. The CEL C abundance, C(CEL), is an expected value.

Download table as:  ASCIITypeset image

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 Section 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 ${{\rm{N}}}^{+,2+,3+}$ and adopt ICF(N) ≈ ICF(O), which is equal to the ${\rm{O}}/({{\rm{O}}}^{+}+{{\rm{O}}}^{2+}$) ratio. Similarly, we assume that the Si abundance is the sum of ${\mathrm{Si}}^{+,2+,3+}$ 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 ${\mathrm{Cl}}^{+,2+,3+}$ and ${\mathrm{Ar}}^{+,2+,3+}$, respectively. Then, we adopt ICF(Cl) ≈ICF(Ar) ≈ICF(S), which corresponds to the ${\rm{S}}/({{\rm{S}}}^{2+}+{{\rm{S}}}^{3+}$) 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 $4.06(-3)\pm 1.19(-3)$. Note that we do not include the CEL ${{\rm{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 ${{\rm{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 (Z) of the object is close to the solar metallicity (see also Section 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 solar metallicity ($Z\sim 0.02$ Z_☉). Hence, the derived RL C abundance of $4.06(-3)\pm 1.19(-3)$ 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 ${{\rm{C}}}^{2+}$(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 ${{\rm{C}}}^{2+}$(RL)/C²⁺(CEL) ratio of 4.41 ± 0.81 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 ${{\rm{C}}}^{2+}$(RL)/C²⁺(CEL) ratio of 4.10 ± 0.49 found among 58 PNe in the Milky Way and Magellanic Clouds (Otsuka et al. 2011) to obtain the scaled expected CEL C of $9.89(-4)\pm 3.14(-4)$.

This expected CEL C of $9.89(-4)\pm 3.14(-4)$ ($\epsilon ({\rm{C}})=9.00$) would be more reasonable than the measured RL C abundance of $4.06(-3)\pm 1.19(-3)$ 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 ${{\rm{C}}}^{2+}$(RL)/C²⁺(CEL) ratio of 3.63 for NGC 6720, which possesses central star and nebula properties very similar to those of NGC 6781 (see Section 3.4.1).

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 ${\mathrm{Cl}}^{2+,3+}$ abundances of $1.07(-7)$ and $1.57(-8)$ and the corresponding ICF(Cl) value of 1.17, while Liu et al. (2004a) used the Cl²⁺ abundance of $7.92(-8)$ only with the corresponding ICF(Cl) of 3.394. In addition, the adopted ${T}_{{\rm{e}}}$ could contribute to the discrepancy because the Cl ionic abundances are determined using their CEL lines, whose emissivities are sensitive to ${T}_{{\rm{e}}}$. Overall, we would argue again that our (Cl) value is more improved than the previous estimate because we have more robust ${T}_{{\rm{e}}}$ for the ionic Cl abundances and we derive a Cl³⁺ abundance that would reduce uncertainties in ICF(Cl).

The discrepancy in RL (C) is due to different values of I(C ii λ4267) (which might be caused by different adopted c(Hβ)) and adopted ICF(C): our (C) and ICF(C) values are $2.0(-3)$ and 2.03, whereas theirs are $9.05(-4)$ 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 9.89(−4) ± 3.14(−4) ((C) = 9.00) adopted in the previous section, in comparison with the observed O(CEL) of $5.81(-4)\pm 2.19(-5)$ ((O) = 8.76), would suggest a slightly C-rich nature for NGC 6781 (C/O number density ratio of 1.70 ± 0.54). Indeed, the Spitzer/IRS mid-IR spectrum (Figure 1, inset) shows polycyclic aromatic hydrocarbon (PAH) emission at 6–9 μ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 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 ${{\rm{C}}}^{2+}$ of $9.05(-4)$ 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 ${{\rm{C}}}^{2+}$(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 NGC 6781 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 I(C ii λ4267) measurements to the ICF(C) and ${{\rm{C}}}^{2+}$(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 $\epsilon ({\rm{C}})=8.56\mbox{--}9.00$ (correspondingly, $[{\rm{C}}/{\rm{H}}]=0.17\mbox{--}0.61$) based on the arguments presented above.

We compare the derived $\epsilon ({\rm{X}})$ with the values predicted by AGB nucleosynthesis models. As for the metallicity Z 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 Z indicators. The average between the observed [Cl/H] and [Ar/H] values of −0.05 corresponds to $Z\sim 0.018$.

However, the S abundance ([S/H] = −0.25) suggests a much lower Z. So far, this S abundance anomaly has been found in many Milky Way and M31 PNe (Henry et al. 2012; see their Figure 1). Henry et al. (2012) concluded that the sulfur deficit in PNe is generally reduced by increasing the ${{\rm{S}}}^{3+}$ 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(α, n)²⁵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 ${{\rm{S}}}^{+}$/${{\rm{H}}}^{+}=8.97(-6)$ in Section 3.1.2, we would obtain $\epsilon ({\rm{S}})=7.20$, which is consistent with ${({\rm{S}})}_{\odot }$).

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

Meanwhile, adopting the upper CEL abundance limit of $\epsilon ({\rm{C}})=9.00$, the ${\chi }^{2}$ values suggest that the observed $\epsilon ({\rm{X}})$ is most consistent with the 2.5 M_☉ model (reduced ${\chi }^{2}=16.15$). The reduced ${\chi }^{2}$ value = 17.5 of the 3.0 M_☉ 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 M_☉.

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. Figure 4 shows the H₂ $v=1-0$ 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 (see Figure 2). Figure 4(a) 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; Figure 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 (Figure 3 of HerPlaNS1). Highly localized distributions of the molecular gas component are apparent from the filamentary appearance of the H₂ emission (Figure 4(b)). These H₂ filaments (and maybe clumps, too) are patches of H₂ that survived in the ionized region.

Zoom In Zoom Out Reset image size

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 (Figure 1, inset). Observations made by Arias & Rosado (2002) show that the intensity of H₂ $v=2-1$ S(1) at 2.248 μm is much fainter than that of H₂ $v=1-0$ 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₂ I(2.248 μm)/I(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 the H₂ I(2.248 μm)/I(2.122 μm) ratio (e.g., Otsuka et al. 2013), even though it is not easy to do with K-band data alone.

Interestingly, the expansion velocity of H₂ (∼22 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₂ $v=1-0$ S(1) surface brightness could be explained by shocks at 10–24 km s⁻¹ heading into the pre-shock region of the H₂ density at 3400–14,900 cm⁻³.

We investigate the conditions in the H₂-emitting regions by comparing the flux ratios of mid-IR H₂ lines to the $v=0-0$ S(3) line at 9.67 μm with the theoretical continuous shock (C-shock) models by Flower & Pineau (2010). The observed I(17.04 μm)/I(9.67 μm) ratio suggests a match for a model with the shock velocity of ${V}_{s}=10$ km s⁻¹ and pre-shock hydrogen density of ${n}_{s}({\rm{H}})={\rm{200,000}}$ cm⁻³, while the observed I(12.29, 8.02, 6.91, 6.11, 5.51 μm)/I(9.67 μm) ratios point to a model with ${V}_{s}=20$ km s⁻¹ and ${n}_{s}({\rm{H}})={\rm{20,000}}$ 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 I(H i 12.3 μm) = 0.971 when I(Hβ) = 100 in the case of ${T}_{{\rm{e}}}$ = 10⁴ K and ${n}_{{\rm{e}}}$ = 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 3D structure modeling using CO velocity maps. A model by Flower & Pineau (2010) with a shock velocity of ${V}_{s}=30$ km s⁻¹ and pre-shock hydrogen density of ${n}_{s}({\rm{H}})={\rm{20,000}}$ cm⁻³ would explain the observed far-IR CO line flux ratios with respect to the CO $J=7-6$ 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 10–30 km s⁻¹ impinging onto the pre-shock region at ${n}_{s}({\rm{H}})\sim {\rm{20,000}}\mbox{--}{\rm{200,000}}$ cm⁻³. These shocks may 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 Section 4.2.8). Together with the filamentary/clumpy appearance of the H₂ emission regions (Figure 4), we would conclude that these structures represent high-density regions delineating the locations of thermal collisional excitation embedded in a 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, and references 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 $\leqslant 10$ 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 need to be resolved in the future, we tentatively conclude that the observed H₂ emission in NGC 6781 has a shock origin based on the arguments presented above.

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 N_u in the upper state is written as

$$\begin{eqnarray}&&{N}_{u}=\displaystyle \frac{4\pi I({{\rm{H}}}_{2})}{A}\cdot \displaystyle \frac{\lambda }{{hc}},\end{eqnarray} \tag{ 5 }$$

where I(H₂) is the H₂ line intensity in erg s⁻¹ cm⁻² sr⁻¹, A is the transition probability taken from Turner et al. (1977), h is the Planck constant, and c is the speed of light. In LTE, the Boltzmann equation relates N_u to the excitation temperature T(H₂) via

$$\begin{eqnarray}&&\mathrm{ln}\,\left(\displaystyle \frac{{N}_{u}}{{g}_{u}}\right)=-\displaystyle \frac{{E}_{u}}{{kT}({{\rm{H}}}_{2})}+\mathrm{ln}\,\left[N({{\rm{H}}}_{2})\cdot \displaystyle \frac{{hcB}}{2{kT}({{\rm{H}}}_{2})}\right],\end{eqnarray} \tag{ 6 }$$

where g_u is the vibrational degeneracy, E_u is the energy of the excited level taken from Dabrowski (1984), k is the Boltzmann constant, and B is the rotational constant (60.81 cm⁻¹).

In Figure 5, we plot the $\mathrm{ln}({N}_{u}/{g}_{u})$ versus ${E}_{u}/k$ for each of the H₂ lines detected in NGC 6781 (Table 4). The ${E}_{u}/k$ of the H₂ $v=2-1$ S(1) (magenta circle) is calculated using the average line intensity of $2.7(-4)$ 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.

Zoom In Zoom Out Reset image size

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 (Figure 5(a)): $T({{\rm{H}}}_{2})=1279\pm 109$ K and $N({{\rm{H}}}_{2})=(2.28\pm 0.49)(18)$ cm⁻². The derived $T({{\rm{H}}}_{2}$) is comparable to 978 K and 880 ± 70 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 and 17.04 μm lines).

Next, we fit all H₂ lines (including 17.04 μm) using two excitation temperatures (Figure 5(b)). The warm component is found to have $T({{\rm{H}}}_{2})=1161\pm 72$ K and $N({{\rm{H}}}_{2})\,=(2.72\pm 0.53)(18)$ cm⁻², whereas the cold component is found to have $T({{\rm{H}}}_{2})=236\pm 50$ K and $N({{\rm{H}}}_{2})=(6.67\pm 4.89)(19)$ 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.

To conclude this subsection, we estimate the mass of the molecular gas component in the nebula by adopting the distance of 0.46 kpc (Section 3.4.1). Based on the H₂ and CO emission maps (Hiriart 2005; 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 ($N({{\rm{H}}}_{2})=2.72(18)$ cm⁻² and $6.67(19)$ cm⁻², respectively), we estimate the H₂ gas mass of $2.5(-3)\,{M}_{\odot }$ and $6.2(-2)\,{M}_{\odot }$ for the warm and cold components, respectively.

Previously, we derived $N(\mathrm{CO})={10}^{14.70\mbox{--}15.08}$ cm⁻² (excitation temperature at ∼60 K) based solely on our Herschel spectra (HerPlaNS1). Bachiller et al. (1997) measured $N(\mathrm{CO})={10}^{16.16}$ cm⁻² (excitation temperature at ∼25 K) based on submillimeter data. Assuming that each of the above $N(\mathrm{CO})$ 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 $4.6(-6)\,{M}_{\odot }$ and $6.6(-5)\,{M}_{\odot }$, respectively. These estimates are combined to yield the total molecular gas mass (of H₂ and CO) of $6.46(-2)\,{M}_{\odot }$.

The empirical N(CO)/N(H₂) ratio turns out to be $2.19(-4)$ and $4.37(-4)$ for the warm and cold temperature regions, respectively. Assuming that the N(CO)/N(H₂) ratio translates roughly to $2\times n({\rm{C}})/n({\rm{H}})$, we can estimate (C) of 8.04–8.34 for the molecular gas component. Compared with the adopted CEL expected (C) of 8.56–9.00 for the ionized gas component, ∼11%–60% of the C atoms were estimated to be locked in as molecules.

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 (Figure 1, Table 11) with the post-AGB evolutionary tracks produced by Vassiliadis & Wood (1994) augmented with a grid of synthesized spectra by Rauch (2003). Although several new evolutionary tracks have been produced since then, there have 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 L_* using a grid of non-LTE line-blanketed plane-parallel hydrostatic atmospheric models generated by Rauch (2003) as templates. We adopt the solar abundance (Z = 0.02) models for the CSPN based on the results of our own nebular abundance analysis presented in Section 3.1.5.

To characterize the stellar atmosphere fully, we also need the effective temperature ${T}_{\mathrm{eff}}$ and surface gravity $\mathrm{log}\,g$ of the CSPN. Previously, Rauch et al. (2004) suggested ${T}_{\mathrm{eff}}={\rm{80,000}}$ K and $\mathrm{log}\,g=6.0$ cm s⁻² based on the stellar absorption line fitting. If this ${T}_{\mathrm{eff}}$ were true, the CSPN would have been still burning hydrogen in a thin surface layer while increasing ${T}_{\mathrm{eff}}$. However, detection of strong He ii λ4686 and [O iv] 25.88 μm lines in the ISIS and Spitzer/IRS spectra, respectively (Figure 1 and Table 12) requires ${T}_{\mathrm{eff}}\gt {\rm{80,000}}$ 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 ${T}_{\mathrm{eff}}$ and $\mathrm{log}\,g$ values in a PN similar to NGC 6781 in terms of nebula and CSPN properties. Among 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).

Table 5.  Similarities between NGC 6781 and NGC 6720

PNe	(He)	(${{\rm{C}}}_{\mathrm{RL}}$)	(${{\rm{C}}}_{\mathrm{CEL}}$)	(N)	(${{\rm{O}}}_{\mathrm{RL}}$)	(${{\rm{O}}}_{\mathrm{CEL}}$)	(Ne)	(S)	(Cl)	(Ar)	${T}_{\mathrm{eff}}$ (kK)	$\mathrm{log}\,g$ (cm s−2)	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 Section 3.1.3); (2) Schwarz & Monteiro 2006, for ${T}_{\mathrm{eff}}$ and $\mathrm{log}\,g$ via photoionization modeling; (3) Rauch et al. 2004, for ${T}_{\mathrm{eff}}$ and $\mathrm{log}\,g$ via stellar absorption fitting; (4) Liu et al. 2004a, for abundances; (5) McCarthy et al. 1997, for ${T}_{\mathrm{eff}}$ and $\mathrm{log}\,g$ via stellar absorption fitting; (6) Napiwotzki 1999, for ${T}_{\mathrm{eff}}$ and $\mathrm{log}\,g$ via stellar absorption fitting; (7) van Hoof et al. 2010, for ${T}_{\mathrm{eff}}$ via Cloudy photoionization modeling.

Download table as:  ASCIITypeset image

The CSPN of NGC 6720 has a ${T}_{\mathrm{eff}}\gt 100$ 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 ${T}_{\mathrm{eff}}=110\mbox{--}140$ kK and $\mathrm{log}\,g=6.9$ 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; see also Section 4).

Then, we scale the synthesized Rauch model spectra of the adopted CSPN characteristics of ${T}_{\mathrm{eff}}=110\mbox{--}140$ kK with a constant 10 kK step with $\mathrm{log}\,g=6.9$ cm s⁻² fixed so that the observed photometry from the WFC u band to WFCAM K band (see Table 11) matches with the model spectra (Figure 6, showing the ${T}_{\mathrm{eff}}=120$ kK case). The scaled spectra are integrated to yield L_*, which is then parameterized with ${T}_{\mathrm{eff}}$ and the distance D in the form of

$$\begin{eqnarray}&&{L}_{* }(D,{T}_{\mathrm{eff}})\\ &&=\,[2.29(-7)\cdot {T}_{\mathrm{eff}}^{2}-4.39(-2)\cdot {T}_{\mathrm{eff}}+2510]\cdot {D}^{2},\end{eqnarray} \tag{ 7 }$$

where D is in kpc and ${T}_{\mathrm{eff}}$ is in K. Note that L_* is not very sensitive to $\mathrm{log}\,g$. For instance, L_* increases only by ∼0.8% when $\mathrm{log}\,g$ is reduced from the adopted 6.9 cm s⁻² to 6.6 cm s⁻². Thus, our choice of single $\mathrm{log}\,g$ value is warranted.

Zoom In Zoom Out Reset image size

Finally, we compute ${L}_{* }(D,{T}_{\mathrm{eff}})$ at ${T}_{\mathrm{eff}}$ = 110–140 kK and for a range of D and plot the resulting (L_*, ${T}_{\mathrm{eff}}$) pairs over the post-AGB evolutionary tracks of the 1.5, 2.25, 2.5, and 3.0 M_☉ initial-mass stars produced by Vassiliadis & Wood (1994), as shown in Figure 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 M_☉ (Section 3.1.5). Also, the previous analysis by Schwarz & Monteiro (2006) suggested a CSPN initial mass of 1.5 M_☉.

Zoom In Zoom Out Reset image size

We find that $D=0.34\mbox{--}0.52\,\mathrm{kpc}$ fits the initially $2.25\mbox{--}3.0\,{M}_{\odot }$ post-AGB evolutionary tracks the best for the adopted ${T}_{\mathrm{eff}}$ range (light-blue box in Figure 7). Therefore, we adopt $D=0.46\,\mathrm{kpc}$, which is the the intermediate value between 0.34 and 0.52 kpc (red circles in Figure 7). Accordingly, we find ${L}_{* }=104\mbox{--}196$ L_☉. 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) ($D=0.46\,\mathrm{kpc}$, using the post-AGB evolutionary tracks for 2.0 and 3.0 M_☉ stars with $Z=0.02;$ see also Figure 12).

Previously, Schwarz & Monteiro (2006) concluded that the progenitor of NGC 6781 was a 1.5 ± 0.5 M_☉ initial-mass star based on their derived L_* and ${T}_{\mathrm{eff}}$ values, provided that $D=0.95\,\mathrm{kpc}$ as suggested from their photoionization model fitting (black triangle in Figure 7; also suggesting that NGC 6781 was in cooling phase). At $D=0.95\,\mathrm{kpc}$, our L_* estimates would be consistent with the 1.5 M_☉ evolutionary track (blue squares in Figure 7). However, the progenitor CSPN mass of NGC 6781 would most likely exceed 1.5 M_☉ because of its empirically determined elemental abundances (Section 3.1.3) and H₂ detection in this object (Section 3.2.2).

With a survey of H₂ $v=1-0$ 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 ($\geqslant 1.5$ M_☉) 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 2.25–3.0 M_☉ initial mass for the CSPN of NGC 6781 and the distance of 0.46 kpc based on the ${L}_{* }(D,{T}_{\mathrm{eff}})$ fitting.

The filamentary appearance of the nebula (Figure 4) and low ${n}_{{\rm{e}}}$ even in the central ionized regions (Section 3.1.1) are also suggestive that NGC 6781 is a highly evolved PN. Referring back to the similarity to NGC 6720, comparisons between L_* and ${T}_{\mathrm{eff}}$, where L_* 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 M_☉ stars of Z = 0.02, 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 $\gtrsim 3.6\,\mathrm{kpc}$. According to Vassiliadis & Wood (1994), L_* is nearly constant at ∼6300 L_☉ along the horizontal part of the post-AGB track for a 2.5 M_☉ initial-mass star with Z = 0.02. In this case, the number of the ionizing photons is 4.25(+47)s⁻¹ for ${T}_{\mathrm{eff}}={\rm{120,000}}$ K and $\mathrm{log}\,g=6.9$ cm s⁻¹. The Strömgren radius for this radiation field in a constant hydrogen density of 300 cm⁻³ (see Table 2, Figure 8) with a filling factor (f) of unity would be ∼0.41 pc. This corresponds to the apparent radius of 237 at $D=3.6\,\mathrm{kpc}$, which disagrees with the observed ionization radius of $\sim {55}^{{\prime\prime} }$. Because the Strömgren radius is proportional to ${f}^{-1/3}$, it would be consistent with the observed ionization radius at $D=3.6\,\mathrm{kpc}$ if f were 0.12. However, according to the empirical method introduced by Mallik & Peimbert (1988), the f value of NGC 6781 is estimated to be ∼0.4 at $D=3.6\,\mathrm{kpc}$ and almost unity at 0.46 kpc. Therefore, we conclude that the CSPN of NGC 6781 already evolved off to the cooling track currently with ${L}_{* }=104\mbox{--}196\,{L}_{\odot }$ and ${T}_{\mathrm{eff}}=110\mbox{--}140$ kK at $D=0.46\,\mathrm{kpc}$.

Zoom In Zoom Out Reset image size

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 (V − I and I − J) and the grid of theoretical color indices as a function of ${T}_{\mathrm{eff}}$, Douchin et al. (2015) argued that CSPN of NGC 6781 shows near-IR excess owing to an M1- to M7-type companion star. However, we do not observe any IR excess in the SED of the CSPN (Figures 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 c(Hβ) = 1.007, the extinction-corrected V–I and I–J colors of the CSPN were −0.44 and −0.90, respectively. If we used $E(B-V)=0.56$ (corresponding to c(Hβ) = 0.82) as adopted by Douchin et al. (2015), the respective V–I and I–J colors would become redder, −0.27 and −0.61, which would be in perfect agreement with Douchin et al. (2015). This would negate the necessity for an M1- to 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 $D=0.46\,\mathrm{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.

As the incident SED from the CSPN, we adopt the theoretical atmospheric model grid by Rauch (2003) for a star with Z = 0.02 and $\mathrm{log}\,g=6.9$ cm s⁻² (see Figure 6 for the case of Z = 0.02, $\mathrm{log}\,g=6.9$ cm s⁻², and ${T}_{\mathrm{eff}}=120$ kK). We keep the distance of 0.46 kpc to NGC 6781 and vary ${T}_{\mathrm{eff}}$ and L_* within the possible ranges, ${L}_{* }=104\mbox{--}196\,{L}_{\odot }$ and ${T}_{\mathrm{eff}}=110\mbox{--}140$ kK, as determined in Section 3.4, during the iterative model fitting to search for the best-fit model parameters that would reproduce the observational data.

For the elemental abundances of the nebula, we adopt the empirically determined abundances (Table 3; Section 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 Section 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 $\epsilon ({\rm{C}})$ of 8.56–9.00 (Table 3). Hence, for the purpose of the present modeling, we set $\epsilon ({\rm{C}})$ 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 M_☉ (see Section 3.1.5).

The expected CEL $\epsilon ({\rm{C}})$ of 8.56 is also consistent with the AGB nucleosynthesis model for the $2.25\,{M}_{\odot }$ stars (Karakas 2010). As we demonstrated in Section 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 $\epsilon ({\rm{C}})$ 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 M_☉ initial-mass star with Z = 0.02 (Karakas 2010). However, the Fe abundance is another exception, because we overpredict the Fe lines when setting $\epsilon (\mathrm{Fe})=7.53$ as determined by Karakas (2010). The model I([Fe ii] 17.9 μm) and I([Fe iii] λ4880) line fluxes turn out to be 31.2 and 2.6 (with respect to $I($ Hβ$)=100$), 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 $\epsilon (\mathrm{Fe})=6.20$ in NGC 6720. Thus, we adopt $\epsilon (\mathrm{Fe})=6.20$, following the same similarity argument between NGC 6781 and NGC 6720 as in Section 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.

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 1D 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 temperatures by the projection effect (Figure 4(a)). While clumps/filaments of H₂ surviving in the ionized region would be plausible (Figure 4(b)), we simply adopt this "onion skin" configuration for the sake of 1D 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 1D 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 J = 3−2 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). Figure 8 shows a schematic of the adopted geometry.

The input radial hydrogen density profile, ${n}_{{\rm{H}}}(R)$ (where R is the distance from the CSPN), is adopted from our previous analysis (HerPlaNS1). In the central cavity surrounded by the barrel wall structure ($R\lt 54^{\prime\prime}$) ${n}_{{\rm{H}}}(R)=300$ cm⁻³, whereas in the barrel wall ($54^{\prime\prime} \leqslant R\lt 58^{\prime\prime}$) ${n}_{{\rm{H}}}(R)=960$ cm⁻³ (Figure 8).

Unfortunately, ${n}_{{\rm{H}}}(R)$ beyond $58^{\prime\prime}$ cannot be determined directly from the observed data, because this radial region is where the surface brightness of the nebula decreases sharply to the detection limit in the narrow- and broadband images of the object (and hence the observational constraints are scarce). Hence, as discussed in Section 3.2.3, we simply adopt a constant density of ${n}_{{\rm{H}}}(R)={10}^{4}$ cm⁻³ beyond $58^{\prime\prime}$. The outer radius is then determined iteratively by increasing the thickness of this dusty PDR layer until the model flux at 170 μm would reproduce the observed value, which is one of our model calculation termination criteria. In the end, the outer radius is set to $61^{\prime\prime}$. The radial density profile of the nebula is also provided in Figure 8.

One might surmise that the adopted ${n}_{{\rm{H}}}(R)$ radial profile would allow for a constant gas pressure model. Therefore, we test a constant gas pressure model, for which we adopt the average log₁₀(${T}_{{\rm{e}}}$ ${n}_{{\rm{e}}}$) = 6.81 K cm⁻³ based on the radial ${T}_{{\rm{e}}}$ and ${n}_{{\rm{e}}}$ profiles measured previously (HerPlaNS1). The result is similar to the nonconstant gas pressure model, except for He ii and [O iv] lines. In order to avoid a collapse of the nebula, the inner radius of the nebula has to be set larger. This correspondingly results in underestimates of the line fluxes of these high I.P. lines. Also, NGC 6781 does not seem to be embedded in a dense ISM region. Because of these reasons, we conclude that the nonconstant gas pressure model that we adopt in the present investigation is a better approximation to NGC 6781 than a constant gas pressure model.

As we summarized in Section 3.3, NGC 6781 is determined to be a PN rich in amorphous carbon. Thus, the nebula's dusty PDR is expected to consist largely of amorphous carbon (AC) plus neutral (and possibly ionized) PAHs, even though the C-richness of the nebula remains uncertain (see Section 3.1.4). Rouleau & Martin (1991) provided two types of optical constants measured from samples "BE" (soot produced from benzene burned in air) and "AC" (soot produced by striking an arc between two amorphous carbon electrodes in a controlled Ar atmosphere). We test both of these BE and AC amorphous carbon grain models, and we find that the AC-type grain models yield generally better overall fit to the observed mid-IR to far-IR dust continuum. Thus, we adopt the AC-type grain optical constants by Rouleau & Martin (1991). We assume spherical grains and adopt the modified interstellar size distribution (i.e., $n(a)\propto {a}^{-3.5}$; Mathis et al. 1977) with $a=0.005\mbox{--}0.50$ μm, which are divided into 20 bins in model calculations.

For PAHs, we adopt the radius a in the range of 0.0004 μm (30 C atoms) to 0.0081 μm (250 C atoms) with the same size distribution as dust (${a}^{-3.5}$; Mathis et al. 1977), approximating the overall shape by a sphere (separated into the same 20 size bins). We include both the neutral and charged PAH grains. The optical constants for PAH-carbonaceous grains are adopted from the theoretical work by Draine & Li (2007). We permit the stochastic heating mechanism of PAH molecules in model calculations.

Figure 9 shows the SED of the CSPN plus PN based on the observed photometry from GALEX 0.22 μm to radio 1.4 GHz (Table 1; Figure 1). Using this empirical SED, we measure the integrated luminosity of 114 L_☉ at $D=0.46\,\mathrm{kpc}$ for the CSPN plus PN. The contribution to this SED only from the CSPN for the wavelength range of $\gtrsim 0.2$ μm is estimated to be 4.6 L_☉. Hence, the remainder has to come from the nebula, i.e., ${L}_{\mathrm{Neb}}\approx 110\,{L}_{\odot }$.

Zoom In Zoom Out Reset image size

As for the luminosity of the CSPN, we already determined the empirical value of ${L}_{* }=104\mbox{--}196\,{L}_{\odot }$ based on Equation (7) (Section 3.4.1). Thus, NGC 6781 could be a density-bounded PN (i.e., ${L}_{\mathrm{Neb}}\lt {L}_{* }$) as previously claimed by Schwarz & Monteiro (2006). However, the fact that NGC 6781 possesses massive molecular gas and dust components indicates that it is more likely an ionization-bounded PN (i.e., ${L}_{\mathrm{Neb}}\approx {L}_{* }$). Realistically speaking, whether a PN is density- or ionization-bounded is not necessarily straightforward, because both situations could be present in one PN. In bipolar PNe such as NGC 6781, both ionization- and density-bounded conditions are expected to be present in the nebula along the equatorial and polar directions, respectively.

Based on the resemblance between the observed spatial distribution of the ionized gas and that of the other (molecules and dust) components (Figure 4; Zuckerman et al. 1990; Hiriart 2005; HerPlaNS1), the transition from the ionized region to the PDR must be happening quite rapidly over a small radial range. Hence, we start model calculations with a nebula that is ionization bounded at around $R=55^{\prime\prime}$, which corresponds to the outer radius of the central ring structure of the nebula and also the intensity peak of H₂ and CO emission (see Section 3.2.4). The use of the cylinder option is also corroborated by the density-bounded nature of the nebula expected in the polar directions of the nebula.

We introduce a high-density PDR wall beyond the ionization front in the model geometry (Figure 8) to explain the observed molecular emission. However, this causes significant underestimates of the observed H₂ and high-J CO line fluxes, as well as their column densities. This failure suggests the presence of an additional heating source in the PDR.

An obvious extra PDR heating source is the interstellar radiation field (ISRF). However, no meaningful heating of the PDR can be achieved by the ISRF in the present model for NGC 6781: only ≲1% of the observed H₂ flux is reproduced by the nominal Galactic ISRF. Hence, it is unrealistic to expect to generate enough heating to reproduce all of the observed H₂ flux by the ISRF alone unless it is unrealistically enhanced. Thus, it is reasonable to expect something other than the ISRF for a PDR heating source to explain the observed H₂ fluxes. By the same token, the Galactic background cosmic ray cannot possibly work as a PDR heating source unless it is unrealistically enhanced.

Soft X-rays.—Another extra heating source is soft X-ray emission from a high-temperature CSPN as suggested by the presence of PNe in which X-rays were detected (e.g., Chu et al. 2001; Kastner et al. 2012; Montez et al. 2015). Soft X-rays (50 ev–10 keV) from a CSPN of ${T}_{\mathrm{eff}}\gt 100$ kK can strengthen H₂ line emission, because such high-energy photons would penetrate into the PDR beyond the ionization front (Natta & Hollenbach 1998). Using data from the Chandra X-ray Observatory, Montez et al. (2015) examined the X-ray luminosities for a group of Galactic PNe including NGC 6781. They found that no X-rays were detected from NGC 6781 in the $0.3\mbox{--}8.0\,\mathrm{keV}$ energy band, while a simple blackbody of ${T}_{\mathrm{eff}}\sim 120\mbox{--}130$ kK at 0.46 kpc is sufficient for detectable X-ray fluxes in the 0.3–8.0 keV energy band (their Figure 14). Hence, the nondetection of X-ray emission in NGC 6781 is indicative of strong interstellar extinction or metal line blanketing, either of which can suppress the X-ray emission to below the detection limit.

We examine whether X-ray emission possible from the CSPN of NGC 6781 can result in a better fit to the observed H₂ line fluxes under the following two scenarios: (1) if the X-ray luminosity (${L}_{{\rm{X}}}$) of the CSPN were to power the entire observed mid-IR H₂ luminosity ($\sim 5.59\times {10}^{33}$ erg s⁻¹ at $D=0.46\,\mathrm{kpc};$ Table 4), but were to be suppressed completely by the extinction, and (2) if the CSPN possessed an atmosphere of subsolar metallicity to circumvent metal line blanketing. The predicted H₂ line fluxes with these X-ray emission enhancements would not reproduce the observed line fluxes even if we adopted (1) an extra blackbody emitting in the range of 0.27–10.4 keV with a luminosity of $\sim {10}^{33}$ erg s⁻¹ at 10³ kK or (2) an atmosphere of Galactic halo metallicity for the CSPN. Therefore, we conclude that extra soft X-ray would not possibly produce the observed H₂ line fluxes in NGC 6781.

Shock heating in the PDR.—Yet another extra heating source is a mechanical heat input by shocks as suggested from the H₂ excitation diagram analysis (Sections 3.2.2 and 3.2.3). This idea, previously used in a study of the C-rich PN NGC 7027 by Hasegawa et al. (2000), can work to excite H₂ lines in regions far enough away from the CSPN. As Cloudy does not handle shocks, the desired extra heating by shocks is achieved by invoking the "temperature floor" option, which forces the predetermined value of the electron temperature ${T}_{{\rm{e}}}$ over a specific region (see Section 4.3). We iteratively search for the optimum floor temperature in the PDR ($R\geqslant 58^{\prime\prime}$) between 800 and 1600 K. This temperature range is suggested by the H₂ excitation temperatures derived from the excitation diagram analysis (Section 3.2.3).

While the use of a "temperature floor" helps to reproduce the observed warm H₂ lines (except for 17.04 μm), as well as high-J CO and OH lines,¹⁶ the adaptation of the "temperature floor" also introduces negative side effects such as (1) suppression of molecular lines with lower excitation temperatures and (2) overestimation of atomic gas line fluxes such as far-IR [O i] and [C ii] lines that have low-excitation energy at the upper levels. These side effects would make the mass fraction of the atomic and molecular gas with respect to the neutral (atomic + molecular) gas highly uncertain, primarily because the model would fail to account for the cold molecular component while introducing the corresponding amount of extra atomic gas (as the total amount of neutral gas was practically set by the input hydrogen density profile; Figure 8). However, the proper amount of the warm and cold molecular components, as well as the atomic gas component, can be recovered (Section 4.3.4).

To find the best-fit model, we vary 13 parameters—${T}_{\mathrm{eff}}$, L_*, the inner radius of the shell (${R}_{\mathrm{in}}$), elemental abundances ((He/N/O/Ne/Si/Cl/Ar), except for $\epsilon ({\rm{C}})$, which was fixed), dust and PAH mass fraction, and the floor temperature of the PDR—within a given range by using the optimize command available in Cloudy. We terminate iterative calculations when any one of the predicted flux densities, ${F}_{\nu }$(170 μm), ${F}_{\nu }$(250 μm), or ${F}_{\nu }$(350 μm), reaches the corresponding observed value. Practically, the terminating conditions would determine the maximum ${R}_{\mathrm{out}}$, i.e., the thickness of the dense PDR beyond the inner ionized region, by setting the amount of far-IR continuum emission. The flux densities at 170, 250, 350 μm are selected as constraints because there are no strong emission lines in these bands and they can be compared with measurements made in the PACS 160 and SPIRE 250 and 350 μm bands. In this sense, ${R}_{\mathrm{out}}$ is not a free parameter.

The best-fit model is determined by the minimum ${\chi }^{2}$ (16 for the best-fit model) calculated from the following 136 observational constraints: 37 broadband fluxes, 78 gas emission line fluxes relative to Hβ as well as I(Hβ), 19 flux densities in mid-IR, far-IR, and radio wavelengths, and the ionization boundary radius (${R}_{\mathrm{IB}}$). We define ${R}_{\mathrm{IB}}$ as the radial distance from the CSPN at which ${T}_{{\rm{e}}}$ drops below 4000 K; below such a temperature, no ionized gas emission lines except for [C ii] and [S ii] would be measurable.

In Table 6, we summarize the best-fit parameters. The SED of the best-fit model, in comparison with the observational data, is presented in Figure 10. Figure 11 is also provided to show the quality of the best-fit model with blow-ups of various wavelength ranges with major emission lines. In Table 15, we list the best-fit model versus observed quantities of gas emission line fluxes relative to Hβ, broadband fluxes relative to Hβ, and flux densities.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Table 6.  Characteristics of the Best-fit Cloudy Model of NGC 6781

Parameters of the CSPN	Values
L*/${T}_{\mathrm{eff}}$/log g	121 L☉/120 870 K/6.9 cm s−2
D	0.46 kpc
Parameters of the Nebula	Values
(X)	He:11.02, C:8.56, N:8.10, O:8.64,
	Ne:8.00, Si: 6.25, S:6.82, Cl:5.01,
	Ar:6.22, Fe:6.20
	Others: Karakas (2010)
Geometry (Figure 8)	"Cylinder" with height = 90'' (0.201 pc)
	Inner radius (${R}_{\mathrm{in}}$) = 0.52'' (0.001 pc)
	Ionization boundary (${R}_{\mathrm{IB}}$) = 55'' (0.123 pc)
	Outer radius (${R}_{\mathrm{out}}$) = 61'' (0.135 pc)
Adopted ${n}_{{\rm{H}}}$ (Figure 8)	Inner cavity ($R\lt 54^{\prime\prime}$): 300 cm−3
	Barrel wall ($54^{\prime\prime} \leqslant R\lt 58^{\prime\prime}$): 960 cm−3
	PDR ($58^{\prime\prime} \leqslant R\lt 61^{\prime\prime}$): 104 cm−3
Temperature ${T}_{{\rm{e}}}$	Inner cavity ($R\lt 54^{\prime\prime}$): 23 820–10 260 K
	Barrel wall ($54^{\prime\prime} \leqslant R\lt 58^{\prime\prime}$): 10 260–2 750 K
	PDR ($58^{\prime\prime} \leqslant R\lt 61^{\prime\prime}$): 2 750–1 420 K
Filling factor (f)	1.0
${\mathrm{log}}_{10}I$(Hβ)	−9.890 erg s−1 cm−2 (de-reddened)
Temperature floor	1420 K
Mass	Ionized gas: 0.094 M☉
	Neutral (atomic + molecular) gas: 0.31 ${M}_{\odot }$ a
Parameters of the Dust	Values
and PAHs
Particle size	PAH (neutral & ionized): 0.0004–0.011 μm,
	AC: 0.005–0.50 μm
Temperature	PAH (neutral): 71–515 K,
	PAH (ionized): 72–367 K,
	AC: 22–299 K
Mass	PAH (neutral): 3.30(−7)M☉
	PAH (ionized): 2.46(−6)M☉
	AC: 1.53(−3)M☉
GDR	268

Note.

^aWe corrected the molecular gas mass of 0.11 M_☉ and the atomic gas mass of 0.20 M_☉. See Section 4.3.4 and Table 10.

Download table as:  ASCIITypeset image

Here, we can retroactively check whether the empirical estimates and adaptation of certain quantities in determining the input model parameters are actually corroborated by the best-fit model. In Section 3.1.1, we used the empirical formulae to estimate the amount of RL contributions to [O iii] λ4363, [O ii] λ7320/30, and [N ii] λ5755 lines in deriving ${T}_{{\rm{e}}}$. The best-fit model yields ${I}_{{\rm{R}}}$/I([O iii] λ4363) = $0.67 \%$, ${I}_{{\rm{R}}}$/I([O ii] λ7320/30) = $1.13 \%$, and ${I}_{{\rm{R}}}$/I([N ii] λ5755) = $0.31 \%$, which are consistent with the empirical determinations adopted ($0.73 \%$, $2.19 \%$, and $0.54 \%$, respectively).

As for the ICFs used in determining the elemental abundances, we can compare the adopted ICFs based on I.P. and the ICFs calculated by the best-fit Cloudy model based on the ionization fraction of each element in the volume average in Table 7. While the values turn out to be consistent in general, discrepancies are found in Cl from the uncertain Cl⁺ fraction and in Si from the largely uncertain (Si) and ICF(I.P.). According to the best-fit model, the fraction of Cl⁺ to Cl is 0.38 and of Si⁺ to Si is 0.668.

As mentioned in the previous section (Section 4.2.8), the best-fit model is achieved by forcing the region of constant temperature at 1420 K in the PDR. This constant-temperature region is established from 5806 to 61'', that is, the radial temperature drops precipitously from 2750 K at 58'' to 1420 K at 58.06'', but is maintained at 1420 K from 5806 to 61'' to reproduce the observed molecular (H₂, CO, and OH) line fluxes. In this region, the relative proportion of molecular gas is maintained. So is the relative proportion of atomic gas.

In reality (of the presumed shocked H₂ scenario), however, shocked molecular regions are highly localized, and hence the relative proportion of molecular gas would keep increasing radially while that of atomic gas would keep decreasing. Therefore, with the presence of this constant-temperature PDR, the amount of the atomic gas component is bound to be overestimated in the PDR, i.e., the [C ii] and [O i] line fluxes are overpredicted (by a factor of 3–9; Figure 11, Table 15).

While our Cloudy model extends as far out as ${R}_{\mathrm{out}}=61^{\prime\prime}$, the optical ISIS and far-IR Herschel/PACS observations do not detect these [C ii] and [O i] lines with a sufficient signal level this far out in the PDR (i.e., the detection limit is reached at $R\approx 55^{\prime\prime}$). If we stopped model calculation at ${R}_{\mathrm{IB}}$ of 55'', we would obtain reasonable predictions of atomic line fluxes: for instance, I([O i] 63 $\mu {\rm{m}})=25.07$ (33.18, observed), I([O i] 145 $\mu {\rm{m}})=2.21$ (2.90, observed), and I([C ii] 157 $\mu {\rm{m}})=8.25$ (15.9, observed). However, of course, we would not be able to fit molecular lines at all (e.g., I(H₂ 9.67 $\mu {\rm{m}})=8(-5)$ for the model versus 25.79 observed).

In the present work, we adopt the average [C ii] and [O i] line fluxes measured in the entire PACS IFU field of view (over both of the "center" and "rim" positions; Figure 2), and the model-predicted [C ii] and [O i] line fluxes are deemed overestimated as a result. However, we actually measure fluxes as high as I([O i] 63 μm) = 103, I([O i] 145 μm) = 8.69, and I([C ii] 157 μm) = 27.24 in individual PACS spaxels over the barrel wall. Because there are no more data available to fit the model, especially the atomic component of the PDR, we have to leave these remaining discrepancies as issues to be resolved in the future when we obtain more sensitive data of the PDR and beyond. We will discuss the molecular component in detail later in Section 4.3.4.

To explore the possible O-rich nature of NGC 6781 (Section 3.1.4), we also construct the other "best-fit" model with amorphous silicate grains, adopting spherical grains of 0.05–0.50 μm radius (Figure 13). Overall, the best-fit model with amorphous carbon grains fit the observed continuum much better than the best-fit model with amorphous silicates. To fit the observed dust continuum with amorphous silicate grains, we have to reduce the amount of small grains in order not to produce any recognizable 10 μm silicate feature while achieving reasonable continuum fluxes in the far-IR. It is almost impossible to fit the dust continuum both in the mid-IR (10–40 μm) and in the far-IR ($\gt 70\,\mu {\rm{m}}$) simultaneously with amorphous silicate grains because amorphous silicates emit continuum only weakly beyond 70 μm. Therefore, we conclude that NGC 6781 was more likely C-rich in terms of the circumstellar dust composition.

Figure 12 shows how the best-fit model compares with the adopted post-AGB evolutionary tracks of Vassiliadis & Wood (1994). In the same plot, the best-fit model of NGC 6720 by van Hoof et al. (2010) is also displayed to confirm the similarity between the two in terms of the evolutionary status. A comparison between the evolutionary tracks implies that the progenitor of both NGC 6781 and NGC 6720 is a $\sim 2.5\,{M}_{\odot }$ star of Z = 0.02 and that the post-AGB age (i.e., the time since the cessation of AGB mass loss) is ∼9400 yr for NGC 6781.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

In addition, we plot in Figure 12 the evolutionary tracks of Miller Bertolami (2016, orange tracks of 2.0 and 3.0 M_☉ stars). These newer tracks are computed to address the shorter-than-expected timescales for Galactic bulge PNe. Their models with Z = 0.01 would take ∼3000, ∼2700, and ∼8000 yr to reach ${T}_{\mathrm{eff}}={\rm{120,870}}$ K for 2.0, 2.5, and 3.0 M_☉ stars, respectively, while models with Z = 0.02 would take ∼2600 to ∼12,000 yr to reach the same temperature for the 2.0 and 3.0 M_☉ models of Z = 0.02 for the same mass stars, respectively (no model track is given for 2.5 M_☉). Thus, the post-AGB age of a $2.5\,{M}_{\odot }$ progenitor with Z = 0.02 would be ∼3000 yr.

Following the method suggested by O'Dell et al. (2007), the empirical dynamical age of a PN can be approximated simply by

$$\begin{eqnarray}&&{t}_{\mathrm{dyn}}\simeq \displaystyle R/\tfrac{{V}_{\exp }(\mathrm{today})+{V}_{\exp }(\mathrm{AGB})}{2},\end{eqnarray} \tag{ 8 }$$

where ${V}_{\exp }(\mathrm{today})$ is the present-day shell expansion velocity and ${V}_{\exp }(\mathrm{AGB})$ is the shell expansion velocity at the beginning of the AGB phase. In this formulation, the shell expansion velocity is taken to be the rough "average" between the AGB wind velocity and the fast wind velocity. Assuming ${V}_{\exp }(\mathrm{AGB})=16\mbox{--}22$ km s⁻¹ (corresponding to the observed expansion velocity of the cold CO gas; Bachiller et al. 1993; Bergstedt 2015), ${V}_{\exp }(\mathrm{today})=12$ km s⁻¹ (from the [N ii] line; Arias & Rosado 2002), and ${R}_{\mathrm{IB}}=55^{\prime\prime}$ (the ionization front radius), the dynamical age would be roughly 7100–8600 yr.

Gesicki et al. (2016) suggested ${t}_{\mathrm{dyn}}\simeq (5/7)\times ({R}_{\mathrm{IB}}/{V}_{\exp })$ based on hydrodynamical model calculations. Adopting ${V}_{\exp }=12$ km s⁻¹ as above, the hydrodynamical age would be 7140 yr. Thus, the theoretical post-AGB age inferred from the Cloudy best-fit model and the evolutionary tracks by Vassiliadis & Wood (1994) is comparable to these (hydro)dynamical age estimates. Meanwhile, the much shorter post-AGB evolutionary time suggested by the evolutionary tracks by Miller Bertolami (2016) is more problematic to reconcile because the observed PN size would not be consistent with the observed expansion velocity, provided that the best-fit distance is 0.46 kpc (Section 3.4.1)

Here, we look into the molecular component of the best-fit model, especially into the PDR. We begin by comparing the model-predicted and empirically derived molecular column densities of H₂, CO, and OH⁺ (Table 8). The model-predicted results are derived by taking into account all of the gas components (i.e., molecular, atomic, and ionized), self-consistently allowing molecular formation processes (e.g., formation on dust grain surfaces and in the gas phase, and so on).

Table 8.  Comparison between the Best-fit Model-predicted and Empirically Derived Molecular Column Densities

Molecule	${\mathrm{log}}_{10}N$(Model)	${\mathrm{log}}_{10}N$(Obs)	Obs. References
	(cm−2)	(cm−2)
H2	18.18	18.36 ± 0.09	This work
CO	15.13	14.70–15.08	HerPlaNS1
OH+	13.00	10.54	Aleman et al. (2014)

Download table as:  ASCIITypeset image

As discussed above (Section 4.2.8), we introduced the warm-temperature component in the PDR as a necessary extra heating source to reproduce the observed H₂, CO, and OH lines. However, the achieved general agreement between the model and empirical column densities (Table 8) and line intensities (Figure 11; Table 15) permit qualitative characterization of the PDR in NGC 6781.

The best-fit floor temperature of 1420 K is consistent with the empirical estimates of $T({{\rm{H}}}_{2})=1279\pm 109$ K and 1161 ± 72 K by the single- and two-temperature excitation diagram fitting, respectively (Section 3.2.3). This suggests that H₂ is most likely in LTE and its kinetic temperature is about 1420 K. With this kinetic temperature, CO and OH lines are fit reasonably well. If we are to fit just the high-J CO lines, the best-fit floor temperature for CO would be 680 K. Either way (fitting with or without H₂), the (kinetic) temperature of CO gas would still be very much higher than an excitation temperature of ∼60 K (HerPlaNS1). This discrepancy can be mitigated if CO is assumed to be in non-LTE. Given the difference in the number density between H₂ and CO, CO could yet be being thermalized while H₂ already is.

Thus, we examine the excitation temperature of each CO line using the 1D non-LTE radiative transfer code RADEX (van der Tak et al. 2007). In RADEX calculations, we adopt the kinetic temperature of 1420 K, a constant n(H) = 10⁴ cm⁻³, and log₁₀ N(CO) = 15.13 cm⁻² as in the Cloudy model. The RADEX results (Table 9) suggest that the excitation temperature of high-J CO lines is 70–80 K on average, supporting the non-LTE condition for CO. We therefore conclude that the best-fit Cloudy model properly accounts for the presence of the warm component.

Table 9.  RADEX Non-LTE CO Model Results, and Comparison with the Observed Line Intensities

J	(μm)	${T}_{\mathrm{ex}}$ (RADEX)	Intensity (RADEX)	Intensity (Obs)
		(K)	(erg s−1 cm−2 sr−1)	(erg s−1 cm−2 sr−1)
4–3	650.3	209	3.51(−7)	3.73(−7) ± 6.35(−8)
5–4	520.2	85	7.71(−7)	7.67(−7) ± 2.03(−8)
6–5	433.6	70	1.19(−6)	1.17(−6) ± 1.51(−7)
7–6	371.7	70	1.47(−6)	1.99(−6) ± 2.49(−7)
8–7	325.2	74	1.58(−6)	9.71(−7) ± 1.38(−8)
9–8	289.1	82	1.55(−6)	1.08(−6) ± 2.91(−7)

Download table as:  ASCIITypeset image

The best-fit model predicts the amount of molecular gas in the PDR to be $4.15(-3)$ M_☉, which accounts only for the warm H₂ component (i.e., there are no other "cold" molecular components in the best-fit model). Meanwhile, this model prediction is actually consistent with the empirical estimate of $2.5(-3)$ M_☉ for the warm component (Section 3.2.4). However, the presence of the cold molecular component is very much expected based on the excitation diagram analysis (Section 3.2.3), as well as the non-LTE analysis we just saw above. In reality, there is probably a temperature gradient in the PDR along the polar direction, which empirically manifests itself as the multitemperature fit of the excitation diagram analysis and the non-LTE nature of the CO distribution.

Now, given that the best-fit model already properly accounts for the amount of ionized and neutral (atomic + molecular) gas, the cold molecular component that should exist must have been treated as part of the atomic gas component, as mentioned earlier (Section 4.2.8). Here, by adopting the ratio of the empirically determined cold H₂ mass to warm H₂ mass ($24.8=6.2(-2)\,{M}_{\odot }/2.5(-3)\,{M}_{\odot };$ Section 3.2.4), we can infer the amount of the cold molecular component to be expected in the best-fit model, $1.12(-1)\,{M}_{\odot }$ $(=4.15(-3)\,{M}_{\odot }\times 24.8)$. From this, we conclude that the modified best-fit model estimates of the mass of the cold molecular, warm molecular, and atomic gas components are $1.12(-1)\,{M}_{\odot }$, $4.15(-3)\,{M}_{\odot }$, and $1.99(-1)\,{M}_{\odot }$ $(=3.11(-1)\,{M}_{\odot }-1.12(-1)\,{M}_{\odot }$), respectively (see also Table 10).

Table 10.  Comparison of Each of the Mass Components between This Work and the HerPlaNS1 Results

Parameters	This Work	HerPlaNS1	HerPlaNS1
			(Scaled)
D (kpc)	0.46	0.95	0.46
Filling factor	1.0	0.5	1.0
Total gas (M☉)	0.41	0.86	0.40
Ionized gas (M☉)	0.09	0.54	0.25
Atomic gas (M☉)	0.20	0.12	0.05
Total molecular gas (M☉)	0.11	0.20	0.09
Warm molecular gas (M☉)	$4.15(-3)$	⋯	⋯
Cold molecular gas (M☉)a	$1.12(-1)$	⋯	⋯
Total dust mass (M☉)	$1.53(-3)$	⋯	⋯
Dust mass beyond ${R}_{\mathrm{IB}}$ (M☉)	$1.06(-3)$	$4.0(-3)$ b	$9.4(-4)$ b
GDR	268	335 (median)	335 (median)

Notes.

^aThe model-predicted cold molecular mass was scaled from the model-predicted warm molecular mass in this work. In HerPlaNS1, the molecular component was not subdivided by temperature. ^bThe empirical dust mass estimate was for the cold dust of 20–40 K.

Download table as:  ASCIITypeset image

We end the discussion on the molecular component in NGC 6781 by pointing out two lesser issues to be resolved that are beyond the scope of the present work. One is obviously the presence of the extra heating source. We incorporated the warm-temperature component in the model PDR assuming that shock interactions between the slower AGB wind and faster PN wind would provide sufficient extra heating to the PDR at the required level. Nonetheless, this extra heating source should be identified and self-consistently incorporated in the future. The other issue is the discrepancy in the OH⁺ column densities. This may well be due to a relatively more uncertain chemical network around OH⁺ and/or outdated reaction parameters in the astrochemistry network installed in Cloudy. However, the cause of the OH⁺ column density discrepancy is also unclear at this moment.

It is of interest to compare the amount of mass ejected during the AGB phase that is empirically accounted for with the adopted panchromatic data set (observational detection + model fitting via the present analyses) to our previous estimates based on an incomplete data set and to a theoretical prediction. As summarized in Table 10, the total gas mass empirically accounted for in NGC 6781 was $0.41\,{M}_{\odot }$, comprising $0.09\,{M}_{\odot }$ of ionized gas, $0.20\,{M}_{\odot }$ of atomic gas, and $0.11\,{M}_{\odot }$ of molecular gas. These values are based on the adopted volume filling factor f of unity (Section 3.4.1).

Previously, using almost exclusively far-IR line data and under the assumption of $D=0.95\,\mathrm{kpc}$, the total gas mass was estimated to be 0.86 M_☉, which consisted of 0.54 M_☉ of ionized gas (only ${{\rm{H}}}^{+}$, He⁺, and He²⁺), 0.12 M_☉ of atomic gas, and 0.20 M_☉ of molecular gas (only H₂ based on N(H₂) calculated from the excitation diagram), while adopting f = 0.5 (HerPlaNS1). With the updated distance of $D=0.46\,\mathrm{kpc}$ and f = 1, these previous estimates correspond to the total gas mass of 0.40 M_☉. While the total mass turns out to be consistent with the present result, the relative proportion of the individual gas components in the previous result is very different. This is, of course, because of the fact that we have to scale the relative proportion to fill gaps of the absence of sufficiently constraining observational data.

According to Karakas & Lattanzio (2007) and Karakas (2010), a 2.5 M_☉ initial-mass star with Z = 0.02 would experience 25 AGB thermal pulse (TP) episodes while ejecting the total mass of ∼1.25 M_☉. However, the predicted amount of the mass-loss ejecta would remain small ($\lt 0.01$ M_☉) until the 22nd TP episode. Over the last three TP episodes, the amount of the ejecta would increase precipitously, reaching $\simeq 0.70$ M_☉ during the last TP episode. Hence, our best-fit model accounts for roughly 60% of the amount of mass theoretically predicted to have been ejected during the last TP episode.

Meanwhile, the total gas mass within the ionization bound, ${R}_{\mathrm{IB}}=55^{\prime\prime}$, is $0.12\,{M}_{\odot }$ (consisting of $0.09\,{M}_{\odot }$ and $0.03\,{M}_{\odot }$ ionized and atomic gas, respectively), accounting for about 23% of the total gas mass. This proportion is consistent with a previous theoretical prediction made by Villaver et al. (2002), in which the evolution of the ejecta was modeled based on the stellar evolution tracks by Vassiliadis & Wood (1993). They concluded that the bright ionized shell would contain about 0.5 M_☉ of gas for a 2.5 M_☉ initial mass (their Figure 25), which roughly translates to 25% of the total ejecta mass.

Comparisons among these quantities indicate that the bulk of the nebular mass is found to be in the PDR of the nebula beyond ${R}_{\mathrm{IB}}$ in the form of neutral (atomic/molecular) gas. This finding is quite intriguing given the fact that PNe are generally known as the hallmark of the presence of ionized gas as ${{\rm{H}}}^{+}$ regions. The present work also demonstrates that PNe would provide a unique window of opportunities to investigate the mass-loss history of the progenitor star, because PNe should allow (1) access to a significant fraction of the AGB mass-loss ejecta when observed with sufficiently sensitive instruments (as opposed to AGB stars themselves) and (2) spatially resolved investigations more into the past (i.e., regions of larger radii) due to much larger energy input by the central star to illuminate the PDR of the nebula (as opposed to proto-PNe).

The best-fit model yields the dust mass (${m}_{\mathrm{dust}}$) of $1.53(-3)$ M_☉, while the empirically determined value obtained by fitting far-IR broadband images (HerPlaNS1), scaled to the present distance estimate of $D=0.46\,\mathrm{kpc}$, is $9.4(-4)$ M_☉. In both estimates, dust grain properties are the same (i.e., AC grains). This discrepancy is expected because the previous empirical estimate considered only the cold dust component detected in the far-IR (∼20–40 K; HerPlaNS1), missing the higher-temperature component emitting mainly in the shorter wavelength (e.g., mid-IR). The present best-fit model includes the entire (warm + cold) dust component (∼22–299 K).

To assess the consistency between the best-fit model and the empirical measurements, we estimate the mass of the cold/far-IR dust component in the best-fit model. Similar to the discussion in the previous section, we consider the cold dust component existing in the PDR beyond the IB, over which the model-predicted dust temperature would be 23–38 K. In the best-fit model, the dust mass beyond ${R}_{\mathrm{IB}}$ is $1.06(-3)$ M_☉, which is consistent with the empirical cold dust mass of 9.4(−4) M_☉.

The circumstellar dust mass is typically estimated via SED fitting of the thermal dust excess in the near- and mid-IR wavelengths. However, the present study reveals that there is a larger amount of cold dust (of $1.06(-3)$ M_☉) than warm dust (of $4.61(-4)$ M_☉) around NGC 6781. This finding suggests that the far-IR/cold dust component could take up a significant portion of the circumstellar dust in PNe (∼69% for the case of NGC 6781), and hence far-IR fluxes must always be incorporated in studying PNe, especially when considering the energetics in the whole volume of the nebula (especially the PDR and beyond).

In Cloudy model calculations, the presence of dust is scaled with the hydrogen density profile by the gas-to-dust mass ratio (GDR). The dust radiative transfer is done at each radial bin taking into account all the radiation available locally for dust heating (i.e., radiation from the ambient gas as well as from the CSPN). However, there is no mechanism to produce/destroy dust grains in the code. The best-fit model yields the "mean" GDR of 268 over the entire volume. The derived GDR is comparable to the average GDR of 386 ± 90 among 18 C-rich evolved stars (Knapp 1985) based on the direct comparison between the gas component (via CO $J=1-0$ observations in the radio, i.e., the cold gas component) and the dust component (via SED fitting of IR excess in the N band, i.e., the warm dust component). From our discussion in the previous section, it is likely that the Knapp work may have missed the cold dust component and hence their GDR may have been overestimated.

In our previous empirical estimate (HerPlaNS1), the GDR distribution in NGC 6781 shows a 10-fold decrease of the GDR from around 500 near the inner radius of the barrel wall to around 50 beyond the IB into the PDR with a median of 335. Caution needs to be exercised to compare these numbers because the empirical GDR distribution is susceptible to the projection effect (i.e., the gas and dust components being ratio-ed may not be present at the same location along the line of sight). Nevertheless, the median value is certainly consistent with the modeling results.

Gesicki et al. (2016) reported that 3D photoionization models could reproduce the observed emission line fluxes with ionized gas mass that is several times less than 1D models may suggest. This is because in 3D models there is usually a greater amount of "surfaces" at which ionization can happen. In 1D models, radiation would always have to be attenuated before penetrating into the next/outer radial layer of the nebula. However, in 3D models, attenuation may not even occur along some lines of sight (e.g., along the polar direction vs. equatorial directions in the case of a bipolar nebula), providing means to ionize the outer parts of the nebula to a greater extent. Indeed, we already saw some indication of the 3D effects especially in the PDR based on the multitemperature fit of the excitation diagram analysis and the non-LTE nature of the CO distribution, suggesting a temperature gradient along the polar direction of the nebula.

While 3D photoionization codes are available, we adopt the 1D Cloudy code because at this point no 3D photoionization codes would satisfactorily incorporate lower-temperature components (i.e., the dusty PDR) to be fit with the broad array of the adopted constraining observational data. For the case of NGC 6781 in particular, this 1D/3D issue implies that there could be a distribution of ionized gas extending along the polar directions (i.e., the regions of the polar caps and beyond), which would alter the overall proportion of the ionized gas in terms of the total mass of the nebula. However, these 3D effects on the ionized gas mass are considered to be minor in the present work. This is because model parameters that are critical in determining line fluxes and hence masses, such as the hydrogen density profile ${n}_{{\rm{H}}}$(R), D, L_*, ${T}_{\mathrm{eff}}$, nebular elemental abundances, and spatial distributions of various gas/dust components, were fixed to empirically derived values based on the spatially resolved data and not treated as free parameters, which is often the case in typical 1D models based on spatially unresolved data.