Physical properties of the fullerene C₆₀-containing planetary nebula   SaSt2-3★

ABSTRACT

We perform a detailed analysis of the fullerene C₆₀-containing planetary nebula (PN)  SaSt2-3 to investigate the physical properties of the central star (B0-1II) and nebula based on our own Subaru/HDS spectra and multiwavelength archival data. By assessing the stellar absorption, we derive the effective temperature, surface gravity, and photospheric abundances. For the first time, we report time variability of the central star’s radial velocity, strongly indicating a binary central star. Comparison between the derived elemental abundances and those predicted values by asymptotic giant branch (AGB) star nucleosynthesis models indicates that the progenitor is a star with initial mass of ∼1.25 M_|$\odot$| and metallicity Z = 0.001/α-element/Cl-rich ([α,Cl/Fe] ∼ +0.3–0.4). We determine the distance (11.33 kpc) to be consistent with the post-AGB evolution of 1.25 M_|$\odot$| initial mass stars with Z = 0.001. Using the photoionization model, we fully reproduce the derived quantities by adopting a cylindrically shaped nebula. We derive the mass fraction of the C-atoms present in atomic gas, graphite grain, and C₆₀. The highest mass fraction of C₆₀ (∼0.19  per cent) indicates that  SaSt2-3 is the C₆₀-richest PN amongst Galactic PNe. From comparison of stellar/nebular properties with other C₆₀ PNe, we conclude that the C₆₀ formation depends on the central star’s properties and its surrounding environment (e.g. binary disc), rather than the amount of C-atoms produced during the AGB phase.

1 INTRODUCTION

Mid-infrared (mid-IR) spectroscopic observations made by the Spitzer/Infrared Spectrograph (IRS, Houck et al. 2004) have recently detected fullerene C₆₀ and C₇₀ in a variety of space environments such as R Coronae Borealis stars (García-Hernández, Kameswara Rao & Lambert 2011a), reflection nebulae (Sellgren et al. 2010), young stellar objects (Roberts, Smith & Sarre 2012), post-asymptotic giant branch (AGB) stars (Gielen et al. 2011a,b), proto-planetary nebula (PNe; Zhang & Kwok 2011), and PN (Cami et al. 2010; García-Hernández et al. 2010, 2011b, 2012; Otsuka et al. 2013, 2014, 2016). At the moment, PNe represent the largest fraction of fullerene detection; since the first detection of the mid-IR C₆₀ and C₇₀ bands in the C-rich PN  Tc1 by Cami et al. (2010), 24 fullerene-containing PNe have been identified in the Milky Way and the Large and Small Magellanic Clouds (LMC and SMC, respectively).

In general, C₆₀ PNe show very similar IR dust features and stellar/nebular properties; their mid-IR spectra display broad 6−9, 11, and 30 |$\mu$|m features in addition to C₆₀ bands at 7.0, 8.5, 17.4, and 18.9 |$\mu$|m, and they have cool central stars and low-excitation nebulae, indicating that their age after the AGB phase is very young (e.g. Otsuka et al. 2014). The excitation mechanisms (e.g. Bernard-Salas et al. 2012) and the formation paths (e.g. Duley & Hu 2012; Berné, Montillaud & Joblin 2015) are not well understood and are still a subject of debate. However, it remains unclear why these objects exhibit the C₆₀ features – is the span of time during which spectral features of C₆₀ are present a short-lived phase that all C-rich PNe go through, or are C₆₀ PNe distinct objects in terms of their stellar/nebular properties and/or evolution? This is directly linked to the question of how C₆₀ forms in evolved star environments. We would like to answer this fundamental question by investigating the physical properties of C₆₀ PNe and comparing them with non-C₆₀ PNe.

Amongst C₆₀ PNe,  SaSt2-3 (PN G232.0+05.7, Acker et al. 1992) first identified by Sanduleak & Stephenson (1972) is a particularly interesting object to that we should pay more attention. Otsuka et al. (2014) discovered C₆₀ bands in this PN for the first time. Surprisingly, the mid-IR C₆₀ band strengths in  SaSt2-3 and  Tc1 are the strongest amongst all the fullerene-containing objects. This strongly indicates that the fullerene formation in these two PNe is particularly efficient.  Tc1 has been extensively studied since the discovery of C₆₀. However,  SaSt2-3 is not entirely understood due to the lack of available data for the central star and nebula and its uncertain distance (D). The uncertain D towards  SaSt2-3 has led to different estimates of the central star luminosity (L_*) and effective temperature (T_eff); accordingly, this has led to inconsistencies in understanding the evolutionary status of this PN (Gesicki & Zijlstra 2007; Otsuka et al. 2014). What we know from the prior studies is that this object has low-metallicity (ε(S)  = 5.48,¹ Pereira & Miranda 2007) and is (possibly) a Type IV PN (i.e., halo population, Pereira & Miranda 2007).

If we obtain the UV to optical wavelength spectra of the central star as well as the nebula, we can resolve issues raised and verify conclusions from previous studies of  SaSt2-3; by so doing, we can hope to gain insights into the C₆₀ formation. Fortunately, the UV-optical photometry data from the AAVSO Photometric All Sky Survey (APASS, Henden et al. 2016) can rigorously constrain L_*, and Frew, Parker & Bojičić (2016) improved its distance estimate (D = 14.86 ± 4.26 kpc). Therefore, we perform a comprehensive analysis on our own high-dispersion spectra of  SaSt2-3 taken using the 8.2 m Subaru telescope/high-dispersion spectrograph (HDS, Noguchi et al. 2002) and archived multiwavelength data.

We organize the next sections as follows. In Section 2, we describe our HDS spectroscopy and near-IR imaging using the NOAO WIYN 3.5 m/WIYN High-Resolution Infrared Camera (WHIRC, Meixner et al. 2010) and the reduction of this data. The Spitzer/IRS observation and its data reduction are described in Otsuka et al. (2014). In Fig. 1, we plot all the data used in the present work. In Section 3, we perform plasma diagnostics and derive ionic/elemental abundances. In Section 4, we derive photospheric elemental abundances, T_eff, and surface gravity |$\log \, g$| by fitting the stellar absorption using the theoretical stellar atmosphere code tlusty (Hubeny 1988). In Section 5, we compare the derived nebular and stellar elemental abundances with those values predicted by AGB nucleosynthesis models in order to infer the initial mass of the progenitor star. In Section 6, we build the spectral energy distribution (SED) model using the photoionization code Cloudy (v.13.05, Ferland et al. 2013) to be consistent with all the derived quantities based on our determined D. In Section 7, we discuss the origin and evolution of  SaSt2-3 and the C₆₀ formation in PNe by comparison of the derived nebular/stellar properties with other non-C₆₀ and C₆₀-containing PNe. Finally, we summarize the present work.

2 DATA SET AND REDUCTION

2.1 Subaru/HDS observation

We secured high-dispersion Echelle spectra using the HDS located at one of the Nasmyth loci of the 8.2 m Subaru Telescope at the top of Mauna Kea in Hawai’i. We summarize our observations in Table 1. We selected the 2 × 2 on-chip binning pattern. We set the slit-width to be 1.2″. We used the blue cross disperser for the 3640−5390 Å observation and the red one for the 4740−7490 Å and 7190−9960 Å observations, respectively. We utilized the atmospheric dispersion corrector (ADC) during the observations. In all the observations, we observed the standard star Hiltner 600 for correcting Echelle blaze functions and flux density simultaneously. In the 7190−9960 Å observation, we observed the telluric standard stars HD 61017 (B9III, m_V = 6.68) and HD 62217 (B9V, m_V = 8.26) at similar airmass.

We reduced the data using iraf² in a standard manner, including over-scan subtraction, scattered light subtraction between Echelle orders, and telluric absorption removal. We adopted the atmospheric extinction correction function measured by Buton et al. (2013) at Mauna Kea. We measured the actual spectral resolution (⁠|$R = \lambda /\delta {\lambda } =$| 32 300−33 500, see Table 1) using >300 Th–Ar comparison lines. The signal-to-noise ratio (S/N) for continuum reaches ∼40 at ∼3640 Å and ∼11 at ∼9950 Å. In ∼3700−4800 Å (this range is important in stellar absorption fittings), S/N is >70. We scaled both spectra to the average flux density in the overlapping regions (4740−5390 Å and 7190−7490 Å), and we connected these scaled spectra into a single 3640−9960 Å spectrum. The resultant spectrum is presented in Fig. 1.

2.2 Flux measurements and interstellar extinction correction

We determine c(Hβ) values by comparing the observed Balmer and Paschen line ratios to Hβ with the theoretical ratios of Storey & Hummer (1995) for the case with an electron temperature T_e  = 10⁴ K and an electron density n_e  = 2000 cm⁻³ under the Case B assumption. We calculate this n_e using the [O ii] F(3726 Å)/F(3729 Å) and the [Cl iii] F(5517 Å)/F(5537 Å) ratios. For the 3640−5390 Å spectrum, we obtain c(Hβ)  = 0.20 ± 0.01, which is the intensity-weight average amongst the Hγ, Hδ, Hε, and Hη to the Hβ ratios. For the 4740−7490 Å spectrum, we obtain c(Hβ)  = 0.40 ± 0.01 from the Hα/Hβ ratio. For the 7190−9960 Å spectrum, we determine c(Hβ)  = 0.24 ± 0.03 from the Paschen H i 9014 Å (P10) to the Hβ ratio. For all HDS spectra, we adopt the average c(Hβ)  = 0.28 ± 0.11 amongst three HDS observations.

Tylenda et al. (1992) reported c(Hβ)  = 1.11 (observation date is unknown). We derive the average c(Hβ)  = 1.11 ± 0.26 using the ratio of F(Hα), F(Hγ), and F(Hδ) to F(Hβ) measured from their spectra.³ Based on the F(Hα) and F(Hβ) reported by Dopita & Hua (1997), we obtain c(Hβ)  = 0.43 ± 0.04 (obs date: 1997 March). Pereira & Miranda (2007) reported E(B − V) = 0.11 ± 0.02, which corresponds to c(Hβ)  = 0.13−0.19 (obs date: 2005 Feb.). Using the line flux table of Pereira & Miranda (2007), we obtain the average c(Hβ)  = 0.41 ± 0.23 calculated from F(Hα), F(Hγ), and F(P10) to F(Hβ). Using the archived ESO Faint Object Spectrograph and Camera (EFOSC) spectrum taken on 2000 April,⁴ we obtain a c(Hβ)  = 0.68 ± 0.11 measured from the F(Hα)/F(Hβ) ratio. A time variation of c(Hβ) seen between 1992 and 2016 might be due to the effect of stellar H i absorption to corresponding nebular H i and also orbital motion of the binary central star (Section 4.2).

We scale the Spitzer/IRS spectrum to match the Wide-field Infrared Survey Explorer (WISE) W3/W4 band flux densities of Cutri (2013) and the L18W AKARI/IRC mid-infrared all-sky survey of Ishihara et al. (2010) (see Section 2.4). For this scaled Spitzer/IRS spectrum, we do not correct interstellar extinction because the interstellar extinction is negligibly small in the mid-IR wavelength. It is common practice in nebular analyses to scale all line intensities in such a way that Hβ has a line flux of 100. To achieve this, we first normalize the line fluxes with respect to the complex of the H i 7.46 |$\mu$|m (n = 5−6, n is the quantum number) and 7.50 |$\mu$|m (n = 6−8) lines. F(7.48/7.50 |$\mu$|m) is (4.87 ± 0.30) × 10⁻¹⁵ erg s⁻¹ cm⁻² (A(−B) means A × 10^−B hereafter). According to Storey & Hummer (1995) for the Case B assumption with T_e  = 10⁴ K and n_e  = 2000 cm⁻³, the ratio of H iI(7.48/7.50 |$\mu$|m)/I(Hβ)  = 3.102/100. Finally, we multiply all the normalized line fluxes by 3.102 to express them relative to Hβ with I(Hβ)  = 100.

In Appendix Table A1, we list the identified emission lines in the Subaru/HDS and Spitzer/IRS spectra. The first column is the laboratory wavelength in air. Here, I(Hβ) is 100. The last column δ I(λ) corresponds to 1-σ.

2.3 NOAO/WHIRC near-IR imaging observation

We took the high-resolution images using NOAO WIYN 3.5 m/WHIRC. We summarize the observation log in Table 1. We took the two narrowband images using the Brγ (λ_c = 2.162 |$\mu$|m, effective bandwidth (W_eff)  = 0.210 |$\mu$|m), and Brγ45 (λ_c = 2.188 |$\mu$|m, W_eff = 0.245 |$\mu$|m) filters.⁵ We selected a 5 pts dithering pattern. We followed a standard manner for near-IR imaging data reductions using iraf, including background sky and dark current subtraction, bad pixel masking, flat-fielding, and distortion correction. Finally, we obtained a single averaged image for each band.

For the flux calibration, we utilized the SED of the standard star 2MASS07480394−1407155 based on its photometry between the 2 micron all sky survey (2MASS, Cutri et al. 2003) JHKs and WISE bands W1/W2. The SED of this star can be well fitted with a single blackbody temperature of 3510 K. Then, we derived the flux density in each band by taking each filter transmission curve into account. Next, we measured the respective count of the standard star in the Brγ and Brγ45 images. Thus, we obtained the conversion factor from the counts in ADU to flux density in erg s⁻¹ cm⁻² |$\mu$|m⁻¹. The measured flux density in each band is listed in Appendix Table A2 and plotted in Fig. 1 (green circles).

2.4 Photometry data

To support the present work, we collected the data taken from APASS, 2MASS, WISE, and AKARI/Infrared Camera (IRC). In Table A2, we list the observed and reddening corrected flux densities F_λ and I_λ, respectively. We obtain I_λ using equation (1) and the average c(Hβ)  = 0.28 ± 0.11 amongst three HDS observations (Section 2.2). Due to negligibly small reddening effect, we do not correct F_λ in the longer wavelength than WISE W1 band (3.35 |$\mu$|m).

3 NEBULAR LINE ANALYSIS

3.1 Systemic nebular radial velocity

We obtain the average heliocentric radial velocity of  +166.6 km s⁻¹ measured from the identified 128 nebular lines in the HDS spectrum (the standard deviation is 3.2 km s⁻¹ amongst all these lines and that of each radial velocity is 0.52 km s⁻¹ in the average). The LSR radial velocity v_r(LSR) of  +149.1 km s⁻¹ is much faster than v_r(LSR) in other Galactic PNe towards l ∼220−240° and b ≲ ±10° (≲ +80 km s⁻¹; Quireza, Rocha-Pinto & Maciel 2007). v_r(LSR) of  +105 km s⁻¹ in the C₆₀ PN M1-9 (PN G214.0+04.3, Otsuka et al. 2014) is the closest to  SaSt2-3’s v_r(LSR) (Quireza et al. 2007). Peculiar velocity relative to Galactic rotation is calculated using v_r(LSR) and D in order to classify PNe into Type I-IV (i.e. thin/thick disc and halo; see e.g. Peimbert 1978). We discuss classification of  SaSt2-3 in Section 7.1. We do not find a time-variation of the radial velocity measured by the nebular lines. We report the radial velocity measurements from the stellar absorption in Section 4.2.

3.2 H β flux of the entire nebula

From the measured F(H i 7.48/7.50 |$\mu$|m) and the theoretical I(H i 7.48/7.50 |$\mu$|m)/I(Hβ) ratio  = 3.102/100 (see Section 2.2), we obtain I(Hβ) of the entire nebula to be (1.57 ± 0.28)(−12) erg s⁻¹ cm⁻². F(Hα) and F(H β) using the 5″ wide slit observation by Dopita & Hua (1997) and c(Hβ)  = 0.43 ± 0.04 (Section 2.2) yields I(Hβ)  = (1.97 ± 0.20)(−12) erg s⁻¹ cm⁻², which is consistent with ours. In the present work, we adopt our own calculated I(Hβ) for the entire nebula because our I(Hβ) is based on interstellar extinction-free and stellar H i absorption effect less mid-IR H i 7.48/7.50 |$\mu$|m.

3.3 Plasma diagnostics

We determine n_e and T_e using diagnostic line ratios listed in Table 2, with the resulting n_e – T_e diagnostic curves for the collisionally excited lines (CELs) shown in Fig. 2. The roughly vertical (blue) lines can be used to determine n_e; more horizontal (grey) lines T_e. Although the other diagnostic curves (red) yield both n_e and T_e, we use them as T_e indicators here. Since the emission of each ion originates from regions of different T_e and n_e, we need to determine both parameters for each ion independently so that we can determine accurate ionic abundances later on. This involves several steps. First, we note that T_e ∼9000 K from T_e([Ar iii]), T_e([S iii])_opt, and T_e([N ii]) curves. Next, we adopt T_e  = 9000 K to solve each equation of population at ≥ 5 multiple energy levels for each n_e sensitive ions; from this, we then calculate n_e from the corresponding diagnostic line ratios for [N i], [S ii], [O ii], [Cl iii], and [S iii] (see Table 2). Note that the precise T_e we assume here does not matter much, since these n_e diagnostic line ratios are fairly insensitive to T_e. With the n_e values established, we then determine T_e by adopting the derived n_e values corresponding to each ion. We adopt n_e([O ii]) for T_e([N ii]) derivation. Our derived values are in agreement with those by Pereira & Miranda (2007) who found n_e([S ii])  = 2100 ± 600 cm⁻³ and T_e([N ii])  = 9600 ± 930 K.

We compute T_e(He i) using singlet He i lines. To calculate T_e(PJ) from the Paschen continuum discontinuity by utilizing the equation (7) of Fang & Liu (2011), first we determine the He⁺ abundance of 1.09(−2) ± 2.28(−4) under the obtained T_e(He i). Eventually, we utilize T_e(He i) for both He⁺ and C²⁺ abundance calculations due to higher T_e(PJ) uncertainty.

3.4 Ionic abundance derivations

We calculate the CEL ionic abundances by solving an equation of population at multiple energy levels under the adopted T_e and n_e as listed in Table 3; T_e  = 9270 K is the average value amongst two T_e([S iii]) and T_e([Ar iii]), T_e  = 9790 K is the average value amongst two T_e([S iii]), T_e  = 9440 K is the average between T_e([N ii]) and T_e([O ii]), and n_e  = 1690 cm⁻³ is the average between n_e([Cl iii]) and n_e([S iii]). For the recombination line (RL) He⁺ and C²⁺, we adopt T_e(He i) and n_e  = 10⁴ cm⁻³. Our choice of the T_e–n_e pair of each ion depends on the potential (IP) of the targeting ion. Except for the CEL N⁺, O^{+, 2+}, and S⁺ which Pereira & Miranda (2007) already measured, the first measurements of all the ionic abundances are done by us. We summarize the resultant CEL and RL ionic abundances in Appendix Table A3. We calculate each ionic abundance using each line intensity. Then, we adopt the weight-average value as the representative ionic abundance as listed in the last line of each ion. We give 1-σ uncertainty of each ionic abundance, which accounts for the uncertainties of line fluxes (including c(Hβ) uncertainty), T_e, and n_e.

The He⁺ abundance of 9.72(−3) in SaSt2-3 is 10 times smaller than in evolved PNe (e.g. |$T_{\rm eff} \gtrsim 50\, 000$| K). For instance, in the C₆₀ PN  M1-20 (⁠|$T_{\rm eff} = 45\, 880$| K, Otsuka et al. 2014), Wang & Liu (2007) find He⁺ abundance of 9.50(−2). Moreover, the He⁺ abundance is also significantly lower than in other Galactic C₆₀ PNe with |$T_{\rm eff} \lesssim 40\, 000$| K where He⁺ abundances have been determined: 6.99(−2) in  IC418 (Hyung, Aller & Feibelman 1994), 6.57(−2) in  M1-6 (Otsuka in prep), 3.93(−2) in  M1-11 (Otsuka et al. 2013), 3.5(−2) in  M1-12 (Henry et al. 2010), and 6.0(−2) in  Tc1 (Pottasch, Surendiranath & Bernard-Salas 2011). Similar to the C₆₀ PN  Lin49 in the SMC (Otsuka et al. 2016), the low He⁺ abundance is due to the smaller number of ionizing photons for He⁺ ( ≥ 21 eV): using the spectra synthesized by tlusty (with L_* = 7000 L_|$\odot$|, |$\log \, g$|  = 3.11 cm s⁻², metallicity Z = 1/10 Z_|$\odot$|, see below), we estimate the number of photons with energy ≥ 21 eV to be 8.3(+45) s⁻¹ in a |$T_{\rm eff} = 28\, 100$|K star-like  SaSt2-3 (see Section 4) and 4.8(+46) s⁻¹ in |$T_{\rm eff} = 32\, 000$|K stars like  M1-11 and  M1-12. Thus, the majority of the He atoms in  SaSt2-3 are in the neutral state.

The higher multiplet C ii lines are generally reliable because these lines are less affected by resonance fluorescence. However, the higher C²⁺ abundances from the C ii 3918.98/20.69 Å (4s²S − 3p²P) and 7231.32/36.42 Å (3d²D − 3p²P) are likely due to the enhancement by resonance from the 635.25/636.99 Å (4s²S − 2p²P) and the 687 Å (3d²D − 2p²P), respectively. Thus, we exclude the C^{2 +} abundances from these C ii lines and C ii 6451.95 Å⁶ in the representative RL C²⁺ determination.

Our N⁺ and O^{+, 2+} are comparable with Pereira & Miranda (2007), who calculated N⁺ = 2.42(−5), O⁺ = 1.87(−4), O²⁺ = 1.22(−6), and S⁺ = 3.0(−7) (they note that their derived ionic abundances has ±30|${{\ \rm per\ cent}}$| uncertainty) under T_e([N ii])  = 9600 K and n_e([S ii])  = 2100 cm⁻³. The discrepancy between their and our S⁺ (5.89(−7)) is caused by T_e selection; if we adopt T_e  = 9600 K and n_e  = 2100 cm⁻³, we obtain S⁺ = 3.48(−7).

3.5 Elemental abundance derivations using the ICFs

To obtain the elemental abundances using the derived ionic abundances, we introduce the ionization correction factors (ICFs, see e.g. Delgado-Inglada, Morisset & Stasińska 2014, for details). Here, the number density ratio of the element X with respect to the hydrogen, n(X)/n(H) is equal to ICF(X) · ∑_{m = 1} n(X^m+)/n(H⁺). The ICFs have been empirically determined based on the fraction of observed ion number densities with similar ionization potentials to the target element, and have also been determined based on the fractions of the ions calculated by photoionization models. Since  SaSt2-3 is very low-excitation PN, the ICFs (He in particular) adopted for more highly excited PNe do not work well. Therefore, we need a special treatment for  SaSt2-3. Thus, in addition to the ICFs established by photoionization grid models of Delgado-Inglada et al. (2014), we refer to ICFs used in  Lin49 by Otsuka et al. (2016).

In  Lin49, the C²⁺/C ratio is similar to the Ar²⁺/Ar ratio. For  SaSt2-3, we adopt equation (A6) of Otsuka et al. (2016) for C and Ar. Both ICF(C) and ICF(Ar) are 2.7 · (Cl/Cl²⁺)  = 5.12 ± 1.41. For He derivation, we adopt two ICF(He) calculated using the equation (42) of Peimbert & Costero (1969, 11.57) and from the ratio of Ar/Ar^{2 +} (5.12). ICF(N)  = 1.04 ± 0.13 is from Delgado-Inglada et al. (2014). ICF(Fe)  = 1.31 ± 0.16 is from Delgado-Inglada & Rodríguez (2014) based on the observation results. ICF of the other elements is unity. We verify whether the adopted ICFs here are proper by comparing with Cloudy photoionization model (Section 6).

In the second and third columns of Table 4, we present the resultant elemental abundances with 1-σ uncertainty, except for ε(He), where we adopt its range. The two columns are the relative abundance to the solar value by Lodders (2010) and the ε(X) by Pereira & Miranda (2007). Our ε(N) and ε(O) are consistent with Pereira & Miranda (2007). As explained in Section 3.4, ε(S) discrepancy between theirs and ours is attributed to the S⁺ abundance. By the Cloudy model under D  = 6 kpc, Otsuka et al. (2014) derived ε(N/O/Ne/S/Ar)  = 7.49, 8.23, 7.68, 6.17, and 5.93 based on the optical spectrum of Pereira & Miranda (2007) and the Spitzer/IRS spectrum. Otsuka et al. (2014) estimated an expected CEL ε(C)  = 8.72 using a [C/H]−[C/Ar] relation established amongst 115 Galactic PNe. Delgado-Inglada & Rodríguez (2014) reported that the RL C²⁺ to the CEL C²⁺ ratio in  IC418 is 2.4. Applying this value to  SaSt2-3, we obtain an expected CEL ε(C)  = 8.83 ± 0.12, which is consistent with Otsuka et al. (2014). We attempt to obtain more plausible expected CEL ε(C) using the stellar ε(C) and ε(O) in Section 4.1. The [Ne/H] is comparable with the [O/H] because Ne together with O had been synthesized in the He-rich intershell during the AGB phase. The Ne enhancement would be due to the increase of ²²Ne.

The [S,Cl,Ar/H] abundances are low, and if these represent the stellar abundances, then  SaSt2-3 is the lowest metallicity object amongst the Galactic C₆₀ PNe, and we infer Z ∼0.1 Z_|$\odot$| from the average [S,Cl,Ar/H]. While most Ar is probably in the gas phase in this object, S could be incorporated into dust grains (e.g. MgS, suggested to be a candidate for the carrier of the broad 30 |$\mu$|m feature that is observed in the C₆₀ PNe). Fe is even more depleted, but it is unlikely that this represents the initial abundance given the other elemental abundances. Rather, a fraction of the Fe will be incorporated into dust grains. We discuss further the elemental abundances in Section 5.

4 STELLAR ABSORPTION ANALYSIS

4.1 Stellar parameter derivations

We perform stellar absorption analysis of the HDS spectra taken on 2013 Oct 6 and Dec 10 using the non-local thermodynamic equilibrium (non-LTE) stellar atmosphere modelling code tlusty. We detect strong Si iii,iv and He ii absorption lines. From our tlusty modelling, T_eff of the central star is 28 100 K (Table 5), which is cooler than |$T_{\rm eff} \ge 30\, 000$|K in the O-type stars. Thus, we classify the stellar spectrum of  SaSt2-3 into early B-type giant B0-1II rather than O-type. Thus, we use comprehensive grid of 1540 metal line-blanketed, non-LTE, plane-parallel, hydrostatic model atmospheres of B-type stars BSTAR2006⁷ by Lanz & Hubeny (2007).

We find the average nebular [Cl,S,Ar/H] of −0.89 ± 0.12 (Section 3.5). Assuming that the metallicity Z of the central star and the nebula is roughly the same as seen in the case of IC 418 by Morisset & Georgiev (2009), we adopt the Z = 0.1 Z_|$\odot$| model grid from BSTAR2006. All of the initial abundances in this model grid are set to ε(He) = 11.00 and [X/H]  = −1 except for He. Based on the Z = 0.1 Z_|$\odot$| model grid, we vary ε(X) to yield each equivalent width (EW) of element X to compare with each EW(X) measured from the observed HDS spectra. Throughout our tlusty synthesis analysis, we do not set [He,C,N,O,Si/H]  = −1 and we do not adopt the derived nebular He,C,N,O/H] as the stellar photospheric ones. Based on the measured EW of the identified 9 He i,ii, 4 C iii,iv, 2 N ii,iii, 13 O ii, and 5 Si iii,iv absorption, we derive the photospheric He/C/N/O/Si abundances, microturbulent velocity (v_t), rotational velocity (⁠|$v\, \sin (i)$|⁠; i is the angle between the rotation axis and the line of sight), T_eff, and |$\log\, g$| of the central star. These absorption lines are lesser affected by the nearby nebular lines and absorption lines of the other elements. As we report later, the central wavelength of the stellar absorption lines changes between observing dates whereas those of the nebular lines remain constant. Before analysis, we convert heliocentric wavelength frame of the HDS spectrum into rest frame using the radial velocity determined by the He ii4686 Å for the 2013 Oct 6 data (+154.8 ± 2.3 km s⁻¹) and the He ii5411 Å for the 2013 Dec 10 data (+182.3 ± 4.0 km s⁻¹).

First, we set the basic parameters characterizing the stellar atmosphere, i.e. Z, v_t, T_eff, and |$\log\, g$|⁠. By setting |$T_{\rm eff} = 28\, 000$| K and |$\log \, g = 3.10$| cm s⁻², we investigate ε(O) versus selected 8 O ii lines’ EW to determine v_t using synfit. For each absorption, we set instrumental line broadening determined by measuring Th-Ar line widths. Since v_t ≥ 10 km s⁻¹ gives minimization of the scatter in ε(O) versus EW, we adopt v_t = 10 km s⁻¹. As a reference, Morisset & Georgiev (2009) adopted v_t = 10 km s⁻¹ for  IC418.

We determine T_eff and |$\log\, g$| using the |$T_{ \mathrm{ eff}} - \log\, g$| curves generated by the model atmosphere with v_t = 10 km s⁻¹, Z  = 0.1 Z_|$\odot$|, and ε(He)  = 10.90. The |$T_{ \mathrm{ eff}} - \log\, g$| curves are generated by the following process; for a fixed T_eff, we vary |$\log \, g$| from 2.90 − 3.30 cm s⁻² in a constant 0.01 cm s⁻² step to find the best fit value for each observed He i,ii’s EW. We test the range of T_eff from  27 500 − 28 500 K (200 K step).

Based on the determined v_t, T_eff, and |$\log\, g$|⁠, we further constrain ε(He) and calculate ε(C/N/O/Si) abundances by comparing the observed and model predicted EWs of each line by synabund.⁸ We summarise the result in Table 5. In Appendix Table A4, we list the elemental abundances using each line. We adopt the average value as the representative abundance as listed in the last line of each element. The uncertainty of elemental abundances includes errors from the measured EWs, T_eff, |$\log\, g$|⁠, and the uncertainty when we adopt the model atmosphere with the [Z/H]  = −0.90 or −1.10 and when we assume the uncertainty of v_t of 2 km s⁻¹. We determine |$v\, \sin (i)$| by line-profile fittings of the selected He i,ii and H i in  4000−4700 Å using synfit⁹

In Fig. 3, we show the synthetic stellar spectrum generated using synspec.¹⁰ We identify the absorption lines (except for H i) with the model predicted EW ≥ 10 mÅ by the blue lines. The synthetic spectrum in 3700−4750 Å is presented in Appendix Fig. A1. Stellar Ne, S, Ar, Mg, Ca, and Ti (α-elements) and Ni, Fe, and Zn are not derived in optical HDS spectra of  SaSt2-3. These abundances are not small and also they are very important in characterizing the spectrum of the central star and its radiation hardness (in particular, X-ray to UV wavelength). We know that the central star radiation is suppressed by the metal line-blanket effect and also it is very important in subsequent Cloudy modelling. Thus, it is worth simulating these elements, too. We adopt the nebular Ne, S, Cl, and Ar abundances due to no detection of stellar absorption of these elements. We adopt ε(Fe) = 6.38 ([Fe/H]  = −1.1, see Section 3.5). Based on the discussion in Section 5, for the other elements up to Fe except for α-elements Mg, Ca, and Ti, we adopt the predicted values by the AGB nucleosynthesis model of initially 1.25 M_|$\odot$| and Z = 0.001 stars by Fishlock et al. (2014). For Mg, Ca, and Ti, we adopt ε(Mg) = 6.80, ε(Ca) = 5.43, and ε(Ti) = 4.05, respectively (i.e. [Mg,Ca,Ti/H]  = −0.7).

The stellar ε(He/N/O) is in agreement with the nebular ε(He/N/O) within their uncertainties. The stellar C/O ratio (2.93 ± 1.33) indicates that  SaSt2-3 is definitely a C-rich PN. Based on the consistency between the stellar and the nebular elemental abundances, we obtain an expected CEL ε(C)  = 8.58 ± 0.20 using the CEL ε(O) and the stellar C/O ratio.

4.2 Time variation of line profile and radial velocity; Evidence of a binary central star

Our important discovery is that the central wavelength of the stellar absorption lines varies from date to date whereas there is no wavelength shift of the nebular emission lines.

We compute the heliocentric radial velocities v_r of the central star via Fourier cross-correlation between the observed spectra and the synthetic tlusty spectrum using fxcor in iraf. FXCOR calculates the velocity shift between two different spectra in the selected wavelength regions.¹¹ Here, we select good S/N regions. In Table 6, we list v_r and v_r − v_sys, where v_sys is the systemic radial velocity measured from the 128 nebular emission lines ( +166.6 km s⁻¹, see Section 3.1). In Fig. 4, we show the singlet He i 5015/7281 Å absorption and the tlusty synthetic spectrum as the guide.

We interpret that the radial velocity time-variation is caused by orbital motion in a binary system. Méndez, Forte & López (1986) reported the photometric and radial velocity variations of the CSPN of C₆₀ PN  IC418. They measured the radial velocities using the stellar C iii 5695 Å and C iv 5801/11 Å. The systemic radial velocity was derived using the nebular [N ii] 5755 Å line. Later, Méndez (1989) concluded that the central star is not likely to be a binary because the orbital motion alone (if present) would not be enough to explain the observed variations. We note that the C iii 5695 Å and C iv 5801/11 Å lines are good indicators of the stellar activity (e.g. wind velocity) and these lines would be unlikely to give more accurate radial velocity of the central star. Thus, as far as we know, this would be the firm detection case of the binary central stars amongst all the C₆₀ PNe. Since we have only three periods of the binary motion, we do not determine any parameters of the binary central star yet.

We expected near-IR excess from the binary circumstellar disc from Otsuka et al. (2016) who detected near-IR excess in most of the SMC C₆₀ PNe and discussed possible links between near-IR excess, disc, and fullerene formation; since the ejected material from the central star can be stably harboured for a long time, even smaller molecules could aggregate into much larger molecules. However, in  SaSt2-3, we do not find near-IR excess in the observed SED (Fig. 1). No near-IR excess might mean a possibility of a nearly edge-on disc rather an inclined disc.

5 COMPARISON WITH AGB MODEL PREDICTIONS

In Table 7, we compile the derived abundances. The nebular CEL ε(C) is an expected value by our analysis (Section 4.1). As the comparisons, we list the AGB nucleosynthesis model predictions by Fishlock et al. (2014) for initially 1.25 M_|$\odot$| stars with Z = 0.001 and Karakas (2010) for initially 1.50 M_|$\odot$| stars with Z = 0.004. Note that Fishlock et al. (2014) and Karakas (2010) set the initial [X/H] to be −1.1 and −0.7, respectively. We calculate reduced chi-squared values (⁠|$\chi _{\nu }^{2}$|⁠, ν is degree of freedom) between the nebular ε(X) and the AGB model predicted values for each of 1.00 − 3.25 M_|$\odot$| star with Z = 0.001 (nine models in total). We use |$\chi _{\nu }^{2}$| as the guide to find out which AGB model’s predicted abundances is the closest to the derived abundances. The aim of this analysis is to infer the initial mass of the progenitor. We should note that these AGB grid models do not aim to explain the observed elemental abundances of  SaSt2-3. We exclude Fe in |$\chi _{\nu }^{2}$| evaluation. We adopt the nebular ε(X) values. For ε(He) and ε(C), we adopt 10.96 ± 0.17 (intermediate value, 9.10(−2) ± 3.52(−3)) and 8.58 ± 0.20 (an expected nebular CEL C value, 3.80(−4) ± 1.76(−4)), respectively.

Since the reduced-χ² for the 1.25 M_|$\odot$| model marks the minimum (⁠|$\chi _{7}^{2}$| = 14 (= 99/(8 – 1)) in 8 elements), this model is the closet to the derived ε(X). |$\chi _{4}^{2}$| is 17 (= 66/(5 – 1)) limited to ε(He/C/N/O/Ne). Next, we compare the AGB models for the same mass stars with Z = 0.004 because these models could account for the derived abundances except for S (no predictions for Cl and Ar, however). The model for 1.5 M_|$\odot$| initial mass stars with Z = 0.004 gives the closest fit to the observation (⁠|$\chi _{5}^{2}$| = 360 (=1800/(6 – 1))in 6 elements). Limited to ε(He/C/N/O/Ne), |$\chi _{4}^{2}$| is 7 (=28/(5 – 1)).

The B-type central star indicates that  SaSt2-3 is an extremely young PN and just finished the AGB phase. The presence of H absorption lines (Fig. 3) indicates that  SaSt2-3 did not experience very late thermal pulse evolution, so this PN is probably in the course of H-burning post-AGB evolution. According to the H-burning post-AGB evolution model of Vassiliadis & Wood (1994), stars with initially 1.5 M_|$\odot$| and Z = 0.004 would evolve into hot stars with the core mass (M_*) of 0.64 M_|$\odot$|. L_* of such stars is ∼7380 L_|$\odot$| when T_eff is ∼28 100 K in ∼1050 years after the AGB phase. Whereas, we infer that 1.25 M_|$\odot$| stars with Z = 0.001 would evolve into stars with M_* of 0.649 M_|$\odot$|; their L_* and T_eff is ∼7765 L_|$\odot$| and ∼28 100 K, respectively in ∼1770 yr after the AGB-phase based on the models of Fishlock et al. (2014) and Vassiliadis & Wood (1994). The main difference in the post-AGB evolution of 1.50 M_|$\odot$|/Z = 0.004 stars and 1.25 M_|$\odot$|/Z = 0.001 stars is evolutionary time-scale.

Through these discussions, we summarise as follows. The Z = 0.001 model gives the closest values to the derived elemental abundances, although there is systematically ∼0.3–0.4 dex discrepancy of ε(O,Cl,Ar). The Z = 0.004 model shows excellent fit to the derived nebular and stellar ε(He/C/N/O/Si). Considering initial settings of [X/H] in the models, we conclude that the progenitor of  SaSt2-3 would be a ∼1.25 M_|$\odot$| star with initially Z∼0.001 ([Fe/H] ∼ −1.1) and [α,Cl/Fe] ∼ +0.3–0.4. This is consistent or comparable with the Galaxy chemical evolution model of Kobayashi, Karakas & Umeda (2011); in the Galactic thick disc, the predicted [Si/Fe], [S/Fe], [Cl/Fe], and [Ar/Fe] are ∼+0.6, ∼+0.4, ∼−0.3, and ∼ +0.3 in [Fe/H] < −1, respectively.

6 PHOTOIONIZATION MODEL

In the previous sections, we characterized the central star and dusty nebula. In this section, we build the photoionization model using Cloudy and tlusty to be consistent with all the derived quantities, AGB nucleosynthesis model, and post-AGB evolution model. Below, we explain how to set each parameter in the model, and then we show the result.

6.1 Modelling approach

6.1.1 Distance

Since the distance D is an important parameter, we estimate it by our own method as explained below. We first extract the stellar spectrum from the observed HDS spectrum because the observed spectrum is the sum of the nebular emission lines and continuum and the central star’s continuum. For this purpose, we scale the HDS spectrum flux density to match the APASS Bg′Vr′i′ bands. Then, we subtract the theoretically calculated nebular continuum from the scaled HDS spectrum. We utilize the nebcont code in the Dispo package of starlink v.2015A¹² to generate the nebular continuum under adopting I(Hβ)  = 1.57(−12) erg s⁻¹ cm⁻² (Section 3.2), T_e  = 10⁴ K, n_e  = 2000 cm⁻³ (Table 2), and n(He⁺)/n(H⁺)  = 1.09(−2) (Table A3). In Fig. 5(a), we show the scaled HDS spectrum and the synthetic nebular continuum. Fig. 5(b) displays the tlusty synthetic spectrum of the central star (in the case of |$T_{\rm eff} = 28\, 100$| K) scaled to match the residual spectrum generated by subtracting the nebular continuum from the HDS spectrum.

Assuming that the progenitor is an initially 1.25 M_|$\odot$|/Z = 0.001 star and its luminosity is currently 7765 L_|$\odot$| (Section 5), we obtain D = 9.90 − 12.76 kpc. If we assume an initially 1.5 M_|$\odot$| progenitor star with Z = 0.004, D is 9.66 − 12.44 kpc.

In our Cloudy model, we adopt D = 11.33 kpc, which is the intermediate value of D when we assume that the central star evolved from a star with initially 1.25 M_|$\odot$| and Z = 0.001. Adopting our measured Galactocentric distance of 17.35 kpc, the predicted ε(O/Ne/Cl/S/Ar) from the Galaxy ε(O/Ne/Cl/S/Ar) gradient established amongst Galactic PN nebular abundances by Henry, Kwitter & Balick (2004) are 8.33 ± 0.21, 7.61 ± 0.35, 6.22 ± 0.25, 4.67 ± 0.34, and 6.06 ± 0.25, respectively. These values are in line with the derived nebular abundances (Table 7). Our adopted D = 11.33 kpc is in agreement with Frew et al. (2016), who reported 14.86 ± 4.26 kpc. Our derived D is also comparable with the value (14.31 ± 8.54 kpc) determined from the parallax measured using Gaia DR2 (σ_π = 0.0699 ± 0.0417 mas; Gaia Collaboration et al. 2018). Thus, we simultaneously justify our estimated D and nebular ε(O/Ne/Cl/S/Ar).

6.1.2 Central star

As input to Cloudy, we use the tlusty synthetic spectrum of the central star, adopting the parameters from Table 5. In our iterations here, we only vary T_eff in the range of |$27\, 800-28\, 500$| K and L_* in the range of 7300 − 8300 L_|$\odot$|.

6.1.3 Nebula geometry and boundary condition

We plot the observed data and its interpolated curve in Fig. 6. The continuum spectrum in the wavelength ≲ 0.36 |$\mu$|m corresponds to the sum of (1) the nebular continuum (green line in Fig. 5(a)) and (2) the synthetic spectrum of the central star (red line in Fig. 5(b)). The far-IR flux density at 65, 90, 100, 120, and 140 |$\mu$|m is an expected value obtained by fitting for the 15−40 |$\mu$|m Spitzer/IRS spectrum.

We use equations (2) and (3) by Otsuka et al. (2014), who fitted the Spitzer/IRS spectra of Galactic C₆₀ PNe in   15−40|$\mu$|m with the synthetic absorption coefficient Q_λ value. For  SaSt2-3, we set the minimum dust temperature  = 20 K and adopt p = q = 2 and α = 0 which are the same values used in Otsuka et al. (2014). From the fitting, we derive the maximum dust temperature of 136.2 ± 0.4 K, the expected F_ν at 65, 90, 100, 120, and 140 |$\mu$|m are 105.3, 63.5, 51.6, 35.0, and 24.6 mJy, respectively.

From integrating this SED (i.e., the orange region in Fig. 6(a)), we find a total luminosity of ∼8215 L_|$\odot$|. The luminosity of each component is

• Central star: ∼7765 L_|$\odot$|,

• Nebular continuum  + dust continuum: ∼392 L_|$\odot$|,

• Nebular emission line: ∼58 L_|$\odot$|.

Here, the luminosity of nebular emission line is the sum of all the detected emission lines in the HDS and Spitzer spectra. We obtain L_* of ∼7765 L_|$\odot$| by integrating flux density of the central star’s SED within the same wavelength range, indicated by the red region in Fig. 6(b). Thus, we find that only ∼6 per cent of the central star’s radiation (= (392 + 58)/7765) seems to be absorbed by the nebula.

To check whether nebula boundary is determined by the power of the stellar radiation (i.e. ionization bound) or material distribution (i.e. material bound), we tested both the ionization boundary model and the material boundary model. The WHIRC Brγ and Brγ45 images in Fig. 7 display the central bright region and the compact nebula. From the Brγ – Brγ45 image (Fig. 7(c)), we measure the radius of the ionized nebula extended up to be 1.2″.

The material boundary model (the model calculation is stopped at the outer radius of 1.2″) gives a good fit except for underestimates of the [O i] and [N i] line fluxes. However, when we adopt an open geometry such as a cylinder in the ionization boundary model (the model calculation is stopped when T_e is dropped down to ∼4000 K where most of the ionized species are not emitted), we explained well the balance between the input energy from the central star and the output energy from the nebula plus dust (i.e. L_* ≫ L_{neb + dust}) and the observed [O i] and [N i] line fluxes. Although the WHIRC images do not clearly show a cylinder or bipolar nebula, it is plausible judging from the [O i] 6300 Å line-profile. Fig. 8 shows the [O i] 6300 Å line-profile fitting by two Gaussian components with v_r  = +154.6 and  +170.6 km s⁻¹ at the peak intensity of each component. The [O i] lines emitted from the most outer part of the nebula show blue-shifted asymmetry. Such asymmetric profiles are seen in e.g. PN Wray 16-423 (Otsuka 2015); Wray16-423 has a bright cylindrical structure surrounded by an elliptically extended nebula shell.

From these discussions, we adopt the cylinder geometry with the height  = 0.8″. We determined this scale height through a small grid model, and we found that the cylinder height ≥ 0.8″ is necessary. We adopt ionization bounded condition, assuming the ionization front radius of 1.2″.

6.1.4 Elemental abundances and hydrogen density

We adopt the nebular value of ε(N/O/Ne/S/Cl/Ar/Fe) (Table 7) as the initial value and then refine via model iterations within 0.2 dex of the input values so that the best-fit abundances would reproduce the observed emission line intensities. We adopt the nebular ε(He) = 10.96 as the first guess, and vary it in range from 10.75 to 11.10. We keep an expected CEL ε(C) = 8.58 (c.f. stellar ε(C) = 8.55) and stellar ε(Si) = 6.81 through the model iterations. For α-elements Mg, Ca, and Ti (not derived, though), we fix the [Mg,Ca,Ti/Fe]  = +0.3, where [Fe/H] is −1.1 (Section 5). We adopt a constant hydrogen number density (n_H) radial profile. We first guess that n_H is equal to n_e; we adopt the average n_e amongst the measured n_e except for n_e([N i]) (Table 2), then we vary n_H to get the best fit.

6.1.5 Dust grains

We assume that the underlying continuum is due to graphite grains based on the fact that  SaSt2-3 shows the spectral signature of carbon-rich species.

We use the optical data of Martin & Rouleau (1991) for randomly oriented graphite spheres, and we assume the ‘1/3 − 2/3’ approximation (Draine & Malhotra 1993). We adopt the grain radius a = 0.05 − 0.25 |$\mu$|m and a^−3.5 size distribution. If we set the smallest a = 0.005 |$\mu$|m, the maximum grain temperature is over the sublimation temperature of 1750 K. Thus, we set the smallest a = 0.05 |$\mu$|m. We resolve the size distribution into 20 bins (the smallest is 0.05 − 0.054 |$\mu$|m and the largest is ∼0.23 − 0.25 |$\mu$|m). We do not attempt to fit the broad   6−9|$\mu$|m and 11 |$\mu$|m features because the carriers of these features are not determined yet and also these profiles are different from typical band profile of the polycyclic aromatic hydrocarbons (PAHs) in the same wavelengths.

6.2 Modelling results

The input parameters of the best fitting result and the derived quantities are summarized in Table 8. In total, we varied 12 parameters within a given range; L_*, ε(He/N/O/Ne/S/Cl/Ar/Fe), inner radius (r_in), n_H, and grain abundance until |$\chi _{\nu }^{2}$| calculated from I(Hβ), 76 emission line fluxes, 25 broadband fluxes, 4 mid-IR flux densities, and ionization-bound radius (i.e. outer radius r_out). Since there is no observed far-IR data, we stop the model calculation at the ionization front, where T_e is dropped down to ∼4000 K. To evaluate the goodness of the model fitting, we refer to |$\chi _{\nu }^{2}$|⁠. For the [O i] and [N i] lines, we adopt 30 per cent relative uncertainty because these lines are mostly from the PDRs. For the higher order Balmer lines H i (B24 - B14), we set 10 per cent relative uncertainty by considering into account the uncertainty of these lines largely affected by the stellar absorption. The reduced-χ² value in the best model is 12. The relatively large reduced-χ² value even in the best fitting would be due to the uncertainty of the atomic data which we cannot control. Therefore, we conclude that our best fitting result reproduces observations very well. The predicted line fluxes, broadband fluxes, and flux densities are compiled in Appendix Table A5. For references, we list expected F_ν at 65, 90, 100, 120, and 140 μm obtained by fitting for the   15−40|$\mu$|m Spitzer/IRS spectrum (Section 6.1.3). In Figs. 9 and 10, we compare the model SED with the observed one.

In Fig. 11, we show the location of the CSPN predicted by our Cloudy model on the post-AGB evolutionary tracks for initially Z = 0.001 and 1.0, 1.25, and 1.50 M_|$\odot$| (Vassiliadis & Wood 1994). We generate this 1.25 M_|$\odot$| track by linear interpolation between the 1.00 and 1.50 M_|$\odot$| tracks of (Vassiliadis & Wood 1994). Our Cloudy model predicts L_* = 7400 L_|$\odot$|, |$T_{\rm eff} = 28\, 170$| K, and m_* = 0.611 M_|$\odot$|. L_* and m_* are justly close to predicted values (∼7765 L_|$\odot$| and 0.649 M_|$\odot$|) based on the models of initially 1.25 M_|$\odot$| stars with Z = 0.001 (Vassiliadis & Wood 1994; Fishlock et al. 2014). In Fig. 11, we plot the model results of Gesicki & Zijlstra (2007) and Otsuka et al. (2014) as well. These two models show a large discrepancy from the predicted post-AGB evolution track, but our model completely improves this.

Our model succeeds in reproducing the derived ε(X), the volume average T_e (9150 K, while 8930 K in the observation), and I(Hβ). In Table 9, we present the fraction of each ion in each element. Except for C and Ar, the model predicted ICF is well consistent with the empirically determined ICF.

The gas mass m_g is the sum of the ionized and neutral atomic/molecular gas species. Note that we stopped the model calculation at the ionization front. Our m_g is ∼18 per cent of the ejected mass at the last thermal pulse (TP) of 1.25 M_|$\odot$| initial mass stars with Z = 0.001 (0.334 M_|$\odot$|, Fishlock et al. 2014). If we increase the emitting volume by adopting a closed-geometry such as a spherical nebula, the situation is slightly improved (we trace ∼35 per cent of the ejected mass) but the fitting model becomes worse as we explained. According to Fishlock et al. (2014), such stars experienced the superwind phase in the final few TPs during which the mass-loss rate reaches a plateau of ∼10⁻⁵ M_|$\odot$| yr⁻¹. One might think that our underestimated m_g might be caused by excluding the neutral gas and molecular gas regions. However, it is unlikely that  SaSt2-3 has the molecular gas rich envelope because the molecular hydrogen H₂ lines in K-band are not detected (Lumsden et al. 2001). Indeed, we confirm this fact by analysis of the UKIRT CGS4 HK-band spectrum of  SaSt2-3 (Fig. 1). Since our m_g is greater than the ejected mass 2.7(–3) M_|$\odot$| at the last TP of initially 1.00 M_|$\odot$| stars with Z = 0.001, we can conclude that the progenitor should be a >1.00 M_|$\odot$| initial mass star. Fishlock et al. (2014) predicts that 1.00 M_|$\odot$| stars lost the majority of its stellar envelope before it reaches the superwind phase. The 1.25 M_|$\odot$| progenitor that we infer might have experienced such mass loss. If this is true, non-detection of the H₂ lines might be because thin circumstellar envelope does not shield UV radiation from the CSPN and H₂ is dissociated. As an other explanation for the estimated small m_g, the ejected mass during the AGB phase might be efficiently transported to the stellar surface of a companion star. Our small m_g would be largely improved by taking cold gas/dust components that can be traced by far-IR observation.

As presented in Fig. 10(b), our Cloudy model predicts an emission line around 7 |$\mu$|m. This line is the complex of the [Ar ii] 6.99 |$\mu$|m and H i 6.95/7.09 |$\mu$|m lines. In the Spitzer/IRS spectrum, we measure the total line-flux of the [Ar ii] 6.99 |$\mu$|m, H i 6.95/7.09 |$\mu$|m, and C₆₀ 7.0 |$\mu$|m to be 28.67 ± 3.23, where I(Hβ)  = 100. The model predicts I([Ar ii] 6.99 |$\mu$|m  + H i 6.95/7.09 |$\mu$|m)/I(Hβ)  = 4.173 (I(Hβ)  = 100). The contribution of the atomic line complex to the C₆₀ 7.0 |$\mu$|m band (14.6 |${{\ \rm per\ cent}}$|⁠) is not as significant as Otsuka et al. (2014) expected (30.3 |${{\ \rm per\ cent}}$|⁠). We estimate I(C₆₀ 7.0 |$\mu$|m)/I(Hβ) to be 24.5 ± 2.8.

7 DISCUSSIONS

7.1 Evolution of  SaSt2-3

Using the Galactic rotation velocity based on the distance scale of Cahn, Kaler & Stanghellini (1992), van de Steene & Zijlstra (1995), and Zhang (1995), v_r(LSR)  = +149.1 km s⁻¹, and D = 11.33 kpc, we obtain ΔV = 68.2 − 77.9 km s⁻¹. The height above the Galactic plane |z| is 1.13 kpc. These results are in agreement with the Type III PN classification of Quireza et al. (2007). ΔV and |z| do not exceed 120 km s⁻¹ and 1.99 kpc for the Type IV PN classification of Quireza et al. (2007), respectively. The metallicity of  SaSt2-3 is much richer than typical halo PNe such as  K648,  BoBn1, and  H4-1 showing [Ar/H] ≲ –2 (Otsuka et al. 2010; Otsuka, Hyung & Tajitsu 2015; Otsuka & Tajitsu 2013). Thus, we conclude that  SaSt2-3 belongs to the thick disc younger population and a Type III PN rather than a Type IV PN (Pereira & Miranda 2007). Note that classification of PN type does not matter whether the central star is binary or not.

SaSt2-3 would have evolved from a binary composed of a ∼1.25 M_|$\odot$| initial mass star and an companion star. However, any parameters on binary motion are unknown yet. According to the simulation using the binary_c ¹³ code by Izzard et al. (2004), an initially 1.25 M_|$\odot$| single star with Z = 0.001 will enter the PN phase within 3.5 Gyr after the progenitor was in the main-sequence. Perhaps, evolutionary time required to reach the AGB phase would be shortened by binary interaction (i.e. < 3.47 Gyr).

SaSt2-3 composes of a ∼0.61 M_|$\odot$| B-type cool central star with |$T_{\rm eff} = 28\, 100$| K and a companion star. We find the Ca I absorption centred at 6616.62 Å and 6689.28 Å in heliocentric wavelength (6613.13 and 6685.6 Å in rest wavelength, respectively). v_r using these two absorption lines is +165.5 km s⁻¹ and +152.1 km s⁻¹, respectively. These v_r is close to v_r of the central star (+183.6 km s⁻¹, Table 6). Thus, we assume that these lines could be originated from the envelope in the companion star. However, since there is only one spectrum covering 6600 − 6690 Å we have (Table 1), we do not yet find radial velocity variations of this absorption. We suppose that the companion star might be a F-K spectral type star in the main-sequence in terms of the initial mass by referring to De Marco et al. (2013).

7.2 Comparisons with non-C₆₀ and C₆₀-containing PNe

Our study can fully characterize the physical properties of  SaSt2-3. Thus, we are able to compare nebular elemental abundances and central star properties with those of other C₆₀ PNe, and we attempt to gain insights into the C₆₀ formation in PNe. For this purpose, we select C₆₀ PNe  Tc1 and  IC418 because they were previously modelled using Cloudy, and they have been extensively studied. In Table 10, we compile their properties. Due to the lack of existing UV spectra, the CEL ε(C) in  SaSt2-3 is not determined yet. However, since we adopt an expected CEL ε(C) (= (C/O)_* × O_CEL) from the stellar C/O ratio, the reliability of the discussion here is not compromised. We adopt the results of  Tc1 and  IC418 by Pottasch et al. (2011) and Morisset & Georgiev (2009), respectively. For  Tc1, since Pottasch et al. (2011) calculated ε(Ar) by the sum of the Ar⁺ and Ar²⁺ abundances without subtracting C₆₀ 7.0 |$\mu$|m and H i lines from the complex line at 7 |$\mu$|m, their calculated ε(Ar) is certainly overestimated. Therefore, we compute ε(Ar) from Ar²⁺ (7.0(–7)) and ICF(Ar)  = S/S²⁺ (1.71). ε(He) in  Tc1 is predicted by their Cloudy model.

Comparisons with theoretical AGB nucleosynthesis models of Karakas (2010) indicate that the initial mass (M_ini.) is ∼1.90 − 2.10 M_|$\odot$| for  IC418 and ∼1.50 − 1.90 M_|$\odot$| for  Tc1, respectively. From plots of L_* and T_eff on the post-AGB evolutionary tracks based on Vassiliadis & Wood (1994) (Fig. 12), we have the same estimate for  IC418 (∼1.8 M_|$\odot$|) and  Tc1 (∼1.5 M_|$\odot$|). Based on our estimates for the initial mass and metallicity, the age of  IC418 and  Tc1 after the main sequence is ∼2 − 3 Gyr. Thus, we attest that  SaSt2-3 is the most metal-deficient and oldest Galactic C₆₀ PN.

Our findings in elemental abundances are as follows; (i) the values of ε(C) in these C₆₀ PNe are not peculiar values that can be explained by the AGB nucleosynthesis models, (ii) despite that the C₆₀ band strength in  SaSt2-3 and  Tc1 is much stronger than that in  IC418 (see Table 10 of Otsuka et al. 2014), and ε(C) in  Tc1 and  SaSt2-3 is smaller than that in  IC418, and (iii) the C/O ratio (an indicator of the amount of the C-atoms unlocked in dust and molecules) in  Tc1 is smaller than that in  IC418. Supporting (i), non C₆₀-containing C-rich PNe  IC2165 (Miller et al. 2019) and  Me2-1 (Pottasch & Bernard-Salas 2010) (they are selected based on metallicity and initial mass) show similar abundances (including ε(C)) to C₆₀ PNe. The difference between non-C₆₀ and C₆₀ PNe is T_eff only, supporting that the weak radiation field from the central star is in favour of the C₆₀ formation. The time during which C₆₀ is present might be a short-lived phase that C-rich PNe go through, but we do not yet have firm observational evidence of this. These findings would result in our conclusion that the C₆₀ formation does not largely depend on the amount of the C-atoms produced during the AGB phase, and  Tc1 and  SaSt2-3 efficiently produced the C₆₀ molecule by some mechanisms not present in  IC418.

It is noticeable that Fe abundance is highly depleted in all C₆₀ PNe. The highly deficient Fe is possibly due to selective depletion in a binary disc (e.g. Otsuka et al. 2016, reference therein). If C₆₀ PNe have a disc around the central star, they can harbour mass-loss including AGB products for a long time and shield C₆₀ molecules from the intense central star’s UV radiation. Then, most of the Fe-atoms might be tied up in dust grains (e.g. FeO) within a disc. Accordingly, an environment suitable for large carbon molecule formation might be created. This would be in the case of  SaSt2-3. The central star of  IC418 would be not a binary (e.g. Méndez 1989). There are no reports on the binary central star of  Tc1 so far.

If it is the case, is C₆₀ formation dependent on the central star’s properties and its surrounding environment, such as a binary disc? To answer this question, we need to calculate the mass of the C-atoms present in atomic gas, dust and C₆₀ by a fair means. From a more global perspective, the fraction of the C-atoms would be a critically important parameter to understand how much mass carbon PN progenitors had returned to their host galaxies. However, we find that the excitation diagram based on the observed C₆₀ 8.5/17.4/18.9 |$\mu$|m fluxes and the expected C₆₀ 7.0 |$\mu$|m flux (I(C₆₀ 7.0 |$\mu$|m)/I(Hβ)  = 24.5 ± 2.8, see Section 6.2) in  SaSt2-3 indicates non-LTE conditions. The same situation exists in other C₆₀ PNe including  Tc1 (Cami et al. 2010) and  IC418 (Otsuka et al. 2014). Therefore, we sought other ways not involving the excitation diagrams. Of these, the method proposed by Berné & Tielens (2012) is seemingly suitable to our aim; this method requires the observed/modelled IR SED, all four mid-IR C₆₀ band fluxes, and dust mass. These three parameters are already determined in  SaSt2-3 and  IC418.

In Table 11, we summarize the mass of the C-atoms present in atomic gas, dust, and C₆₀ in  SaSt2-3,  IC418, and  Lin49 in the SMC (as a comparison of the possible binary C₆₀ PN). For  SaSt2-3 and  Lin49, we obtain I(IR) by integrating the atomic gas emission free SED (generated by the Cloudy best model) in the range from 4 to 200 |$\mu$|m. Then, we calculate the total C₆₀ flux I(C₆₀) (Table A1). Finally, assuming that the dust and C₆₀ are emitted in the same regions, we estimate the C₆₀ mass (=m_d × I(C₆₀)/I(IR)). The column density and total number of C₆₀ in  SaSt2-3 are 8.57(+12) cm⁻² and 7.53(+47), respectively. By applying the method to  IC418 (but using the combined Infrared Space Observatory (ISO)/Short Wavelength Spectrometer (SWS)/Long-Wave Spectrometer (LWS) and Spitzer combined spectrum), we obtain their C₆₀ mass.

From this analysis, we quantitatively demonstrate that  SaSt2-3 and  Lin49 produced C₆₀ more efficiently than  IC418. Thus, we conclude that the C₆₀ formation depends upon the central star’s properties and its surrounding environment (e.g. a binary disc), rather than the amount of C-atoms.

8 SUMMARY

We have studied the fullerene-containing PN  SaSt2-3 in order to investigate its physical properties and gain insights into the C₆₀ formation in PN progenitors. We derived the nine and four elemental abundances from nebular line and stellar absorption analysis, respectively. The derived elemental abundances indicate that the progenitor is an initially ∼1.25 M_|$\odot$| star with Z = 0.001 and α-element and Cl enhanced ([α,Cl/Fe] ∼ +0.3–0.4). The distance D of 11.33 kpc is determined by comparing the derived luminosity as a function of D and T_eff with the predicted luminosity by the post-AGB evolution model for the 1.25 M_|$\odot$| stars with Z = 0.001.  SaSt2-3 is classified as part of the thick disc population with an early B-type central star with |$T_{\rm eff} = 28\, 100$| K, |$\log \, {g} = 3.11$| cm s⁻², and the core-mass of ∼0.61 M_|$\odot$|. We discovered the binary central star of  SaSt2-3 from time-variation of the stellar radial velocity. Further observations are necessary to understand the binary system. We built the comprehensive photoionization model. The calculated gas mass is much smaller than the AGB model prediction for the single 1.25 M_|$\odot$| stars with Z = 0.001. The lower gas mass could be due to the short duration time of the superwind phase or efficient mass-transfer into the stellar surface of the companion star during the AGB phase. From the simple analysis, we quantitatively demonstrate that  SaSt2-3 produced C₆₀ more efficiently than other C₆₀ PNe. The C₆₀ formation would depend on the central star’s properties and its surrounding environment. There might be a link between the C₆₀ formation efficiency and the binary central star. Spatially resolved spectral maps of the atomic carbon, carbon dust, and C₆₀ are necessary to identify the locations of C₆₀ and investigate the abundance distribution of these species within dusty nebula. We succeed to demonstrate what type of and how much mass of stars can produce how much C₆₀ molecules. We will further investigate stellar/nebular properties and C₆₀ in order to find what parameters are of critical importance in the C₆₀ formation.

ACKNOWLEDGEMENTS

I am grateful to the anonymous referee for carefully reading and the useful suggestions which greatly improved this article. I learned a lot of things from his/her comments. I was supported by the research fund 104-2811-M-001-138 and 104-2112-M-001-041-MY3 from the Ministry of Science and Technology (MOST), R.O.C. I thank Dr. Akito Tajitsu for supporting my Subaru HDS observations. I sincerely thank Drs. Benjamin Sargent, Peter Scicluna, and Toshiya Ueta for critically reading the paper and giving suggestions. This work was partly based on archival data obtained with the Spitzer Space Telescope, which is operated by the Jet Propulsion Laboratory, California Institute of Technology under a contract with NASA. This research is in part based on observations with AKARI, a JAXA project with the participation of ESA. A portion of this work was based on the use of the ASIAA clustering computing system.

Footnotes

REFERENCES

APPENDIX A: SUPPORTING RESULTS

The following tables and figure support our works.