2023/07/31 2023/09/23

planetary nebulae: individual (H4-1) – ISM: abundances – dust, extinction –
stars: Population II

#  Seimei KOOLS-IFU mapping of the gas and dust distributions in Galactic PNe:
Unveiling the origin and evolution of Galactic halo PN H4-1

Masaaki Otsuka11affiliation: Okayama Observatory, Kyoto University, Honjo,
Kamogata, Asakuchi, Okayama 719-0232, Japan Corresponding author
<EMAIL_ADDRESS>Toshiya Ueta 22affiliation: Department of
Physics and Astronomy, University of Denver, 2112 E Wesley Ave., Denver, CO
80208, USA<EMAIL_ADDRESS>Akito Tajitsu33affiliation: Subaru Telescope
Okayama Branch Office, National Astronomical Observatory of Japan, National
Institutes of Natural Sciences, 3037-5 Honjo, Kamogata, Asakuchi, Okayama
719-0232, Japan<EMAIL_ADDRESS>
###### Abstract

H4-1 is a planetary nebula (PN) located in the Galactic halo, and is notably
carbon-rich and one of the most metal-deficient PNe in the Milky Way. To
unveil its progenitor evolution through the accurate measurement of the gas
mass, we conducted a comprehensive investigation of H4-1, using the newly
obtained Seimei/KOOLS-IFU spectra and multiwavelength spectro-photometry data.
The emission line images generated from the KOOLS-IFU datacube successfully
resolve the ellipsoidal nebula and the equatorial flattened disk that are
frequently seen in bipolar PNe evolved from massive progenitors. By a fully
data-driven method, we directly derived the seven elemental abundances, the
gas-to-dust mass ratio, and the gas and dust masses based on our own distance
scale. By comparing the observed quantities with both the photoionization
model and the binary nucleosynthesis model, we conclude that the progenitors
of an initial mass of 1.87 M⊙ and 0.82 M⊙ are second generation stars formed
$\sim 4$ Gyrs after the Big Bang, and underwent mass-transfers, binary merger,
and ultimately evolved into a PN showing unique chemical abundances. Our
binary model successfully reproduces the observed abundances and also explains
evolutionary time scale of H4-1.

##  1 Introduction

Table 1: Observed nebular elemental abundances of the halo PNe.

PN | He | C (CEL) | N | O | Ne | S | Cl | Ar | Fe | Ref.
---|---|---|---|---|---|---|---|---|---|---
H4-1 | 11.02 | 8.86 | 7.74 | 8.17 | 6.42 | 5.14 | 3.88 | 4.61 | $\cdots$ | (1), (2)
K648 in M15 | 11.02 | 8.97 | 6.36 | 7.73 | 7.44 | 5.40 | 3.58 | 4.60 | 5.02 | (3)
BoBn1 | 11.07 | 9.02 | 8.03 | 7.74 | 7.96 | 5.32 | 3.39 | 4.33 | 5.08 | (4)
Average | 11.04 | 8.95 | 7.38 | 7.88 | 7.27 | 5.29 | 3.62 | 4.51 | 5.05 |
DdDm1 | 11.01 | 7.02 | 7.36 | 8.07 | 7.45 | 6.35 | 4.61 | 5.79 | 6.58 | (5)
NGC2242 | 11.00 | 8.39 | 7.72 | 8.03 | 7.80 | 6.33 | 4.45 | 5.89 | $\cdots$ | (6),(7)
NGC4361 | 11.03 | 8.13 | 7.41 | 7.82 | 7.54 | $\cdots$ | $\cdots$ | 6.01 | $\cdots$ | (7)
PN G135.9+55.9 | 10.95 | 7.84 | 7.15 | 6.82 | 6.83 | 5.50 | $\cdots$ | 4.50 | $\cdots$ | (8)
PRMG1 | 10.96 | $\cdots$ | $\cdots$ | 8.10 | 7.50 | $\cdots$ | $\cdots$ | 5.80 | $\cdots$ | (9)
PRMT1 | 11,03 | 7.60 | 8.00 | 8.30 | 7.90 | $\cdots$ | $\cdots$ | 6.30 | $\cdots$ | (10)
M2-29 in M22 | 10.97 | $\cdots$ | 6.98 | 7.31 | 6.72 | 5.91 |  | 5.26 | $\cdots$ | (11)
JaFu1 in Pal6 | 11.15 | $\cdots$ | 7.97 | 8.49 | 8.00 | 6.62 | $\cdots$ | $\cdots$ | $\cdots$ | (12)
JaFu2 in NGC6441 | 11.06 | $\cdots$ | 7.96 | 7.73 | 6.79 | 6.77 | $\cdots$ | 5.46 | $\cdots$ | (12)
Average | 11.02 | 7.80 | 7.57 | 7.85 | 7.39 | 6.25 | 4.53 | 5.63 | 6.58 |
Sun | 10.93 | 8.43 | 7.83 | 8.69 | 7.93 | 7.12 | 5.50 | 6.40 | 7.50 | (13)

Note – The C abundances are derived from collisional excitation lines (CELs).
The S abundance in NGC2242 is calculated by us using the data of Kwitter et
al. (2003). There are no reports of the nebular elemental abundances of GJJC1
in M22.
References – (1) This work; (2) Otsuka & Tajitsu (2013); (3) Otsuka et al.
(2015); (4) Otsuka et al. (2010); (5) Otsuka et al. (2009); (6) Kwitter et al.
(2003); (7) Torres-Peimbert et al. (1990); (8) Stasińska et al. (2010); (9)
Peña et al. (1989); (10) Pena et al. (1990); (11) Pena et al. (1991); (12)
Jacoby et al. (1997); (13) Asplund et al. (2009).

Planetary nebulae (PNe) are the final stage of evolution of low to
intermediate mass stars with initial mass of $1-8$ M⊙ (e.g., Kwok, 2000). In
our Galaxy, more than 3000 PNe have been identified so far (Parker et al.,
2016), mostly belong to the Galactic disk. However, there are also 13 halo
PNe111Otsuka et al. (2015) reported 14 objects as halo PNe by following
Pereira & Miranda (2007) reported that SaSt2-3 is a new halo PN. Later, Otsuka
(2019) concluded that this PN is not halo but thin disk population.,
comprising eight field PNe in the Galactic halo and five members of globular
clusters.

Halo PNe are classified into type-IV on the basis of kinematical/chemical
characteristics which was originally proposed by Peimbert (1978) and later
redefined by Quireza et al. (2007). Halo PNe have low metallicity, with the
average [Ar/H] of $-1.33$ among 12 halo PNe (table 1). Despite their scarcity,
halo PNe have been considered special because (i) they were formed in the
early stages of Galaxy formation, and (ii) they allow us to directly measure
both chemical abundances and core masses of the central star after the
asymptotic giant branch (AGB) star nucleosynthesis. To better understand halo
PNe, we have investigated physical properties of K648 in the globular cluster
M15 (Otsuka et al., 2015), BoBn1 (Otsuka et al., 2008, 2010), DdDm1 (Otsuka et
al., 2009), and H4-1 (Otsuka et al., 2003; Otsuka & Tajitsu, 2013). In this
paper, we discuss the origin and evolution of H4-1.

H4-1 is one of the most metal-poor ($\mathrm{[Ar/H]}=-1.78\pm 0.15$,
corresponding to the metallicity of $Z\sim 10^{-4}$) and carbon-rich
($\mathrm{C/O}=4.90$) PNe in the Milky Way, and hence, is of special interest
along with K648 and BoBn1, which share similar chemical properties (table 1).
Also, H4-1 is remarkably $r$-process rich ($\mathrm{[Xe/H]}>+0.51$; Otsuka &
Tajitsu, 2013). Given the low metallicity, H4-1 is believed to have formed
$\sim 10-13$ Gyrs ago. However, our understanding of its evolution remains
incomplete. More specifically, there have been longstanding unresolved issues
on its slow evolutionary time scale and high C abundance. Based on the
isochron fitting to the main-sequence turnoff point of M15 (of $Z=10^{-4}$,
Fahlman et al., 1985), H4-1 was considered to be of 0.8 M⊙ initial mass. Post-
AGB evolutionary models of a 0.85 M⊙ and $Z=10^{-4}$ star by Miller Bertolami
(2016) suggest that such stars would take $\gtrsim 15000$ yrs to increase the
effective temperature ($T_{\rm eff}$) from $\sim 5000$ K (when mass loss is
considered to be terminated; Schoenberner, 1983) to $\sim 30000$ K (when
ionization of the circumstellar nebula begins). In addition, such stars are
expected to spend $>60000$ yrs to reach the current $T_{\rm eff}$ of $\sim
120000$ K from $T_{\rm eff}$ of $\sim 5000$ K (table 5). It is highly unlikely
to visually identify such extremely slow-evolving PNe. However, why can we
still observe H4-1?

According to the stellar evolution models, third dredge-up (TDU) is required
during the thermal pulse (TP) AGB phase in order for the low-intermediate mass
to become carbon-rich. The lower mass limit for TDU to happen is dependent on
model assumptions; as far as we know, it is about 0.9 M⊙ (e.g. Lugaro et al.,
2012), which is more massive than an expected initial mass of H4-1. If so, how
could H4-1 become C-rich?

Is H4-1 really evolved from a single star of about $0.8-0.9$ M⊙ initial mass?
Otsuka & Tajitsu (2013) found out that the [C/O] and [$s$-element/$r$-element]
ratios in H4-1 are similar to those in carbon-enhanced metal-poor (CEMP) stars
(e.g., Beers & Christlieb, 2005) which are believed to originate from binary
systems. Following the classification by Suda et al. (2011) for metal-poor
stars, H4-1 is CEMP-$r/s$ or CEMP-no (Otsuka & Tajitsu, 2013). Binary mass-
transfer and merger would drastically shorten the evolutionary timescale,
enabling an explanation of the evolutionary timescale. As such binary
evolution would increase up to the mass required for TDUs, we would explain
the C-richness of H4-1. Consequently, binary evolution would allow us to
explain both the evolutionary timescale and the carbon enhancement without any
contradictions. Indeed, there is a recent report on the detection of molecular
hydrogen emission lines, which implies that H4-1 originated from a massive
star (Tajitsu et al., 2023). Hence, we need to seek more conclusive evidence
of its binary origin based on observational quantities.

One way to determine the origin and evolution of H4-1 is to estimate the total
mass of the circumstellar matter. Miller Bertolami (2016) predicts that a
single star of initially $0.85\,\mathrm{M}_{\odot}$ with $Z=10^{-4}$ would
become a $0.53\,\mathrm{M}_{\odot}$ white dwarf. That is, such a star would
eject the circumstellar mass of $0.32\,\mathrm{M}_{\odot}$ during its
evolution. If the progenitor is a binary system consisting of a $\sim
0.8\,\mathrm{M}_{\odot}$ star and a more massive star, a greater amount of
mass would have to be ejected. The origin and evolution of H4-1 and its
anomalous abundance pattern will also be important to understand the nature of
CEMP stars.

To accurately estimate the total mass of the circumstellar matter, and
subsequently, the mass of the progenitor system, an accurate distance estimate
is critical. However, the distance toward H4-1 has not been fully established
(see subsection 7.2). The distance to H4-1 appears too far to be measured even
with Gaia (González-Santamaría et al., 2021). Therefore, we undertake the
distance determination on our own, which requires multi-wavelength spectro-
photometry data and post-AGB evolutionary tracks.

This paper is organized as follows. In sections 2 and 3, we explain a new set
of data obtained with our own Seimei/KOOLS-IFU observations, the adopted
archival data, and how the data are reduced, before presenting the obtained
spectrum and emission line images. Then, a fully data-driven plasma
diagnostics with self-consistent dust extinction correction is performed in
section 4, followed by a detailed abundance analysis in section 5. The gas and
dust masses and the gas-to-dust mass ratio are calculated in section 6. We
have discussions on the evolution of the progenitor in section 7, and
summarize the current work in section 8.

Figure 1:  One of the dither positions of the packed fiber array in KOOLS-IFU
observations overlaid on the Subaru 8.2-m/Suprime-Cam image presented in
Tajitsu et al. (2023). The intensity scale ranges between the minimum and the
maximum. The spatial resolution of the image is indicated by the gray circle.
Each hexagon indicates the position and dimension of the 110 fibers.

##  2 Dataset and reduction

We explain the used dataset and the data reduction below. The logs of the
spectroscopic observations are summarized in table 2.

###  2.1 Seimei/KOOLS-IFU optical spectra

We performed the mapping observation (figure 1) using the Kyoto Okayama
Optical Low-dispersion Spectrograph with optical-fiber integral field unit
(KOOLS-IFU; Matsubayashi et al., 2019) attached to one of the two Nasmyth foci
of the Kyoto University Seimei 3.8-m telescope (Kurita et al., 2020).

The observations were carried out under clear and stable sky conditions on
three nights in early 2023. The seeing was $\sim 2{\arcsec}$ in $R$-band at
full width at the half maximum (FWHM). The data were taken under the airmass
of $1.10-1.00$ in the VPH-blue and $1.06-1.01$ in the VPH-red settings. We
took $12\times 600$ sec (2023 March) and $4\times 1200$ sec (2023 April) in
the blue and $8\times 900$ sec in the red observations, respectively. The logs
of the spectroscopic observations are summarized in table 2.

For the background sky subtraction, the off-source frames were taken every
three science frames. To perform flux calibration and remove telluric
absorption, we observed Hz44 every 30 min. This standard star was also used to
evaluate light loss in the KOOLS-IFU field and synthesize a point-spread
function (PSF). We also took instrumental flat, twilight sky, and Hg/Ne/Xe arc
frames as the baseline calibration.

We reduced the data using our own codes and IRAF (Tody, 1986) by following the
procedure elucidated by Otsuka (2022). Wavelength calibration was performed
using over 40 Hg/Ne/Xe lines that cover the full range of effective
wavelengths. We adopted a constant plate scale of 2.5 Å. The resultant RMS
errors in wavelength determination were 0.11 Å in the VPH-blue and 0.13 Å in
the VPH-red, respectively. The average spectral resolution ($R$, defined as
$\lambda$/Gaussian FWHM) measured from the Hg/Ne/Xe lines was $\sim 500$ and
$\sim 800$ in the VPH-blue and VPH-red spectra, respectively.

We obtained an RA-DEC-wavelength datacube using our codes and image
reconstruction based on the drizzle technique. We corrected the spatial
displacements due to atmospheric dispersion at each wavelength by following
Owens (1967) and Otsuka (2022). We adopted a constant plate scale of 0.40 per
pixel, which posses good signal-to-noise ratio in each spaxel. We evaluated
the off-center value in RA and DEC relative to a reference position by
measuring the center coordinate of the standard star in the 3-D datacube: the
offset did not exceed 0.02 (0.05 pixel) in both directions over the
wavelengths.

Figure 2: Selected PSF-deconvolved KOOLS-IFU emission-line images. The intensity range is between the minimum and the maximum. The achieved spatial resolution is indicated by the gray circle.  Table 2: Summary of spectroscopic observations.  Date | Telescope/Inst. | Range () | $R$
---|---|---|---
2023/03/13 | Seimei/KOOLS-IFU | $0.420~{}-~{}0.853$ | $\sim 500$
2023/04/10 | Seimei/KOOLS-IFU | $0.420~{}-~{}0.853$ | $\sim 500$
2023/02/11 | Seimei/KOOLS-IFU | $0.560~{}-~{}1.010$ | $\sim 800$
1983/08/03 |  _IUE_ /SWP | $0.115~{}-~{}0.198$ | $\sim 300$
1983/08/03 |  _IUE_ /LWR | $0.185~{}-~{}0.335$ | $\sim 400$
2007/06/17 |  _Spitzer_ /IRS | $5.2~{}-~{}14.5$ | $\sim 100$
2007/06/17 |  _Spitzer_ /IRS | $9.9~{}-~{}37.2$ | $\sim 600$

For flux calibration and telluric removal, we synthesized the non-LTE
theoretical spectrum of Hz44 using the Tlusty code (Hubeny, 1988). We adopted
the basic stellar parameters from Dorsch et al. (2019) who determined surface
gravity $\log_{10}\,g=5.64$ cm s-2 and $T_{\rm eff}=39100$ K. We calculated
the color excess between B and V bands $E(B-V)$ of 0.036 toward Hz44 in the
total-to-selective extinction $R_{\rm V}=3.1$, using a comparison between the
synthesized spectrum and the photometry spectral energy distribution (SED) in
the range from 0.35 to 0.92 generated from Geier et al. (2019) and Aparicio
Villegas et al. (2010). In this procedure, we adopted the reddening law of
Cardelli et al. (1989). The flux calibration and telluric removal were
performed using the function derived by comparing this reddened telluric free
theoretical spectrum and the 1-D spectrum extracted from the observed 3-D
datacube. In the 1-D spectrum extraction, we recovered the lost flux outside
the datacube using the Hz44’s PSF in each wavelength.

We evaluated the achieved spatial-resolution by measuring the Moffat
function’s FWHM of the Hz44’s PSF every wavelength. The Moffat function is
represented by the following:

$\displaystyle I(r)$ $\displaystyle=$
$\displaystyle\left[1+\left(r/{\alpha}\right)^{2}\right]^{-\beta}$ (1)
$\displaystyle{\alpha}$ $\displaystyle=$ $\displaystyle{\rm
FWHM}\left[2^{\left(1/{\beta}\right)}-1\right]^{-0.5},$ (2)

where $I(r)$ is the radial profile of these PSF standard stars, and $\beta$ is
the shape parameter, set to be 4.0 through the wavelengths. To correct the
image quality degradation due to airmass, we adopted the equation established
by Otsuka (2022). Thus, the expected Moffat PSF’s FWHM in arcsec can be
expressed as a function of $\lambda$ in using equation (6).

$\displaystyle{\rm FWHM}=\left\\{\begin{array}[]{ll}(1.70\pm
0.01)\,(\lambda/0.55)^{-0.16\pm 0.01}&\mbox{(Feb)}\\\ (2.28\pm
0.01)\,(\lambda/0.55)^{-0.20\pm 0.01}&\mbox{(Mar)}\\\ (1.88\pm
0.01)\,(\lambda/0.55)^{-0.21\pm 0.01}&\mbox{(Apr)}\end{array}\right.$ (6)

After performing PSF-matching in all spectral bins in the datacube with
equation (6), we carried out spectral analysis.

Figure 3: The KOOLS-IFU 1-D spectrum. Prominent lines among the identified
emission-lines are indicated in blue. The wavelengths are in air and also have
a heliocentric correction. The dark-gray color shade indicates the one-sigma
uncertainty at each wavelength.

###  2.2 _IUE_ UV and _Spitzer_ mid-IR spectra

We downloaded the SWP20599 and LWR16513 spectra (IUE Program ID: FA083, PI: M.
J. Barlow) obtained by the International Ultraviolet Explorer (_IUE_) from the
Mikulski Archive for Space Telescopes (MAST). The size of the slit (9.9
$\times$ 22) covers the entire nebula of H4-1. We combined these two spectra
into a single $1150-3348$ Å spectrum.

We used the mid-IR Spitzer/IRS (Houck et al., 2004) short-low (SL), short-high
(SH), and long-high (LH) spectra used in Otsuka & Tajitsu (2013). The size of
the slit is $3.7\arcsec\times 57\arcsec$ in SL, $4.7\arcsec\times 11.3\arcsec$
in SH, and $11.1\arcsec\times 22.3\arcsec$ in LH. We combined these spectra
into a single spectrum by scaling their flux densities to match with mid-IR
Wide-field Infrared Survey Explorer (_WISE_ ; Cutri et al., 2013) bands 3 and
4. More detailed information on the data and their reduction can be found
therein.

###  2.3 Multiwavelength Photometry

We obtained broadband UV data from _GALEX_ (Bianchi et al., 2017) and _XMM-
Newton_ (Page et al., 2012), optical data from SDSS (Alam et al., 2015), near-
IR data from UKIRT/UKIDSS (Lawrence et al., 2007), mid-IR data from _Spitzer_
/IRAC/MIPS (Fazio et al., 2004; Rieke et al., 2004, for IRAC and MIPS,
respectively) and _WISE_ , and far-IR data from _Herschel_ /PACS (Poglitsch et
al., 2010) as summarized in table A1. The second column provides the
extinction function $f(\lambda)$ at $\lambda$ computed using the reddening law
of Cardelli et al. (1989) for the wavelengths shorter than UKIDSS 2.2 and
Fluks et al. (1994) for the wavelengths longer than IRAC 3.6 band. Herein, we
adopt $R_{\rm V}$ of 3.1. The third and fourth columns in table A1 are the
dust extincted flux density $F_{\lambda}$ and the extinction-free one
$I_{\lambda}$ obtained using equation (7) in section 4, respectively.

##  3 Results of the KOOLS-IFU observations

First, we demonstrate the accuracy of the KOOLS-IFU datacube as photometric
data. We measured the SDSS $r$-band $F_{\lambda}$ using the corresponding band
image synthesized from the KOOLS-IFU data cube and the SDSS-$r$ band filter
transmission. The resultant KOOLS-IFU SDSS $r$-band $F_{\lambda}$ is $(1.75\pm
0.07)\times 10^{-15}$ erg s-1 cm-2 Å-1 while the cataloged SDSS $r$-band
$F_{\lambda}$ is $1.76\times 10^{-15}$ erg s-1 cm-2 Å-1 (table A1). This means
our reduced KOOLS-IFU datacube has sufficient photometric accuracy.

In figure 2, we present the PSF-deconvolved emission-line images. The signal-
to-noise ratio per pixel is $\geq 5$. We generated these images by measuring
the line fluxes in each spectrum every spaxel using the automated line fitting
algorithm code (Alfa; Wesson, 2016). The achieved spatial resolution is
evaluated by the Gaussian FWHM of the PSF deconvolved standard stars.

The ionization stratification in H4-1 appears to be consistent with what is
expected in a gaseous nebula surrounding the central ionization source; that
is, the spatial extent of gas line emission becomes larger as the ionization
energy of the emission becomes lesser. The PSF-deconvolved Gaussian 1-sigma
and the ellipticity of H$\beta$ emission are 0.87 and 0.06, respectively. They
are 1.16 and 0.11 in the [O i] 6300 Å. The nebula is asymmetric and elongated
along position angle of $\sim-25^{\circ}$. These KOOLS-IFU spectral images
successfully resolved the spatially compact yet elongated emission core
structure in the [O ii] 7320/30 Å and [O i] 6300 Å lines. Such structure was
resolved and confirmed only in near-IR H2 images so far (Tajitsu et al.,
2023). The core structure is one of the most distinctive spatial features seen
in H4-1 and is considered to represent the flattened central equatorial disk
in a bipolar PN.

In figure3, we show the spatially-integrated 1-D spectrum generated from the
datacube. The flux density loss is corrected in every wavelength bin using the
synthesized PSF of the standard star. The wavelength is in air and also has a
heliocentric correction; the heliocentric radial velocity measured using all
the detected lines is $-181.8\pm 4.9$ km s-1, which is consistent with Otsuka
& Tajitsu (2013, $-181.35\pm 0.35$ km s-1). In the following sections, we use
the KOOLS-IFU 1-D spectrum.

##  4 Dust Extinction

The observed flux $F$($\lambda$) (and flux density $F_{\lambda}$) is reddened
by not only the interstellar and also the circumstellar dust. We obtain the
extinction-free flux $I$($\lambda$) (and flux density $I_{\lambda}$) using
equation (7):

$\displaystyle I(\lambda)$ $\displaystyle=$ $\displaystyle
F(\lambda)~{}\cdot~{}10^{c({\rm H\beta})(1+f(\lambda))},$ (7)

where $c$(H$\beta$) ($=\log_{10}$ $I$(H$\beta$)/$F$(H$\beta$)) is the line-of-
sight logarithmic reddening coefficient at H$\beta$ 4861.33 Å. Because
$c$(H$\beta$) equals to 1.45 $E(B-V)$ under the assumption of $R_{\rm V}$
being 3.1, one may simply adopt the $E(B-V)$ value measured from the Galaxy
extinction map (e.g., Schlafly & Finkbeiner, 2011); the $E(B-V)$ value toward
H4-1 is 0.0082, which is 0.012 in terms of $c$(H$\beta$). However, this
$E(B-V)$ obtained by the Galactic extinction map is the average global
interstellar medium (ISM) over a large solid angle around H4-1 because the
angular resolution of the map is not high enough. Therefore, we directly
calculate $c$(H$\beta$) by comparing observed H i line fluxes and
corresponding extinction-free H i line-fluxes that are calculated by the H i
emissivity of Storey & Hummer (1995) under the assumption of Case B. It should
be noted that the estimated $c$(H$\beta$) value corresponds to the sum of the
line-of-sight reddening coefficient at 4861.33 Å by the circumstellar and
interstellar media.

To calculate $c$(H$\beta$) using the H i lines, the de-reddened H i line
fluxes must be calculated for comparison. Obviously, the de-reddened H i line
fluxes are computed as functions of $n_{\rm e}$ and $T_{\rm e}$ (Storey &
Hummer, 1995). However, $n_{\rm e}$ and $T_{\rm e}$ are actually the results
of plasma diagnostics, which require extinction-corrected line fluxes. To make
matters even worse, obtaining the accurate $c$(H$\beta$) value in highly
excited nebulae showing He ii lines, we need to subtract the contaminating
neighboring He ii line contribution to H i lines.

This “Chicken or the Egg” dilemma (for details, see Ueta & Otsuka, 2021; Ueta,
2022) and He ii line contamination in $c$(H$\beta$) calculation can be
simultaneously resolved via an iterative loop by equations (8) – (14). The He
ii line subtraction processes are included in equations (8) – (10).

$\displaystyle F({{\rm H}\,{\rm I}}_{\lambda_{i},j})$ $\displaystyle=$
$\displaystyle F({{\rm H}\,{\rm I}}_{\lambda_{i,j}})_{\rm obs}-F({\rm
He}\,{\rm II}\,4686)_{\rm obs}$ (8) $\displaystyle\cdot(I({{\rm He}\,{\rm
II}}_{\lambda_{i,j}})/I({\rm He}\,{\rm II}\,4686))_{(n_{\rm e},T_{\rm e})}$
$\displaystyle\cdot 10^{c({\rm H}\beta)_{i,j}(f(4686)-f({\rm He}\,{\rm
II}_{\lambda_{i,j}}))}$ $\displaystyle F({\rm H}\alpha)$ $\displaystyle=$
$\displaystyle F({\rm H}\alpha)_{\rm obs}-F({\rm He}\,{\rm II}\,4686)_{\rm
obs}$ (9) $\displaystyle\cdot(I({{\rm He}\,{\rm II}}\,6560)/I({\rm He}\,{\rm
II}\,4686))_{(n_{\rm e},T_{\rm e})}$ $\displaystyle\cdot 10^{c({\rm
H}\beta)_{j}(f(4686)-f({\rm He}\,{\rm II}\,6560))}$ $\displaystyle F({\rm
H}\beta)$ $\displaystyle=$ $\displaystyle F({\rm H}\beta)_{\rm obs}-F({\rm
He}\,{\rm II}\,4686)_{\rm obs}$ (10) $\displaystyle\cdot(I({{\rm He}\,{\rm
II}}\,4860)/I({\rm He}\,{\rm II}\,4686))_{(n_{\rm e},T_{\rm e})}$
$\displaystyle\cdot 10^{c({\rm H}\beta)_{i}(f(4686)-f({\rm He}\,{\rm
II}\,4860))}$ $\displaystyle c({\rm H}\beta)_{i}$ $\displaystyle=$
$\displaystyle\frac{1}{f(\lambda_{i})}\log_{10}\frac{(I({{\rm H}\,{\rm
I}}_{\lambda_{i}})/I({\rm H}\beta))_{(n_{\rm e},T_{\rm e})}}{F({{\rm H}\,{\rm
I}}_{\lambda_{i}})/F({\rm H}\beta)}$ (11) $\displaystyle c({\rm H}\beta)_{j}$
$\displaystyle=$ $\displaystyle\frac{1}{f(\lambda_{j})}\log_{10}\frac{(I({{\rm
H}\,{\rm I}}_{\lambda_{j}})/I({\rm H}\alpha))_{(n_{\rm e},T_{\rm e})}}{F({{\rm
H}\,{\rm I}}_{\lambda_{j})}/F({\rm H}\alpha)}$ (12) $\displaystyle c({\rm
H}\beta)$ $\displaystyle=$ $\displaystyle\frac{\sum_{i,j}c({\rm
H}\beta)_{i,j}{\sigma}_{i,j}^{-2}}{\sum_{i,j}{\sigma}_{i,j}^{-2}}$ (13)
$\displaystyle\delta c({\rm H}\beta)^{2}$ $\displaystyle=$
$\displaystyle\sum_{i,j}{\sigma}_{i,j}^{-2},$ (14)

where $F({{\rm H}\,{\rm I}}_{\lambda_{i},j})_{\rm obs}$ is the observed line-
flux of the H i and He ii at $\lambda_{i,j}$, $F({{\rm H}\,{\rm
I}}_{\lambda_{i},j})$ is the sole H i line-flux at $\lambda_{i,j}$,
$I$(HI${}_{\lambda_{i}}$)/$I$(H$\beta$)${}_{(n_{\rm e},T_{\rm e})}$ and
$F$(HI${}_{\lambda_{i}}$)/$F$(H$\beta$) are the He ii contamination-free
theoretical and observed ratios of H i line at $\lambda_{i}$ with respect to
H$\beta$ in $n_{\rm e}$ and $T_{\rm e}$, and
$I$(HI${}_{\lambda_{j}}$)/$I$(H$\alpha$)${}_{(n_{\rm e},T_{\rm e})}$ and
$F$(HI${}_{\lambda_{j}}$)/$F$(H$\alpha$) are the theoretical and observed
ratios of H i line at $\lambda_{j}$ wrt. H$\alpha$ in $n_{\rm e}$ and $T_{\rm
e}$. $I$(He ii${}_{\lambda_{i,j}}$)/$I$(He ii 4686) is the theoretical ratio
of He ii at $\lambda_{i,j}$ wrt. He ii 4686 Å. $\sigma_{i,j}$ is the standard
deviation of $c$(H$\beta$)i,j. We set the searching range of $T_{\rm e}$ and
$n_{\rm e}$ to be $5000-15000$ K and $100-5000$ cm-3.

In practice, due to the lack of the available higher order H i lines which are
sensitive to $n_{\rm e}$, we used [S ii] 6717/31 Å line ratio as the density
constraint. Adopting a genetic algorithm, we obtained the $T_{\rm e}$–$n_{\rm
e}$ pair whose values give the minimum of $\delta$ $c$(H$\beta$)2 after 100
times independent calculations. Finally, we obtained $c$(H$\beta$) = $0.103\pm
0.024$ in $T_{\rm e}$ = $12635\pm 1412$ K and $n_{\rm e}$ = $420\pm 38$ cm-3.

##  5 Chemical Abundances

Table 3: Derived $n_{\rm e}$, $T_{\rm e}$, and ionic/elemental abundances.  Parameter | Value
---|---
$T_{\rm e}$([O iii]) | $12885~{}\pm~{}72\,({\rm K})$
$T_{\rm e}$([N ii]) | $12976~{}\pm~{}46\,({\rm K})$
$T_{\rm e}$(He i) | $7456~{}\pm~{}1774\,({\rm K})$
$n_{\rm e}$([S iii]) | $3074~{}\pm~{}2311\,({\rm cm}^{-3})$
$n_{\rm e}$([S ii]) | $441~{}\pm~{}158\,({\rm cm}^{-3})$
Ion Xi+ ($\\#$ of the used lines) |  $n$(Xi+)/$n$(H+)
He+ (5) | $9.50\times 10^{-2}~{}\pm~{}3.53\times 10^{-3}$
He2+ (3) | $8.81\times 10^{-3}~{}\pm~{}1.74\times 10^{-4}$
C0(CEL) (2) | $2.64\times 10^{-6}~{}\pm~{}1.54\times 10^{-7}$
C+(CEL) (1) | $1.23\times 10^{-4}~{}\pm~{}1.13\times 10^{-5}$
C2+(CEL) (1) | $5.31\times 10^{-4}~{}\pm~{}4.23\times 10^{-5}$
C2+(RL) (1) | $4.07\times 10^{-4}~{}\pm~{}8.71\times 10^{-5}$
C3+(CEL) (1) | $7.63\times 10^{-5}~{}\pm~{}7.04\times 10^{-6}$
N0 (1) | $4.81\times 10^{-7}~{}\pm~{}4.19\times 10^{-8}$
N+ (3) | $1.13\times 10^{-5}~{}\pm~{}2.29\times 10^{-7}$
O0 (1) | $8.78\times 10^{-6}~{}\pm~{}2.85\times 10^{-7}$
O+ (2) | $3.07\times 10^{-5}~{}\pm~{}2.84\times 10^{-6}$
O2+ (3) | $1.01\times 10^{-4}~{}\pm~{}1.67\times 10^{-6}$
O3+ (1) | $7.34\times 10^{-6}~{}\pm~{}3.52\times 10^{-6}$
Ne+ (1) | $4.07\times 10^{-7}~{}\pm~{}5.11\times 10^{-8}$
Ne2+ (1) | $1.91\times 10^{-6}~{}\pm~{}1.29\times 10^{-7}$
S+ (2) | $2.69\times 10^{-8}~{}\pm~{}1.36\times 10^{-9}$
S2+ (3) | $8.25\times 10^{-8}~{}\pm~{}6.96\times 10^{-9}$
S3+ (1) | $2.99\times 10^{-8}~{}\pm~{}4.71\times 10^{-9}$
Ar2+ (1) | $2.42\times 10^{-8}~{}\pm~{}2.81\times 10^{-9}$
Element X |  $n$(X)/$n$(H)
He | $1.04\times 10^{-1}~{}\pm~{}3.53\times 10^{-3}$
C(CEL) | $7.31\times 10^{-4}~{}\pm~{}4.44\times 10^{-5}$
N | $5.45\times 10^{-5}~{}\pm~{}5.46\times 10^{-6}$
O | $1.47\times 10^{-4}~{}\pm~{}4.83\times 10^{-6}$
Ne | $2.60\times 10^{-6}~{}\pm~{}1.89\times 10^{-7}$
S | $1.39\times 10^{-7}~{}\pm~{}8.51\times 10^{-9}$
Ar | $4.09\times 10^{-8}~{}\pm~{}6.38\times 10^{-9}$
Element X | 12 + $\log_{10}$ $n$(X)/$n$(H)
He | $11.02~{}\pm~{}0.01$
C(CEL) | $8.86~{}\pm~{}0.03$
N | $7.74~{}\pm~{}0.05$
O | $8.17~{}\pm~{}0.01$
Ne | $6.42~{}\pm~{}0.04$
S | $5.14~{}\pm~{}0.03$
Ar | $4.61~{}\pm~{}0.07$

We measured the line-fluxes of atomic species appeared in the 1-D KOOLS-IFU,
_IUE_ , and _Spitzer_ spectra using Alfa (Wesson, 2016). For several lines, we
manually measured their fluxes by fitting a multiple Gaussian function to
their line-profiles. In table A2, we summarize the resultant dereddened line-
fluxes normalized to the line-flux of H$\beta$ whose intensity is 100. Here,
the dereddened H$\beta$ line-flux $I$(H$\beta$) is ($4.54\pm 0.25$)$\times
10^{-13}$ erg s-1 cm-2.

To calculate elemental abundances, $T_{\rm e}$, $n_{\rm e}$, and then ionic
abundance derivations are required. The $n_{\rm e}$ and $T_{\rm e}$ values
using collisionally excited lines (CELs) are determined at the intersection of
the $n_{\rm e}$ and $T_{\rm e}$ diagnostic curves in each ionization zone
presented in table 3. In the neutral/low-ionization zone where the neutral or
singly ionized lines would be emitted, $T_{\rm e}$([N ii]) and $n_{\rm e}$([S
ii]) are determined from their curves.

In the high-ionization zone where the doubly or more ionized lines would be
emitted, $T_{\rm e}$([O iii]) and $n_{\rm e}$([S iii]) are derived from their
curves. $T_{\rm e}$([N ii] and $n_{\rm e}$([S ii]) are derived from the [N ii]
($I$(6548 Å)+$I$(6583 Å))/$I$(5755 Å) and [S ii] $I$(6717 Å)/$I$(6731 Å).
$T_{\rm e}$([O iii] and $n_{\rm e}$([S iii]) are derived from the [O iii]
($I$(4959 Å)+$I$(5007 Å))/$I$(4363 Å) and [S iii] $I$(18 )/$I$(9068 Å).
Because of no detection of the $[$N iii$]$ 1744-54 Å lines in the _IUE_
spectrum, we did not consider the recombination contribution of the N2+ to the
[N ii] 5755 Å line. The contributions of the O2+ and O3+ to the [O ii] 7320/30
Å and [O iii] 4363 Å lines are not subtracted because their contributions are
very small as reported by Otsuka & Tajitsu (2013). $T_{\rm e}$(He i) is
derived from the He i $I$(7281 Å)/$I$(6678 Å) ratio and the $T_{\rm
e}$–$n_{\rm e}$ logarithm interpolate function of the He i effective
recombination coefficients by Otsuka (2022). In the calculation, we adopt
$n_{\rm e}$([S ii]) and also subtract the He ii 6680 Å line-flux from the
complex of the He i 6678 Å plus He ii 6680 Å lines based on the theoretical
ratio of the He ii $I$(6680 Å)/$I$(4686 Å) under $T_{\rm e}$([O iii]) and
$n_{\rm e}$([S iii]).

The ionic abundance $n$(Xi+)/$n$(H+) and elemental abundances are summarized
in table 3. In the CEL ionic abundances, the values are obtained by solving
the equations of the populations at multiple energy levels under the obtained
CEL $T_{\rm e}$ and $n_{\rm e}$. When two or more lines are available, we
calculate the ionic abundance using each line. Then, we perform the weight-
average among their values. In calculations of Xi+/H+, $T_{\rm e}$ and $n_{\rm
e}$ of H+ and co-existing Xi+ are the same. $n$(He+)/$n$(H+) is calculated
using the $T_{\rm e}$(He i) and $n_{\rm e}$([S ii]) values. $n$(He2+)/$n$(H+)
is calculated using the $T_{\rm e}$([O iii]) and $n_{\rm e}$([S iii]) values.
$n$(C2+)/$n$(H+) from the recombination line C ii 7236 Å is calculated using
the $T_{\rm e}$([O iii]) and $n_{\rm e}$([S iii]) values. In the calculation
of the CEL ionic abundances, we adopt $T_{\rm e}$([N ii]) and $n_{\rm e}$([S
ii]) for the neutral or singly ionized species, and $T_{\rm e}$([O iii]) and
$n_{\rm e}$([S iii]) for the doubly or more ionized ones.

Finally, we calculate the elemental abundances $n$(X)/$n$(H) by introducing
the ionization correction factor ICF(X) to element X. ICFs recover the ionic
abundances in unobserved ionization stages covered by the obtained spectra.
The elemental abundance is calculated by ICF(X) $\cdot$ $\sum$
$n$(Xi+)/$n$(H+). We determine ICF(X) by referring to Otsuka & Tajitsu (2013).
For He, C, O, and S, their elemental abundances are derived by simply summing
up the derived ionic abundances. ICF(N) is equal to the O/O+ value ($4.82\pm
0.50$). ICF(Ne) corresponds to the O/(O+ \+ O2+) value ($1.13\pm 0.05$).
ICF(Ar) is the S/S2+ value ($1.70\pm 0.19$).

H4-1 is proven to be an extremely C-rich and metal-poor PN; the respective CEL
[C,N,O,Ne,S,Ar/H] are $+0.42\pm 0.05$, $-0.09\pm 0.07$, $-0.52\pm 0.05$,
$-1.50\pm 0.11$, $-1.97\pm 0.04$, and $-1.78\pm 0.15$, where Asplund et al.
(2009) is adopted as the solar abundance. The CEL C abundance (12 +
$\log_{10}n({\rm C})/n({\rm H})$) shows large scatter among previous works,
from 8.68 dex in Kwitter et al. (2003) to 9.02 dex in Otsuka & Tajitsu (2013).
This variation can be attributed to the difference in the slit-entrance size
adopted between the _IUE_ observations and the ground-based long-slit
observations, and the accuracy of aperture correction is compromised. As our
KOOLS-IFU study took the aperture correction into account, our derived CEL C
abundance is likely to be more reliable than previous estimates. We noticed a
large difference in the N abundance (12 + $\log_{10}n({\rm N})/n({\rm H})$) in
comparison to 7.56 dex, previously reported by Otsuka & Tajitsu (2013). This
is probably because the slit-length used in their Subaru/HDS observations was
not long enough to cover the neutral/singly ionized regions where [N ii] lines
are emitted.

##  6 Gas-to-Dust mass ratio, gas and dust masses

Table 4: Comparison of elemental abundances between the models by Lugaro et al. (2012) and our observation. The initial [$\alpha$/Fe] in the 2.0 M⊙ model are the same as the prediction by the Galaxy chemical evolution model by Kobayashi et al. (2011).  Model of initial mass | C | N | O | Ne | S | Ar
---|---|---|---|---|---|---
0.89 M⊙ (w/o PMZ) | $9.07$ | $7.59$ | $7.48$ | $7.33$ | $5.00$ | $4.28$
1.00 M⊙ (w/o PMZ) | $8.15$ | $6.56$ | $6.97$ | $7.33$ | $5.00$ | $4.26$
1.25 M⊙ (w/o PMZ) | $8.94$ | $6.66$ | $7.31$ | $6.95$ | $5.00$ | $4.27$
1.25 M⊙ (w/o PMZ) | $8.90$ | $6.68$ | $7.47$ | $7.65$ | $5.00$ | $4.27$
1.50 M⊙ (w/o PMZ) | $9.26$ | $6.80$ | $7.56$ | $7.68$ | $5.01$ | $4.29$
2.00 M⊙ (w/o PMZ) | $9.55$ | $6.84$ | $7.93$ | $8.66$ | $5.34$ | $4.57$
H4-1 | $8.86$ | $7.74$ | $8.17$ | $6.42$ | $5.14$ | $4.61$

The $c$(H$\beta$) value is the sum of the line-of-sight ISM and circumstellar
logarithmic reddening coefficient at 4861.33 Å. By subtracting the ISM
$c$(H$\beta$) value (0.01) from the derived $c$(H$\beta$) value (section 4),
we obtain the reddening coefficient due only to circumstellar media. Thus, we
derive the circumstellar gas-to-dust mass ratio (GDR). _Spitzer_ /IRS spectrum
exhibits (1) emission bands at $5-9$ , 11.3 , and 12.7 which would be
attributed to poly-aromatic amorphous carbon (PAH) molecules and (2) feature-
less dust continuum which would be attributed to amorphous carbon grains.
Here, we estimate the GDR $\psi$ and the gas/dust masses $m_{\rm g}$/$m_{\rm
d}$ using equations (7)–(13) established by Otsuka (2022).

When gas and dust grains coexist in the same volume, $\psi$ can be expressed
by equation (15):

$\displaystyle\psi$ $\displaystyle=$ $\displaystyle\frac{3(\log\,e)Q_{\rm
ext}({\rm H}\beta,a)C}{4a\rho_{\rm d}c({\rm H}\beta)_{\rm cir.}}\sum\frac{4\pi
I_{g}\mu_{g}}{n_{\rm e}j_{\rm g}(T_{\rm e},n_{\rm e})},$ (15)

where $Q_{\rm ext}$(H$\beta$,$a$) is the extinction cross section of the grain
of radius $a$ at H$\beta$, and $\rho_{\rm d}$ and $a$ are the density and the
radius of the grain, respectively. $C$ is a dimensionless constant: when
assuming the spherical shaped nebula with the radius of $r$ in arcsec, $C$
equals to $\pi^{-1}(4.848\times 10^{-6}~{}r)^{-2}$ (Otsuka, 2022).
$c$(H$\beta$)cir. is the line-of-sight reddening coefficient due to
circumstellar media, corresponding to $c$(H$\beta$)tot. (the sum of the line-
of-sight ISM and circumstellar reddening coefficient, i.e., $0.103\pm 0.024$)
minus $c$(H$\beta$)ISM (due to ISM, i.e., 0.01). $n_{\rm g}$ and $\mu_{\rm g}$
are the number density and atomic weight of the target gas, respectively.
$j_{\rm g}$($T_{\rm e}$,$n_{\rm e}$) is the volume emissivity of the gas
component per $n_{\rm g}n_{\rm e}$. $I_{\rm g}$ is the observed extinction-
free line flux of the gas component that is corrected using equation (7) with
$c$(H$\beta$)tot..

In calculating GDR, we utilize the H i, He i, and He ii lines. $T_{\rm e}$ =
$12634\pm 1412$ K and $n_{\rm e}$ = $420\pm 38$ cm-3 are adopted for the H+
mass, $T_{\rm e}$(He i) and $n_{\rm e}$([S ii]) for the He+ mass, $T_{\rm
e}$([O iii]) and $n_{\rm e}$([S iii]) for the He2+ mass. We adopt $r=1.5$
estimated from the H2 radial profile presented in Tajitsu et al. (2023). We
also adopt an AC-type (a-C:H) amorphous carbon grain ($\rho_{\rm d}$ of 1.85 g
cm-3) with the radius of $0.01-0.50$ and $a^{-3.5}$ size distribution (the
average $a$ of 0.07 ). $\langle Q_{\rm ext}({\rm
H}\beta,0.01-0.50\,{\micron})\rangle$ of 0.88 from the optical data of Rouleau
& Martin (1991). Then, we obtain GDR of $387\pm 86$ including the uncertainty
of line flux, c(H$\beta$), $T_{\rm e}$, and $n_{\rm e}$.

Next, we directly calculate the total atomic gas mass $m_{\rm g}$ of $(1.28\pm
0.24)\times 10^{-3}$ $D_{\rm kpc}^{2}$ M⊙ at the distance of $D_{\rm kpc}$ in
kpc, by simply adding the H+ and He+,2+ gas masses derived from equation (16):

$\displaystyle m_{\rm g}$ $\displaystyle=$ $\displaystyle\sum\frac{4\pi
D^{2}I_{g}\mu_{g}}{n_{\rm e}j_{\rm g}(T_{\rm e},n_{\rm e})}.$ (16)

The respective H and C gas masses are $(9.04\pm 1.70)\times 10^{-4}$ $D_{\rm
kpc}^{2}$ M⊙ and $(7.93\pm 1.57)\times 10^{-6}$ $D_{\rm kpc}^{2}$ M⊙. Adopting
the determined $m_{\rm g}$, we obtain the dust mass $m_{\rm d}$ of $(3.31\pm
0.62)\times 10^{-6}$ $D_{\rm kpc}^{2}$ M⊙. Note that the derived $m_{\rm g}$
and $m_{\rm d}$ here indicate the fully ionized gas mass and the dust mass co-
existing with it. The total mass of the C atom in the gas phase and that
incorporated into the dust grains is estimated to be $(1.12\pm 0.17)\times
10^{-5}$ $D_{\rm kpc}^{2}$ M⊙. This means $\sim$30 $\%$ of the total C atom
synthesized by the progenitor could exist as dust grains. As a reference, 12 +
$\log_{10}$ $n$(C)/$n$(H) would be $9.02\pm 0.11$ when including the C atom
locked in the dust grains. Through the photoionization model using the derived
GDR value ($387\pm 86$) as a constraint, we will estimate the gas mass
including the neutral hydrogen and molecules later.

##  7 Discussions

###  7.1 The possible origin and evolution of H4-1

The initial progenitor mass range can be determined by comparing the observed
elemental abundances with their predicted counterparts taken from AGB
nucleosynthesis models for a single star (Lugaro et al. 2012; table 4). Note
that the values in H4-1 are in the gas phase. As the available AGB
nucleosynthesis models for stars with $Z=10^{-4}$ are only those by Lugaro et
al. (2012), we compare their model results with the observed elemental
abundances.

Compared with the other initial mass models, the 0.89 M⊙ initial-mass model
without a partial mixing zone (PMZ222In low-mass AGB stars, a PMZ between the
bottom of the H-rich convective envelope and the outermost region of the
12C-rich intershell layer leads to a synthesis of extra 13C and 14N. Hence,
the presence of a PMZ would yield more C, O, Ne, and neutron-capture elements
than cases without a PMZ (e.g., Karakas, 2010) ) is very close to the observed
abundances. However, the discrepancy in the Ne and Xe333 For example, the 0.89
M⊙ initial-mass model predicts 12+$\log_{10}$(Xe/H) = 1.54. Otsuka & Tajitsu
(2013) conclude that the observed Xe abundance can be explained by models with
initially $r$-process element rich. abundances is rather large even for this
preferred single-star interpretation, as seen in table 4 and previously
determined ($12+\log_{10}$(Xe/H) $>2.75$; Otsuka & Tajitsu, 2013). We should
also note that the evolutionary time scale since the end of the AGB phase for
the 0.89 M⊙ initial-mass model appears much too long (table 5). Thus, in
practice, we need not consider the AGB nucleosynthesis model for stars of
initially 0.89 M⊙.

The 1.00 M⊙ initial-mass model without PMZ predicts much lower C and O
abundances (8.16 and 6.94 dex, respectively) compared to the observed values.
The 1.25 M⊙ initial-mass model without PMZ shows good agreement in terms of C
and Ne abundances, whereas O and Xe (0.45 dex) abundances are very low (table
4, w/o PMZ). On the other hand, the 1.25 M⊙ initial-mass model with PMZ (of
$2.0\times 10^{-3}$ M⊙; table 4, w/ PMZ) predicts more C, O, Ne, and Xe (1.45
dex). The 1.50 M⊙ initial-mass model without PMZ appears similar to the 0.89
M⊙ initial-mass model except for N and O. Previously, it was found that the
observed elemental abundances agreed model predictions from the 2.00 M⊙
initial-mass model with PMZ of $2.0\times 10^{-3}$ M⊙ and initially
$\alpha$-element/$r$-process-element rich ([$\alpha$/Fe]$\sim+0.3$ and
[$r$/Fe]$\sim+0.4$), except for C, N and Ne abundances (Otsuka & Tajitsu,
2013). Based on the reduced-$\chi^{2}$ evaluation between the observed and
predicted C, N, O, Ne, and Xe abundances, the 1.25 M⊙ initial-mass star model
without PMZ corresponds the most favorably to the observed abundances. Even if
adopting the C abundance including the C atom locked in the dust grains
($9.02\pm 0.11$ dex), this initial mass model keeps good fit to the observed
abundances. Hence, the progenitor mass of H4-1 could have been greater than
$0.8-0.9$ M⊙, but not so much more than $\sim 2$ M⊙.

However, all single star models are still largely discrepant in terms of the
N, O, and Ne abundances. The high N abundance in H4-1 could be caused by
extra-mixing during the red giant branch phase and also by helium-flash-driven
deep mixing during the AGB phase. The latter process would produce 14N through
the 13C($p$,$\gamma$)14N reaction by mixing protons into the He-burning shell.
Incorporating extra overshooting, which was not considered by Lugaro et al.
(2012), more C and O could be dredged up from the core to the stellar surface
(M. Lugaro, priv. comm.). However, this mechanism would also cause an Ne
overabundance. To lower the Ne abundance, the progenitor star could start as
Ne poor or extra mixing would have to be invoked (22Ne abundance would be
reduced; M. Lugaro, priv. comm). Henry (1989) interpreted that H4-1 formed in
the low Ne environment. We believed that the low Ne abundance can be improved
by adding the Ne+ abundance re-calculated using the PSF-matched KOOLS-IFU and
Spitzer spectra. However, the Ne abundance remains low and also cannot be
reproduced by single star AGB nucleosynthesis models. The extremely low Ne
abundance in H4-1 seems to be a serious problem in the stellar nucleosynthesis
and the Galaxy chemical evolution (Henry et al., 2004). Because there does not
seem to be any single star models that simultaneously satisfy all of the
observed N, O, and Ne abundances, other possibilities, e.g., binary star
evolution, would have to be considered.

Otsuka & Tajitsu (2013) already speculated that a binary system could explain
the rather anomalous elemental abundances seen from H4-1. The apparently
anomalous N, O, and Ne abundances of H4-1 could originate from its binarity.
According to Miller Bertolami (2016), a single star of $0.85-2.5$ M⊙ initial
mass is predicted to eventually evolve into a white dwarf of $\sim 0.53-0.75$
M⊙. This naively means that a binary system composed of a 2.5 M⊙ primary and a
0.85 M⊙ companion would become a star of about 1.6 M⊙ via a merger event that
takes place after the primary experienced mass loss. Such a merged system
could exhibit abundances that are largely consistent with originally $0.8-2.5$
M⊙ single stars (table 4).

To further determine the origin of H4-1, we quantify the total circumstellar
mass based on the observational data. To that end, we first estimate the
distance to H4-1 by assuming (i) a single 0.85 M⊙ initial-mass star and (ii) a
binary system in which a 0.85 M⊙ main-sequence secondary merging with a white
dwarf primary that has become a single main-sequence star of 1.25 M⊙ or a 1.5
M⊙.

###  7.2 Distance toward H4-1

The distance to H4-1 has not yet been well constrained, from 4.4 to 25.3 kpc
with the average of 15.68 kpc (1-sigma = 6.1 kpc) in the literature (Cahn &
Kaler, 1971; Acker, 1978; Kohoutek & Martin, 1981; Maciel, 1984; Cahn et al.,
1992; Kingsburgh & Barlow, 1992; Zhang, 1995; Mal’Kov, 1997; Stanghellini et
al., 2008; Frew et al., 2016). Assuming that the nebula is optically thick and
absorbs all the energy emitted from the central star and re-emits it, the
integrated flux from the central star corresponds to the sum of all the
observed fluxes of the nebular atomic and molecular gas and dust grains and is
proportional to the square of the distance. Therefore, even if no direct
spectrum and photometry of the central star itself is available, the distance
can still be estimated by comparing all the nebular-line/continuum flux of a
target PN with those of a reference PN whose distance is well known.

We note, however, that this distance estimate is valid only when the target
and reference PNe have similar initial mass/metal abundances and their central
stars are in the same evolutionary stage. As far as we are aware, K648 in M15
is the only candidate as a reference PN for H4-1: both are halo PNe and the
distance to K648 is well calibrated (Otsuka et al., 2015). However, K648 has a
cooler central star ($T_{\rm eff}=36360$ K), still evolving toward the PN
turn-off, while H4-1 has a hotter central star (e.g., $T_{\rm
eff}=66000-99000$ K; Otsuka et al., 2003) already in the white dwarf cooling
track. Therefore, the distance to H4-1 obtained through this flux comparison
method remains uncertain at best.

Table 5: The estimated distance ($D$) toward H4-1 for different initial masses, effective temperatures ($T_{\rm eff}$), luminosities ($L_{\ast}$), surface gravities ($\log_{10}\,g$) of the central star. We define the age as the time since when $T_{\rm eff}=5000$ K. Dagger ($\dagger$; the first row for each initial mass star) indicates that the central star is still becoming hot, while the rest are later epochs on the cooling track.  Init. mass | $T_{\rm eff}$ | $L_{\ast}$ | $\log_{10}\,g$ | $D$ | Age
---|---|---|---|---|---
(M⊙) | (K) | ($L_{\odot}$) | (cm s-2) | (kpc) | ($\times 10^{4}$ yrs)
$0.85^{\dagger}$ | 125000 | $1179.1$ | 6.43 | $25.03$ | 3.66
$0.85$ | 125000 | $406.7$ | 6.89 | $14.70$ | 5.72
$0.85$ | 122500 | $332.1$ | 6.96 | $13.28$ | 6.23
$0.85$ | 120000 | $274.3$ | 7.00 | $12.07$ | 6.59
$1.25^{\dagger}$ | 125000 | $7139.9$ | 5.70 | $62.21$ | 0.34
$1.25$ | 125000 | $174.4$ | 7.31 | $9.63$ | 0.49
$1.25$ | 122500 | $156.8$ | 7.32 | $9.13$ | 0.51
$1.25$ | 120000 | $140.4$ | 7.33 | $8.64$ | 0.54
$1.50^{\dagger}$ | 125000 | $8758.5$ | 5.61 | $68.22$ | 0.19
$1.50$ | 125000 | $136.9$ | 7.41 | $8.53$ | 0.46
$1.50$ | 122500 | $123.4$ | 7.42 | $8.10$ | 0.54
$1.50$ | 120000 | $111.1$ | 7.43 | $7.68$ | 0.66

There is an alternative method to determine the distance to H4-1 based on the
adopted spectro-photometry data. We estimate the distance by comparing the
luminosity derived from the sum of the atomic gas/molecular H2 line emission
and dust continuum emission with the theoretical post-AGB evolution models by
Miller Bertolami (2016). Detailed explanations as to how to derive the
distance are found in Appendix A. Here, we consider 0.85, 1.25, and 1.5 M⊙
initial-mass stars following the discussion given in the previous section.
Based on the evolutionary tracks by Miller Bertolami (2016), 0.85, 1.25, and
1.5 M⊙ initial-mass stars are expected to become white dwarfs with the final
core-mass of 0.53, 0.58, and 0.60 M⊙, respectively. Hence, the total ejected
mass during the evolution of the 0.85, 1.25, and 1.5 M⊙ progenitor star is
0.32, 0.67, and 0.90 M⊙, respectively.

In optically thick conditions, the distance ($D$) to a PN can be recovered if
the luminosity of the central star ($L_{\ast}$) is known, as the total
luminosity should equal to the integrated flux – the sum of nebular continuum
and lines from the ionized and atomic gas, molecular lines (i.e., H2 and
PAHs), and dust continuum – multiplied by $4{\pi}D^{2}$. Under such optically
thick conditions, we can simultaneously fit $D$, $L_{\ast}$, and the surface
gravity ($\log_{10}\,g$) of the central star for different initial masses at
different $T_{\rm eff}$ (representing different evolutionary epochs) as
summarized in table 5. With this method, the uncertainty of $D$ comes out to
be $\sim 10$ $\%$ at 95 $\%$ confidence.

For each star of specific initial mass in table 5, the first row corresponds
to an evolutionary epoch while the central star is still becoming hot prior to
the PN turn-off, whereas the following three rows correspond to the subsequent
epochs along the cooling track after the PN turn-off. While the central star
is still getting hotter as it evolves toward the PN turn-off, the ejected gas
mass can be estimated using the relation $m_{\rm g}$ = $(1.28\pm 0.24)\times
10^{-3}$ $D_{\rm kpc}^{2}$ M⊙ established in section 6: a star of 0.85 M⊙,
1.25 M⊙, and 1.5 M⊙ initial-mass (at $D=25.03$, $62.21$, and $68.22$ kpc)
would have to have ejected the circumstellar mass of $\geq 0.83$ M⊙, $\geq
5.15$ M⊙, and $\geq 6.19$ M⊙, respectively. Since these ejecta masses are
obviously too large, H4-1 must be a PN presently on the white dwarf cooling
track. If so, H4-1 is unlikely to be a $\sim 0.8$ M⊙ initial-mass star unless
its evolution were accelerated. Therefore, the distance toward H4-1, presently
being a $1.25-1.5$ M⊙ white dwarf on the cooling track, would be $\sim 10$ kpc
or less (table 5).

###  7.3 The circumstellar gas mass

Now we want to consider the circumstellar mass. If the detected gas mass is
more than $0.32\,\mathrm{M}_{\odot}$, the progenitor would be most likely a
binary system. Thus, we need to verify the binary nature of H4-1 by estimating
the amount of the circumstellar mass-loss ejecta. To do so, we construct
photoionization models using Cloudy (v17.03; Ferland et al. 2017) by adopting
the distances derived above.

Table 6:  The adopted and derived parameters in the Cloudy models. GDR1 value corresponds to the ratio of the total gas mass to the AC mass. GDR2 value corresponds to the ratio of the total gas mass to the (AC+PAH) mass.  Central star | Values for 1.25 M⊙ model | Values for 1.50 M⊙ model
---|---|---
$T_{\rm eff}$ / $\log_{10}\,g$ / $L_{\ast}$ / $D$ | 127452 K/7.30 cm s-2/ 192 $L_{\odot}$/ 10.117 kpc | 128428 K/7.40 cm s-2/160 $L_{\odot}$/9.225 kpc
Nebula | Values for 1.25 M⊙ model | Values for 1.50 M⊙ model
Elemental abundances | He: 11.04, C: 8.79, N: 7.52, O: 8.20, | He: 11.03, C: 8.79, N: 7.50, O: 8.20,
(12 + $\log_{10}$$n$(X)/$n$(H)) | Ne: 6.28, S: 5.07, Ar: 4.44 | Ne: 6.25, S: 5.08, Ar: 4.46
| Others: Lugaro et al. (2012) | Others: Lugaro et al. (2012)
Geometry | spherical shell nebula | spherical shell nebula
| inner radius = 2399 AU (0.24) | inner radius = 2359 AU (0.20)
| outer radius = 23356 AU (2.31) | outer radius = 30313 AU (2.57)
Ionization boundary radius | 13252 AU (1.31) | 12038 AU (1.31)
Radial hydrogen density profile | 1297 cm-3 in $r\leq 13252$ AU | 1062 cm-3 in $r\leq 12038$ AU
($n_{\rm H}(r)$) | 2722 cm-3 in $r>13252$ AU | 2037 cm-3 in $r>12038$ AU
Deredden H$\beta$ line flux ($I$(H$\beta$)) |  $4.54\times 10^{-13}$ erg s-1 cm-2 |  $4.54\times 10^{-13}$ erg s-1 cm-2
Temperature floor | 727 K | 749 K
H2 column density ($\log_{10}N({\rm H_{2}})$ | 18.69 cm-2 | 18.61 cm-2
Ionized gas mass | 0.194 M⊙ | 0.169 M⊙
Neutral gas mass | 0.186 M⊙ | 0.208 M⊙
Total gas mass | 0.380 M⊙ | 0.377 M⊙
AC and PAHs | Values for 1.25 M⊙ model | Values for 1.50 M⊙ model
AC temp. |  $25-122$ K |  $24-131$ K
PAH temp. |  $83-198$ K |  $80-210$ K
AC mass |  $9.31{\times}10^{-4}$ M⊙ |  $8.78{\times}10^{-4}$ M⊙
Ionized PAH mass |  $3.37{\times}10^{-5}$ M⊙ |  $3.12{\times}10^{-5}$ M⊙
Neutral PAH mass |  $1.82{\times}10^{-5}$ M⊙ |  $1.95{\times}10^{-5}$ M⊙
AC + PAH masses |  $9.83{\times}10^{-4}$ M⊙ |  $9.29{\times}10^{-4}$ M⊙
GDR1 | 405 | 429
GDR2 | 386 | 406
Reduced-$\chi^{2}$ | 15.1 | 17.8

The central star is supposedly on the white dwarf cooling track, and its large
surface gravity suggests that heavy elements remain deep in the photosphere.
This means that the surface stellar abundances would be composed mainly of
light elements. Hence, as the input SED, we adopt theoretical spectra of non
local thermodynamic equilibrium (non-LTE) stellar atmosphere models by Rauch
(2003) for stars whose atmosphere is composed of He and H (the He/H mass ratio
of 0.35).

Metals are significant coolants in ionized nebula. With Cloudy, we adopt
elements up to Zinc as coolants. Starting from the observed elemental
abundances (table 3) as the input values, we converge on the best-fit model by
iteratively adjusting elemental abundances. To constrain carbon abundance, we
use the CEL C lines. The abundances of the undetected elements lighter than Zn
are fixed at the values predicted for models of 1.25 M⊙ and 1.5 M⊙ initial
masses (Lugaro et al., 2012).

To estimate the amount of dust, we consider spherical amorphous carbon
particles of $0.01-0.50$ radius ($a$) following the power-law size
distribution (of $a^{-3.5}$). The best-fit grain size range is determined
based on a grid of models. The broad PAH band at $5-9$ and two emission bands
at 11.3/12.7 are attributed to charged and neutral species, respectively.
Hence, we adopt both charged and neutral PAH molecules with the radius of
$(3.6-11)\times 10^{-4}$ that follow the same $a^{-3.5}$ size distribution. In
computing PAH emission, we adopt the optical constants by Draine & Li (2007)
and use the stochastic heating mechanism.

We aim at reproducing pure-rotational H2 lines at 9.67, 12.27, and 17.03 as
well as all the detected atomic lines. As discussed by Otsuka et al. (2017)
and Otsuka & Hyung (2020), the presence of H2 lines in the mid-IR requires
high hydrogen densities in a warm temperature photodissociation region (PDR).
Thus, we adopt a two-layer spherical nebula with distinct hydrogen density
radial profiles ($n_{\rm H}(r)$, where $r$ is the distance from the central
ionization source) to be fit through a grid of models. We set the lower
temperature limit in PDRs using the temperature floor option. We can regard
the value of the temperature floor as the excitation temperature for these
three H2 lines. We adopt the UMIST chemical reaction network by McElroy et al.
(2013).

Figure 4: (left two panels) Comparison between the observed spectra/photometry
and the predicted SED for the 1.25 M⊙ star. The lower panel focuses on the
$3-40$ range, where the spectral resolution of the model SED is matched with
the Spitzer/IRS spectrum. (right panel) The current location of the central
star plotted on the theoretical post-AGB evolutionary tracks by Miller
Bertolami (2016).

We vary $T_{\rm eff}$ between 120000 and 130000 K. At each $T_{\rm eff}$, we
calculate $L_{\ast}$, $\log_{10}\,g$, and $D$ by the method described in
subsection 7.2, and then, keep $T_{\rm eff}$, $L_{\ast}$, $\log_{10}\,g$, and
$D$. We emphasize that $T_{\rm eff}$, $L_{\ast}$, $\log_{10}\,g$ values for
each initial mass star are completely along to those predicted in the post-AGB
evolutionary tracks of Miller Bertolami (2016) and also $L_{\ast}$,
$\log_{10}\,g$, and $D$ are the function of $T_{\rm eff}$. We vary the
following 12 parameters: 7 elemental abundances (He/C/N/O/Ne/S/Ar), the inner
radius of the nebula, temperature floor in PDRs, GDR determined from
neutral/ionized PAHs and amorphous carbon mass fractions. We vary the
temperature floor in order to reproduce the selected H2 lines. Iterative
calculations are terminated either when the predicted flux density at 100
reaches the observed flux in the _Herschel_ /PACS 100 band or when the
predicted total gas mass amounts to 0.9 M⊙. The quality of fit is based on the
reduced-$\chi^{2}$ calculated from the following 78 constraints: 47 atomic and
3 H2 line fluxes relative to $I$(H$\beta$), $\log_{10}$ $I$(H$\beta$), 3 PAH
band fluxes, 23 mid- and far-IR flux densities, and the ionization boundary
radius of 1.31444This ionization boundary radius corresponds to 1.3 times
Gaussian FWHM of the H$\beta$ emission image, within which the 85 $\%$ of the
total flux is accounted for., which is defined as the radial distance from the
central star at which $T_{\rm e}$ drops below 4000 K.

Table 6 summarizes all the input and derived/best-fit quantities for each of
the two distinct initial mass models. The distance in the best-fit for 1.25
and 1.50 M⊙ stars is 10.12 and 9.23 kpc, respectively. The goodness of fit
according to the reduced-$\chi^{2}$ ($>15$) is dictated by the two [O ii]
lines. Excluding these lines, the reduced-$\chi^{2}$ value diminishes to $\sim
7$, which is comparable to other similarly detailed Cloudy model results for
two PNe, NGC650-1 (van Hoof et al., 2013) and NGC6781 (Otsuka et al., 2017).
Hence, we can say that both 1.25 and 1.5 M⊙ models equally reproduce the
observational results. Table A3 presents a full comparison between the
observed and modeled values, including atomic and molecular hydrogen line
fluxes, PAH band fluxes at $6-9$ and 11 , fluxes for two bands near the peak
dust emission intensity wavelength, and the dust emission flux density at 23
wavelengths. Figure 4 shows comparisons of the SED between the observed data
and 1.25 M⊙ initial mass model results (left panels) as well as where each of
the two models is found on the evolutionary tracks (Miller Bertolami 2016;
right panel) Both models successfully reproduce the observed SED across all
wavelengths. Also, the observed GDR ($387\pm 86$) is reproduced fairly well in
the models (386 and 406 in table 6). Hence, we conclude that the models
represent H4-1 fairly well. We emphasize here that both the ionized and
neutral gas components are separately accounted for in the models, residing
inside and outside of the ionization boundary, respectively. We immediately
see that H4-1 is surrounded by as much neutral gas as ionized gas. Hence, one
must miss the whole picture of a PN if only the ionized gas component is
focused on.

Both models suggest the total gas mass of $0.38\pm 0.08$ M⊙ at 95 %
confidence, which exceeds the upper and lower limit of the circumstellar gas
mass for the single star and binary system scenario, respectively (subsection
7.3). Given the sensitivity of observations and how the circumstellar shell
surface brightness drops off as a function of radial distance from the central
star, combined with our previous experience in modeling PNe, the neutral gas
mass estimated here accounts for probably $\sim 30-50\,\%$ of the existing
neutral gas mass (e.g., Otsuka et al., 2017; Otsuka & Hyung, 2020). This means
that there is a sensitivity limit beyond which tracing the mass-loss history
in the earlier AGB phase is challenging. Based on these quantities, the total
gas mass is scaled to be $0.57-0.81\,\mathrm{M}_{\odot}$ and
$0.59-0.86\,\mathrm{M}_{\odot}$ for the 1.25 M⊙ and 1.50 M⊙ initial-mass
models, respectively. Accordingly, the models suggest that the initial mass of
the progenitor (i.e., the sum of the present mass of the central star and the
total circumstellar gas mass) is
$0.58+(0.57\mathrm{~{}to~{}}0.81)=1.15\mathrm{~{}to~{}}1.39$ M⊙ and
$0.60+(0.59\mathrm{~{}to~{}}0.86)=1.19\mathrm{~{}to~{}}1.46$ M⊙ for the 1.25
M⊙ and 1.50 M⊙ initial-mass models, respectively, resulting in very self-
consistent. From the argument made above, we can confidently conclude that the
binary system scenario is more favorable for H4-1.

###  7.4 The evolutionary history of H4-1 as a binary

Provided that the progenitor binary system of H4-1 now evolves as a single
$\sim 1.25-1.50$ M⊙ after a merger event between a low-mass main-sequence
secondary of $\sim 0.8$ M⊙ initial mass and an evolved primary, we would like
to understand the initial mass of the primary of the progenitor binary system.
We also expect that binary models can explain the observed N, O, and Ne
abundances unexplained by AGB nucleosynthesis models for single stars.

Table 7: Comparison between the elemental abundances predicted by the binary nucleosynthesis models and by our observation. The second row indicates the initial abundances of the primary and secondary stars, which are consistent with the values predicted for $Z=10^{-4}$ by the Galaxy chemical evolution models by Kobayashi et al. (2011).  | He | C | N | O | Ne | S | Ar
---|---|---|---|---|---|---|---
Binary model | 11.04 | 9.03 | 7.47 | 7.87 | 6.54 | 5.35 | 4.59
initial abund. | 10.90 | 6.06 | 4.64 | 6.82 | 6.06 | 5.26 | 4.49
H4-1 | 11.02 | 8.86 | 7.74 | 8.19 | 6.42 | 5.14 | 4.61

For that purpose, we run a large grid of binary models using the binary-c code
(Izzard et al., 2004, ver 2.2.3), for which the parameters varied are the
initial mass of the primary ($M1$) and the secondary ($M2$) and the initial
separation ($R_{\rm sep}$). The adopted mass range for $M1$ and $M2$ is
$1.10-1.89$ M⊙ and $0.80-0.90$ M⊙, respectively, with a constant step of 0.01
M⊙. The range of $R_{\rm sep}$ is varied from 150 $R_{\odot}$ to 500
$R_{\odot}$ with a constant step of 1 $R_{\odot}$. To be consistent with the
Galaxy chemical evolution model by Kobayashi et al. (2011), we adopt
$\mathrm{[}\alpha\mathrm{/Fe]}=+0.3$. The PMZ mass associated with the He-
intershell products is set to be 0, which is in line with the discussion above
in subsection 7.1. Overshooting is not considered. We adopt the Reimer’s mass-
loss law with the fitting parameter $\eta=0.5$ during the red giant stage, and
specific rates of mass loss during the early AGB phase (we set GDR = 400 from
our study; Goldman et al., 2017) and the thermal pulsing AGB phase
(Vassiliadis & Wood, 1993). Note that the models by Lugaro et al. (2012)
referenced earlier adopt specific rates of mass loss (Vassiliadis & Wood,
1993) throughout the AGB phase irrespective of thermal pulses. Adopting high
mass loss during the early AGB phase (Goldman et al., 2017) can suppress the
unwanted Ne enhancement, compared with an exclusive use of the rates of mass
loss during the thermal pulsing AGB phase (Vassiliadis & Wood, 1993). This is
because that the mass loss in the early AGB determines the mass in the
beginning of TP phase and affects the duration of the TP phase. We set the
orbital eccentricity of 0, common envelope efficiency ($\alpha$) of 0.5, and
structure parameter ($\lambda$) of 0.5. $\alpha$ is the efficiency of
converting the orbital energy into the kinetic energy of the envelope and
$\lambda$ depends on the structure of the primary star, conventionally chosen
as a constant, typically $\sim 0.5$ (Zuo & Li, 2014). $\alpha$ and $\lambda$
control the ratio of final to initial orbital separation before and after
common envelope phase. The remaining parameters required to simulate the
binary evolution are set to the default values in the Binary-c code. Out of
the total number of models that turns out to be 308880, favorable candidates
are selected as those satisfying the following criteria; (1) the total
duration of evolution exceeds 10 Gyr, (2) the final mass of the primary is
greater than 0.58 M⊙, and (3) the C/O ratio (by number) is $>1$. Among these
favorable candidates, we find out the best-fit model for the progenitor system
of H4-1 by evaluating the reduced-$\chi^{2}$ between the observed and
predicted He/C/N/O/Ne/S/Ar abundances, assuming that the uncertainty of each
elemental abundance is of 0.1 dex.

In table 7, we compare the best-fit model abundances (the top row) with the
observed abundances (the bottom row), while the initial abundances adopted for
the models are shown in the middle. Adopting the observed C abundance
including the C atom incorporated into the dust grains ($9.02\pm 0.11$ dex),
the binary model gives better fitting. Either way, the binary model fairly
reproduces all the observed abundances. The best-fit model (the
reduced-$\chi^{2}$ of 1.58) suggests the initial mass of 1.87 M⊙ and 0.82 M⊙
for the primary and secondary, respectively, with the initial separation of
248 $R_{\odot}$. Presently, the now-merged star is of 0.62 M⊙ on the white
dwarf cooling track at the age of 10.02 Gyr since the formation of the
progenitor binary system in the Galactic halo.

The evolutionary history of H4-1, according to the best-fit model, is as
follows. About 3.8 Gyrs after the Big Bang, a Pop II binary system was formed
in the young Milky Way. When the primary of 1.87 M⊙ initial mass became an
early AGB star, the primary began to lose its envelope to the Roche-lobe
overflow, rapidly evolving to be a CO white dwarf of $\sim 0.6$ M⊙ without
experiencing TPs and TDUs. As the secondary went into the giant phases, it
experienced the Roche-lobe overflow several times while never losing its outer
envelope. The surface of the secondary was contaminated by the primary, its
He/C/N/O/Ne/S/Ar composition was changed into
10.97/7.65/7.43/7.77/6.33/5.33/4.57 dex, respectively. The N-excess observed
in H4-1 can be explained by the result of mass-transfer from the primary to
the secondary stars. After about 10 Gyrs since the formation, the back and
forth large-scale mass transfer between the white dwarf primary and the main
sequence secondary via the Roche-lobe overflow was repeated several times,
while the binary separation was reduced each time. Ultimately, the system
became a contact binary and then coalesced into a single star of 1.42 M⊙,
showing the surface composition identical to the sum of those from both stars.
The merged star subsequently underwent the TPs and TDUs required to become
C-rich during the AGB phase of its own and became a CO white dwarf of 0.62 M⊙
surrounded by a C- and N-rich and Ne-poor dusty nebulae. The multipolar
structure of the nebula and the presence of the equatorial torus are
consistent with the binarity of the progenitor system.

If adopting mass loss by Vassiliadis & Wood (1993) through the AGB phase, the
best-fit model (the 1.65 M⊙ primary and 0.86 M⊙ secondary stars with the
initial separation of 274 $R_{\odot}$) predicts the binary system ended as a
0.64 M⊙ CO white dwarf after binary mass-transfer and merger, but high C and
Ne abundances555The model predicts that He/C/N/O/Ne/S/Ar composition is
11.00/9.24/7.45/7.88/7.23/5.34/4.58 dex, respectively.. The low Ne in H4-1
would be due to the high mass loss during the early AGB.

As we discussed above, binary evolution completely resolved the long-standing
issues on evolutionary time scale and anormal chemical abundance pattern.
Therefore, we confidently conclude that H4-1 originated from a Pop II binary
star.

##  8 Summary

We have performed a comprehensive study of the Galactic halo PN H4-1 based on
the newly acquired Seimei/KOOLS-IFU data and auxiliary multiband spectral
photometry data to reveal its origin and evolution through accurate
measurements of gas abundances and mass. The spatially-resolved emission line
maps reconstructed from the KOOLS-IFU data cube exhibit both the elliptical
elongation along the polar axis and the central flattening along the
equatorial plane commonly seen as a signature of the presence of an equatorial
disk structure in bipolar PNe evolved from high-mass progenitors.

Fully data-driven analyses of the spectra that are point-spread-function-
matched at each wavelength bin have yielded seven elemental abundances, gas-
to-dust mass ratio, and gas/dust masses, which are also supported by the
distance to the object newly estimated consistently by ourselves. The detailed
photoionization model has reproduced all the observed quantities of the
adopted multiband spectral photometry data. The derived total gas mass
strongly indicates that the progenitor is not a typical single halo star but a
binary system. By modeling the possible binary evolution and two-star
nucleosynthesis, we have been able to reproduce the observed elemental
abundances that cannot be explained by nucleosynthesis of a single star. We
have also successfully estimated the initial mass and the initial orbital
separation of the progenitor binary system. The best-fit model predicts that
the progenitors of 1.87 M⊙ and 0.82 M⊙ initial mass were originally Pop II
stars formed about 4 Gyrs after the Big Bang, which underwent episodes of
mass-transfers and an eventual binary merger, subsequently evolved into a PN
showing unique chemical abundances. The N-excess would be due to a result of
mass-transfer, and the extremely low Ne would be due to the high mass loss
during the early AGB. Mutual follow-up between halo PNe and CEMP stars will
greatly help in understanding their origin and evolution.

We are grateful to the referee for a careful reading and valuable suggestions.
MO thanks Dr. Maria Lugaro for fruitful discussion on chemical abundances in
metal-poor stars. MO thanks Drs. Robert Izzard and David Hendriks to kindly
instruct him about the use of their developed codes. MO was supported by the
Japan Society for the Promotion of Science (JSPS) Grants-in-Aids for
Scientific Research (C) (19K03914 and 22K03675). MO thanks his parents for
their deep love and support. TU was partially supported by the JSPS through
its invitation fellowship program (FY2020, long-term).

## References

  * Acker (1978) Acker, A. 1978, A&AS, 33, 367
  * Alam et al. (2015) Alam, S., Albareti, F. D., Allende Prieto, C., et al. 2015, ApJS, 219, 12
  * Aparicio Villegas et al. (2010) Aparicio Villegas, T., Alfaro, E. J., Cabrera-Caño, J., et al. 2010, AJ, 139, 1242
  * Asplund et al. (2009) Asplund, M., Grevesse, N., Sauval, A. J., & Scott, P. 2009, ARA&A, 47, 481
  * Beers & Christlieb (2005) Beers, T. C., & Christlieb, N. 2005, ARA&A, 43, 531
  * Bianchi et al. (2017) Bianchi, L., Shiao, B., & Thilker, D. 2017, ApJS, 230, 24
  * Cahn & Kaler (1971) Cahn, J. H., & Kaler, J. B. 1971, ApJS, 22, 319
  * Cahn et al. (1992) Cahn, J. H., Kaler, J. B., & Stanghellini, L. 1992, A&AS, 94, 399
  * Cardelli et al. (1989) Cardelli, J. A., Clayton, G. C., & Mathis, J. S. 1989, ApJ, 345, 245
  * Currie et al. (2014) Currie, M. J., Berry, D. S., Jenness, T., et al. 2014, in Astronomical Society of the Pacific Conference Series, Vol. 485, Astronomical Data Analysis Software and Systems XXIII, ed. N. Manset & P. Forshay, 391
  * Cutri et al. (2013) Cutri, R. M., Wright, E. L., Conrow, T., et al. 2013, Explanatory Supplement to the AllWISE Data Release Products, Explanatory Supplement to the AllWISE Data Release Products, by R. M. Cutri et al.
  * Dorsch et al. (2019) Dorsch, M., Latour, M., & Heber, U. 2019, A&A, 630, A130
  * Draine & Li (2007) Draine, B. T., & Li, A. 2007, ApJ, 657, 810
  * Fahlman et al. (1985) Fahlman, G. G., Richer, H. B., & Vandenberg, D. A. 1985, ApJS, 58, 225
  * Fazio et al. (2004) Fazio, G. G., Hora, J. L., Allen, L. E., et al. 2004, ApJS, 154, 10
  * Ferland et al. (2017) Ferland, G. J., Chatzikos, M., Guzmán, F., et al. 2017, Rev. Mexicana Astron. Astrofis., 53, 385
  * Fluks et al. (1994) Fluks, M. A., Plez, B., The, P. S., et al. 1994, A&AS, 105, 311
  * Frew et al. (2016) Frew, D. J., Parker, Q. A., & Bojičić, I. S. 2016, MNRAS, 455, 1459
  * Geier et al. (2019) Geier, S., Raddi, R., Gentile Fusillo, N. P., & Marsh, T. R. 2019, A&A, 621, A38
  * Goldman et al. (2017) Goldman, S. R., van Loon, J. T., Zijlstra, A. A., et al. 2017, MNRAS, 465, 403
  * González-Santamaría et al. (2021) González-Santamaría, I., Manteiga, M., Manchado, A., et al. 2021, A&A, 656, A51
  * Henry (1989) Henry, R. B. C. 1989, MNRAS, 241, 453
  * Henry et al. (2004) Henry, R. B. C., Kwitter, K. B., & Balick, B. 2004, AJ, 127, 2284
  * Herschel Point Source Catalogue Working Group et al. (2020) Herschel Point Source Catalogue Working Group, Marton, G., Calzoletti, L., et al. 2020, VizieR Online Data Catalog, VIII/106
  * Houck et al. (2004) Houck, J. R., Roellig, T. L., van Cleve, J., et al. 2004, ApJS, 154, 18
  * Hubeny (1988) Hubeny, I. 1988, Computer Physics Communications, 52, 103
  * Izzard et al. (2004) Izzard, R. G., Tout, C. A., Karakas, A. I., & Pols, O. R. 2004, MNRAS, 350, 407
  * Jacoby et al. (1997) Jacoby, G. H., Morse, J. A., Fullton, L. K., Kwitter, K. B., & Henry, R. B. C. 1997, AJ, 114, 2611
  * Karakas (2010) Karakas, A. I. 2010, MNRAS, 403, 1413
  * Kingsburgh & Barlow (1992) Kingsburgh, R. L., & Barlow, M. J. 1992, MNRAS, 257, 317
  * Kobayashi et al. (2011) Kobayashi, C., Karakas, A. I., & Umeda, H. 2011, MNRAS, 414, 3231
  * Kohoutek & Martin (1981) Kohoutek, L., & Martin, W. 1981, A&A, 94, 365
  * Kurita et al. (2020) Kurita, M., Kino, M., Iwamuro, F., et al. 2020, PASJ, 72, 48
  * Kwitter et al. (2003) Kwitter, K. B., Henry, R. B. C., & Milingo, J. B. 2003, PASP, 115, 80
  * Kwok (2000) Kwok, S. 2000, The Origin and Evolution of Planetary Nebulae (Cambridge University Press)
  * Lawrence et al. (2007) Lawrence, A., Warren, S. J., Almaini, O., et al. 2007, MNRAS, 379, 1599
  * Lugaro et al. (2012) Lugaro, M., Karakas, A. I., Stancliffe, R. J., & Rijs, C. 2012, ApJ, 747, 2
  * Maciel (1984) Maciel, W. J. 1984, A&AS, 55, 253
  * Mal’Kov (1997) Mal’Kov, Y. F. 1997, Astronomy Reports, 41, 760
  * Matsubayashi et al. (2019) Matsubayashi, K., Ohta, K., Iwamuro, F., et al. 2019, PASJ, 71, 102
  * McElroy et al. (2013) McElroy, D., Walsh, C., Markwick, A. J., et al. 2013, A&A, 550, A36
  * Miller Bertolami (2016) Miller Bertolami, M. M. 2016, A&A, 588, A25
  * Otsuka (2019) Otsuka, M. 2019, MNRAS, 482, 2354
  * Otsuka (2022) —. 2022, MNRAS, 511, 4774
  * Otsuka & Hyung (2020) Otsuka, M., & Hyung, S. 2020, MNRAS, 491, 2959
  * Otsuka et al. (2009) Otsuka, M., Hyung, S., Lee, S.-J., Izumiura, H., & Tajitsu, A. 2009, ApJ, 705, 509
  * Otsuka et al. (2015) Otsuka, M., Hyung, S., & Tajitsu, A. 2015, ApJS, 217, 22
  * Otsuka et al. (2008) Otsuka, M., Izumiura, H., Tajitsu, A., & Hyung, S. 2008, ApJ, 682, L105
  * Otsuka & Tajitsu (2013) Otsuka, M., & Tajitsu, A. 2013, ApJ, 778, 146
  * Otsuka et al. (2010) Otsuka, M., Tajitsu, A., Hyung, S., & Izumiura, H. 2010, ApJ, 723, 658
  * Otsuka et al. (2003) Otsuka, M., Tamura, S., Yadoumaru, Y., & Tajitsu, A. 2003, PASP, 115, 67
  * Otsuka et al. (2017) Otsuka, M., Ueta, T., van Hoof, P. A. M., et al. 2017, ApJS, 231, 22
  * Owens (1967) Owens, J. C. 1967, Appl. Opt., 6, 51
  * Page et al. (2012) Page, M. J., Brindle, C., Talavera, A., et al. 2012, MNRAS, 426, 903
  * Parker et al. (2016) Parker, Q. A., Bojičić, I. S., & Frew, D. J. 2016, in Journal of Physics Conference Series, Vol. 728, Journal of Physics Conference Series, 032008
  * Peña et al. (1989) Peña, M., Ruíz, M. T., Maza, J., & González, L. E. 1989, Rev. Mexicana Astron. Astrofis., 17, 25
  * Peimbert (1978) Peimbert, M. 1978, in Planetary Nebulae, ed. Y. Terzian, Vol. 76, 215–224
  * Pena et al. (1990) Pena, M., Ruiz, M. T., Torres-Peimbert, S., & Maza, J. 1990, A&A, 237, 454
  * Pena et al. (1991) Pena, M., Torres-Peimbert, S., & Ruiz, M. T. 1991, PASP, 103, 865
  * Pereira & Miranda (2007) Pereira, C. B., & Miranda, L. F. 2007, A&A, 467, 1249
  * Poglitsch et al. (2010) Poglitsch, A., Waelkens, C., Geis, N., et al. 2010, A&A, 518, L2
  * Quireza et al. (2007) Quireza, C., Rocha-Pinto, H. J., & Maciel, W. J. 2007, A&A, 475, 217
  * Rauch (2003) Rauch, T. 2003, A&A, 403, 709
  * Rieke et al. (2004) Rieke, G. H., Young, E. T., Engelbracht, C. W., et al. 2004, ApJS, 154, 25
  * Rouleau & Martin (1991) Rouleau, F., & Martin, P. G. 1991, ApJ, 377, 526
  * Schlafly & Finkbeiner (2011) Schlafly, E. F., & Finkbeiner, D. P. 2011, ApJ, 737, 103
  * Schoenberner (1983) Schoenberner, D. 1983, ApJ, 272, 708
  * Stanghellini et al. (2008) Stanghellini, L., Shaw, R. A., & Villaver, E. 2008, ApJ, 689, 194
  * Stasińska et al. (2010) Stasińska, G., Morisset, C., Tovmassian, G., et al. 2010, A&A, 511, A44
  * Storey & Hummer (1995) Storey, P. J., & Hummer, D. G. 1995, MNRAS, 272, 41
  * Suda et al. (2011) Suda, T., Yamada, S., Katsuta, Y., et al. 2011, MNRAS, 412, 843
  * Tajitsu et al. (2023) Tajitsu, A., Otsuka, M., & Ueta, T. 2023, in prep.
  * Tody (1986) Tody, D. 1986, in Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series, Vol. 627, Instrumentation in astronomy VI, ed. D. L. Crawford, 733
  * Torres-Peimbert et al. (1990) Torres-Peimbert, S., Peimbert, M., & Pena, M. 1990, A&A, 233, 540
  * Ueta (2022) Ueta, T. 2022, Galaxies, 10, 30
  * Ueta & Otsuka (2021) Ueta, T., & Otsuka, M. 2021, PASP, 133, 093002
  * van Hoof et al. (2013) van Hoof, P. A. M., Van de Steene, G. C., Exter, K. M., et al. 2013, A&A, 560, A7
  * Vassiliadis & Wood (1993) Vassiliadis, E., & Wood, P. R. 1993, ApJ, 413, 641
  * Wesson (2016) Wesson, R. 2016, MNRAS, 456, 3774
  * Zhang (1995) Zhang, C. Y. 1995, ApJS, 98, 659
  * Zuo & Li (2014) Zuo, Z.-Y., & Li, X.-D. 2014, MNRAS, 442, 1980

##  Appendix A The method to calculate the nebula radiation

Figure A1:  The spectral energy distribution (SED) plots are generated from
the observed dereddened spectra and photometry, shown as gray lines and
circles respectively. The broken line represents the sum of the H i/He i/He ii
bound-free (B-F), free-free (F-F), and two-photon continuum calculated using
Nebcont in Starlink (Currie et al., 2014). We fitted the dust continuum in the
_gas-emission-free_ _Spitzer_ /IRS spectrum and _Herschel_ /PACS 70/100/160
photometry SED with a modified two-temperature (148 and 69.5 K) blackbody
function (dashed and dot lines, respectively). The total nebular continuum,
indicated by the black line, is the sum of the dust continuum fit and the
B-F/F-F/two-photon continuum.

We calculated the nebular radiation by following steps. First, we calculated
the total line flux from the observed atomic and molecular H2 lines. The
expected fluxes of H i and He ii lines appeared outside wavelengths in the
obtained spectra are calculated using H i and He ii recombination coefficients
in $T_{\rm e}$ = 12634 K and $n_{\rm e}$ = 420 cm-3 for H i and $T_{\rm e}$([O
iii]) and $n_{\rm e}$([S iii]) for He ii. The expected CEL fluxes are
calculated by solving atomic models of each ion in the obtained $T_{\rm e}$([O
iii], [N ii]) and $n_{\rm e}$([S ii], [S iii]). The following atomic lines are
considered: H i (Lyman $\alpha$; Balmer $n-2$, $n=3-40$, $n$: quantum number;
Paschen $n-3$, $n=4-40$; Bracket $n-4$, $n=5-40$; Pfund $n-5$, $n=6-40$), He
ii (Fowler $n-3$, $n=4-50$; Pickering $n-4$, $n=5-40$), CEL $[$C ii$]$ 157 ,
[N ii] 121/205 , [O i] 63/145 , [O iii] 51/88 , [O ii] 3726/29 Å, [Ne iii]
3868/3967 Å. Thus, we obtained the total atomic and H2 line flux of
$3.41\times 10^{-11}$ erg s-1 cm-2 by simply summing up the observed and
expected fluxes of the atomic and molecular H2 lines.

Second, we generated the integrated spectral energy distribution (SED) of the
H i/He i/He ii bound-free, free-free, and two-photon continuum using Nebcont
in Starlink (Currie et al., 2014) by adopting $I$(H$\beta$) = $4.54\times
10^{-13}$ erg s-1 cm-2, $T_{\rm e}$ = 12634 K, $n_{\rm e}$ = 420 cm-3, and
He+,2+ abundances (table 3).

Third, we fit a two temperature (149 and 69.5 K) modified blackbody function
to the dust continuum in the range from 14 to 160 based on the gas-emission-
free _Spitzer_ /IRS spectrum and the _Herschel_ /PACS photometry SED. Here, we
adopt $Q_{\rm ext}$ in the case of amorphous carbon grain (i.e., $Q_{\rm
ext}\propto{\lambda}^{-1}$). The expected flux contributions of the gas
emission to each 70/100/160 band were calculated by taking the PACS throughput
at each wavelength into account. The gas-emission-free PACS 70/100/160 band
flux densities are $9.65\times 10^{-14}$, $3.99\times 10^{-14}$, and
$7.44\times 10^{-15}$ erg s-1 cm-2 -1, respectively. Then, we added the
resultant fit to the nebular H i/He i/He ii continuum calculated in the second
step.

Thus, we obtained a synthesized nebular continuum SED (figure A1). For the
continuum flux in $5.2-36.0$ , we used the _Spitzer_ /IRS spectrum where the
gas emission was removed out but PAHs were not. Finally, we obtained the
nebular continuum-flux integrated between 0.1 and 600 to be $2.61\times
10^{-11}$ erg s-1 cm-2. The breakdown of the integrated nebular-line and
continuum fluxes is summarized in table A4.

##  Appendix B Supporting Tables

The follwings are the supproting tables.

Table A1: Broadband flux densities of H4-1.

$\lambda_{\rm c}$ |  $f$($\lambda$) | $F_{\lambda}$ | $I_{\lambda}$ | Sources | Ref.
---|---|---|---|---|---
() |  |  (erg s-1 cm-2 -1) |  (erg s-1 cm-2 -1) |  |
$0.1549$ | $~{}1.2383$ | $6.95\times 10^{-11}~{}\pm~{}4.04\times 10^{-13}$ | $2.02\times 10^{-10}~{}\pm~{}6.94\times 10^{-11}$ | GALEX | (1)
$0.2042$ | $~{}1.5521$ | $8.24\times 10^{-11}~{}\pm~{}1.09\times 10^{-12}$ | $2.79\times 10^{-10}~{}\pm~{}1.09\times 10^{-10}$ | XMM | (2)
$0.2305$ | $~{}1.4439$ | $5.35\times 10^{-11}~{}\pm~{}1.49\times 10^{-13}$ | $1.72\times 10^{-10}~{}\pm~{}6.44\times 10^{-11}$ | GALEX | (1)
$0.2934$ | $~{}0.5938$ | $2.28\times 10^{-11}~{}\pm~{}2.12\times 10^{-13}$ | $4.88\times 10^{-11}~{}\pm~{}1.19\times 10^{-11}$ | XMM | (2)
$0.3608$ | $~{}0.3455$ | $2.61\times 10^{-11}~{}\pm~{}2.21\times 10^{-14}$ | $4.96\times 10^{-11}~{}\pm~{}1.02\times 10^{-11}$ | SDSS | (3)
$0.4672$ | $~{}0.0538$ | $3.28\times 10^{-11}~{}\pm~{}1.47\times 10^{-14}$ | $5.42\times 10^{-11}~{}\pm~{}8.76\times 10^{-12}$ | SDSS | (3)
$0.6141$ | $-0.2409$ | $1.76\times 10^{-11}~{}\pm~{}7.84\times 10^{-15}$ | $2.53\times 10^{-11}~{}\pm~{}2.95\times 10^{-12}$ | SDSS | (3)
$0.7458$ | $-0.4168$ | $2.68\times 10^{-12}~{}\pm~{}2.16\times 10^{-15}$ | $3.54\times 10^{-12}~{}\pm~{}3.16\times 10^{-13}$ | SDSS | (3)
$0.8923$ | $-0.5828$ | $1.64\times 10^{-12}~{}\pm~{}3.72\times 10^{-15}$ | $2.00\times 10^{-12}~{}\pm~{}1.28\times 10^{-13}$ | SDSS | (3)
$1.0305$ | $-0.6694$ | $9.99\times 10^{-13}~{}\pm~{}1.63\times 10^{-15}$ | $1.17\times 10^{-12}~{}\pm~{}5.93\times 10^{-14}$ | UKIDSS | (4)
$1.2483$ | $-0.7572$ | $5.05\times 10^{-13}~{}\pm~{}1.10\times 10^{-15}$ | $5.67\times 10^{-13}~{}\pm~{}2.11\times 10^{-14}$ | UKIDSS | (4)
$1.6313$ | $-0.8422$ | $5.56\times 10^{-13}~{}\pm~{}2.50\times 10^{-15}$ | $6.00\times 10^{-13}~{}\pm~{}1.48\times 10^{-14}$ | UKIDSS | (4)
$2.2010$ | $-0.9026$ | $1.01\times 10^{-12}~{}\pm~{}2.93\times 10^{-15}$ | $1.06\times 10^{-12}~{}\pm~{}1.61\times 10^{-14}$ | UKIDSS | (4)
$3.6000$ | $-0.9750$ | $2.72\times 10^{-13}~{}\pm~{}2.14\times 10^{-14}$ | $2.75\times 10^{-13}~{}\pm~{}2.16\times 10^{-14}$ |  _Spitzer_ /IRAC | (5)
$4.5000$ | $-0.9830$ | $1.85\times 10^{-13}~{}\pm~{}2.57\times 10^{-14}$ | $1.86\times 10^{-13}~{}\pm~{}2.59\times 10^{-14}$ |  _Spitzer_ /IRAC | (5)
$5.8000$ | $-0.9880$ | $3.07\times 10^{-13}~{}\pm~{}6.99\times 10^{-14}$ | $3.09\times 10^{-13}~{}\pm~{}7.03\times 10^{-14}$ |  _Spitzer_ /IRAC | (5)
$8.0000$ | $-0.9840$ | $3.95\times 10^{-13}~{}\pm~{}3.43\times 10^{-14}$ | $3.98\times 10^{-13}~{}\pm~{}3.46\times 10^{-14}$ |  _Spitzer_ /IRAC | (5)
$11.5600$ | $-0.9720$ | $2.14\times 10^{-13}~{}\pm~{}1.68\times 10^{-15}$ | $2.16\times 10^{-13}~{}\pm~{}1.94\times 10^{-15}$ | WISE | (6)
$24.0000$ | $-0.9880$ | $2.68\times 10^{-13}~{}\pm~{}9.83\times 10^{-15}$ | $2.69\times 10^{-13}~{}\pm~{}9.90\times 10^{-15}$ |  _Spitzer_ /MIPS | (5)
$70.0000$ |  | $9.66\times 10^{-14}~{}\pm~{}1.94\times 10^{-15}$ | $9.66\times 10^{-14}~{}\pm~{}1.94\times 10^{-15}$ |  _Herschel_ /PACS | (7)
$100.0000$ |  | $4.07\times 10^{-14}~{}\pm~{}4.52\times 10^{-15}$ | $4.07\times 10^{-14}~{}\pm~{}4.52\times 10^{-15}$ |  _Herschel_ /PACS | (7)
$160.0000$ |  | $8.14\times 10^{-15}~{}\pm~{}1.83\times 10^{-15}$ | $8.14\times 10^{-15}~{}\pm~{}1.83\times 10^{-15}$ |  _Herschel_ /PACS | (7)

Note – The reddening correction for the flux densities at 70, 100, and 160 is
not performed because the dust extinction at these bands is negligible. The
gas-emission-free PACS 70/100/160 band flux densities are $9.65\times
10^{-14}$, $3.99\times 10^{-14}$, and $7.44\times 10^{-15}$ erg s-1 cm-2 -1,
respectively (see subsection 7.2.)
References – (1) Bianchi et al. (2017); (2) Page et al. (2012); (3) Alam et
al. (2015); (4) Lawrence et al. (2007); (5) Tajitsu et al. (2023); (6) Cutri
et al. (2013); (7) Herschel Point Source Catalogue Working Group et al.
(2020).

Table A2: The detected line-fluxes normalized to the line-flux of H$\beta$ whose intensity is 100.  $\lambda_{\rm lab}$ |  $f(\lambda)$ | Species |  $I$($\lambda$) |  $\delta$ $I$($\lambda$) | $\lambda_{\rm lab}$ |  $f(\lambda)$ | Species |  $I$($\lambda$) |  $\delta$ $I$($\lambda$)
---|---|---|---|---|---|---|---|---|---
1335.70 Å | $~{}1.584$ |  C ii | $107.376$ | $27.150$ | 7319/20 Å | $-0.398$ |  [O ii] | $4.119$ | $0.301$
1549.03 Å | $~{}1.238$ |  C iv | $243.453$ | $46.506$ | 7330/31 Å | $-0.400$ |  [O ii] | $3.514$ | $0.269$
1640.03 Å | $~{}1.177$ |  He ii | $68.000$ | $13.763$ | 8467.25 Å | $-0.541$ |  H i | $0.388$ | $0.101$
1664.03 Å | $~{}1.167$ |  O iii$]$ | $29.920$ | $8.179$ | 8502.48 Å | $-0.544$ |  H i | $0.561$ | $0.124$
1909.49 Å | $~{}1.257$ |  [C iii]+C iii$]$ | $1260.795$ | $243.188$ | 8545.38 Å | $-0.549$ |  H i | $0.625$ | $0.122$
2321/30 Å | $~{}1.382$ |  O iii$]$ | $1.817$ | $0.227$ | 8598.39 Å | $-0.554$ |  H i | $0.661$ | $0.181$
2322-28 Å | $~{}1.382$ |  [C ii] | $305.116$ | $66.494$ | 8665.02 Å | $-0.560$ |  H i | $0.714$ | $0.148$
4340.47 Å | $~{}0.157$ |  H i | $44.191$ | $1.545$ | 8728.90 Å | $-0.566$ |  C ii ? | $1.056$ | $0.123$
4363.21 Å | $~{}0.149$ |  [O iii] | $7.665$ | $0.959$ | 8750.47 Å | $-0.568$ |  H i | $1.182$ | $0.136$
4471.50 Å | $~{}0.115$ |  He i | $6.836$ | $0.760$ | 8862.78 Å | $-0.578$ |  H i | $1.519$ | $0.178$
4685.71 Å | $~{}0.050$ |  He ii | $10.068$ | $0.417$ | 9014.91 Å | $-0.590$ |  H i | $1.620$ | $0.179$
4712.42 Å | $~{}0.042$ |  He i+[Ar iv] | $0.770$ | $0.146$ | 9068.60 Å | $-0.594$ |  [S iii] | $0.856$ | $0.104$
4861.33 Å | $~{}0.000$ |  H i | $100.000$ | $1.121$ | 9229.01 Å | $-0.605$ |  H i | $2.666$ | $0.292$
4921.93 Å | $-0.016$ |  He i | $1.306$ | $0.243$ | 9537.45 Å | $-0.625$ |  [S iii] | $2.296$ | $0.474$
4958.91 Å | $-0.026$ |  [O iii] | $224.438$ | $3.212$ | 9545.97 Å | $-0.626$ |  H i | $3.844$ | $0.484$
5006.84 Å | $-0.038$ |  [O iii] | $670.451$ | $9.689$ | 9824.13 Å | $-0.643$ |  [C i] | $4.142$ | $0.710$
5198/200 Å | $-0.082$ |  [N i] | $1.587$ | $0.127$ | 9850.26 Å | $-0.644$ |  [C i] | $14.058$ | $1.571$
5411.52 Å | $-0.126$ |  He ii | $0.902$ | $0.116$ | 10049.37 Å | $-0.656$ |  H i | $4.361$ | $0.724$
5754.60 Å | $-0.185$ |  [N ii] | $2.406$ | $0.156$ |  5.51  | $-0.987$ |  H2 | $5.976$ | $0.917$
5875.66 Å | $-0.203$ |  He i | $12.539$ | $0.488$ |  6.11  | $-0.989$ |  H2 | $2.324$ | $0.359$
6300.34 Å | $-0.263$ |  [O i]+[S iii] | $9.208$ | $0.408$ |  6.91  | $-0.990$ |  H2 | $10.573$ | $1.621$
6363.78 Å | $-0.271$ |  [O i] | $3.110$ | $0.153$ |  9.67  | $-0.949$ |  H2 | $8.904$ | $1.312$
6548.10 Å | $-0.296$ |  [N ii] | $33.293$ | $3.103$ |  10.52  | $-0.960$ |  [S iv] | $1.010$ | $0.170$
6562.77 Å | $-0.298$ |  H i | $277.865$ | $13.369$ |  12.28  | $-0.980$ |  H2 | $2.254$ | $0.345$
6583.50 Å | $-0.300$ |  [N ii] | $105.668$ | $5.310$ |  12.40  | $-0.980$ |  H i | $0.871$ | $0.153$
6678.16 Å | $-0.313$ |  He i | $4.274$ | $0.226$ |  12.81  | $-0.983$ |  [Ne ii] | $0.325$ | $0.062$
6716.44 Å | $-0.318$ |  [S ii] | $0.957$ | $0.070$ |  15.56  | $-0.985$ |  [Ne iii] | $2.971$ | $0.456$
6730.82 Å | $-0.320$ |  [S ii] | $0.929$ | $0.072$ |  17.03  | $-0.982$ |  H2 | $2.070$ | $0.317$
7065.25 Å | $-0.364$ |  He i | $4.873$ | $0.282$ |  18.71  | $-0.981$ |  [S iii] | $0.764$ | $0.124$
7135.80 Å | $-0.374$ |  [Ar iii] | $0.487$ | $0.056$ |  25.89  | $-0.989$ |  [O iv] | $25.845$ | $3.930$
7236.42 Å | $-0.387$ |  C ii+[Ar iv] | $0.521$ | $0.108$ |  28.22  | $-0.990$ |  H2 | $2.043$ | $0.315$
7281.35 Å | $-0.393$ |  He i | $0.769$ | $0.089$ |  |  |  |  |
Table A3: Comparison between the observed and the Cloudy model predicted
values.

Species |  $\lambda_{\rm lab}$ (Å) |  $I$(Obs) |  $I$(1.25 M⊙) |  $I$(1.50 M⊙) | Species |  $\lambda_{\rm lab}$ (Å) |  $I$(Obs) |  $I$(1.25 M⊙) |  $I$(1.50 M⊙)
---|---|---|---|---|---|---|---|---|---
C iv | $1549.03$ | $242.730$ | $320.547$ | $320.043$ |  [N ii] | $6583.50$ | $110.303$ | $107.786$ | $106.932$
He ii | $1640.03$ | $73.154$ | $87.047$ | $88.555$ |  He i | $6678.16$ | $4.271$ | $3.619$ | $3.610$
$[$C iii$]$+C iii$]$ | $1909.00$ | $1254.077$ | $1250.051$ | $1228.957$ |  [S ii] | $6716.44$ | $0.979$ | $0.971$ | $1.024$
[O iii] | $2320.95$ | $2.098$ | $2.493$ | $2.463$ |  [S ii] | $6730.82$ | $0.948$ | $1.147$ | $1.214$
$[$C ii$]$ | $2326.00$ | $293.310$ | $273.891$ | $277.782$ |  [Ar iii] | $7135.80$ | $0.461$ | $0.428$ | $0.447$
He ii | $4338.55$ | $0.245$ | $0.310$ | $0.316$ |  He i | $7281.35$ | $0.729$ | $0.909$ | $0.909$
H i | $4340.47$ | $47.394$ | $47.224$ | $47.218$ |  [O ii] | $7320.00$ | $4.102$ | $8.109$ | $8.417$
[O iii] | $4363.21$ | $8.851$ | $9.932$ | $9.811$ |  [O ii] | $7330.00$ | $3.458$ | $6.490$ | $6.736$
He ii | $4685.71$ | $10.419$ | $10.954$ | $11.143$ |  H i | $8467.25$ | $0.393$ | $0.443$ | $0.442$
He ii | $4859.18$ | $0.505$ | $0.618$ | $0.628$ |  H i | $8502.48$ | $0.562$ | $0.521$ | $0.520$
H i | $4861.33$ | $100.000$ | $100.000$ | $100.000$ |  H i | $8545.38$ | $0.664$ | $0.621$ | $0.620$
He i | $4921.93$ | $1.457$ | $1.322$ | $1.318$ |  H i | $8598.39$ | $0.634$ | $0.751$ | $0.750$
[O iii] | $4958.91$ | $217.864$ | $211.798$ | $208.263$ |  H i | $8665.47$ | $0.712$ | $0.923$ | $0.922$
[O iii] | $5006.84$ | $651.687$ | $631.914$ | $621.370$ |  H i | $8750.43$ | $1.233$ | $1.154$ | $1.152$
He ii | $5411.52$ | $0.804$ | $0.939$ | $0.955$ |  H i | $8862.78$ | $1.396$ | $1.472$ | $1.470$
[N ii] | $5754.60$ | $2.617$ | $2.644$ | $2.622$ |  H i | $9014.91$ | $1.578$ | $1.694$ | $1.692$
He i | $5875.66$ | $14.077$ | $13.735$ | $13.716$ |  [S iii] | $9068.60$ | $0.879$ | $0.774$ | $0.794$
[O i] | $6363.78$ | $3.577$ | $5.878$ | $6.268$ |  H i | $9229.01$ | $2.568$ | $2.332$ | $2.328$
[N ii] | $6548.10$ | $35.768$ | $36.565$ | $36.275$ |  [S iii] | $9537.45$ | $2.229$ | $1.911$ | $1.961$
He ii | $6559.91$ | $1.382$ | $1.529$ | $1.555$ |  $[$C i$]$ | $9824.13$ | $3.548$ | $2.653$ | $2.799$
H i | $6562.77$ | $278.955$ | $286.983$ | $287.272$ |  $[$C i$]$ | $9850.26$ | $13.766$ | $7.926$ | $8.363$
Species |  $\lambda_{\rm lab}$ () |  $I$(Obs) |  $I$(1.25 M⊙) |  $I$(1.50 M⊙) | Species |  $\lambda_{\rm lab}$ () |  $I$(Obs) |  $I$(1.25 M⊙) |  $I$(1.50 M⊙)
H2 | $9.67$ | $9.722$ | $9.521$ | $9.573$ |  [Ne iii] | $15.56$ | $3.272$ | $2.402$ | $2.225$
[S iv] | $10.52$ | $1.105$ | $0.921$ | $0.908$ |  H2 | $17.03$ | $2.279$ | $2.111$ | $2.121$
H2 | $12.28$ | $2.479$ | $2.543$ | $2.573$ |  [S iii] | $18.71$ | $0.933$ | $0.936$ | $0.962$
[Ne ii] | $12.81$ | $0.358$ | $0.466$ | $0.499$ |  [O iv] | $25.89$ | $28.494$ | $32.949$ | $32.718$
$F_{\nu}$ |  $\lambda_{\rm lab}$ () |  $F_{\nu}$(Obs) |  $F_{\nu}$(1.25 M⊙) |  $F_{\nu}$(1.50 M⊙) | $F_{\nu}$ |  $\lambda_{\rm lab}$ () |  $F_{\nu}$(Obs) |  $F_{\nu}$(1.25 M⊙) |  $F_{\nu}$(1.50 M⊙)
IRS01 | $10.00$ | $0.003$ | $0.003$ | $0.003$ | IRS13 | $25.00$ | $0.063$ | $0.059$ | $0.058$
IRS02 | $14.00$ | $0.012$ | $0.009$ | $0.010$ | IRS14 | $26.00$ | $0.068$ | $0.065$ | $0.063$
IRS03 | $15.00$ | $0.015$ | $0.011$ | $0.011$ | IRS15 | $27.00$ | $0.074$ | $0.070$ | $0.067$
IRS04 | $16.00$ | $0.019$ | $0.015$ | $0.015$ | IRS16 | $28.00$ | $0.078$ | $0.075$ | $0.072$
IRS05 | $17.00$ | $0.023$ | $0.021$ | $0.021$ | IRS17 | $29.00$ | $0.083$ | $0.080$ | $0.077$
IRS06 | $18.00$ | $0.027$ | $0.023$ | $0.023$ | IRS18 | $30.00$ | $0.087$ | $0.085$ | $0.082$
IRS07 | $19.00$ | $0.032$ | $0.027$ | $0.027$ | IRS19 | $31.00$ | $0.091$ | $0.090$ | $0.086$
IRS08 | $20.00$ | $0.037$ | $0.031$ | $0.031$ | IRS20 | $32.00$ | $0.095$ | $0.095$ | $0.091$
IRS09 | $21.00$ | $0.042$ | $0.036$ | $0.036$ | IRS21 | $33.00$ | $0.097$ | $0.099$ | $0.095$
IRS10 | $22.00$ | $0.047$ | $0.042$ | $0.041$ | PAC01 | $70.00$ | $0.158$ | $0.153$ | $0.149$
IRS11 | $23.00$ | $0.053$ | $0.048$ | $0.047$ | PAC02 | $100.00$ | $0.133$ | $0.133$ | $0.134$
IRS12 | $24.00$ | $0.058$ | $0.054$ | $0.052$ | PAC03 | $160.00$ | $0.066$ | $0.072$ | $0.074$
Band |  $\lambda_{\rm lab}$ () |  $I$(Obs) |  $I$(1.25 M⊙) |  $I$(1.50 M⊙) | Band |  $\lambda_{\rm lab}$ () |  $I$(Obs) |  $I$(1.25 M⊙) |  $I$(1.50 M⊙)
PAH1 | $8.85$ | $52.737$ | $52.481$ | $50.710$ | DUST1 | $22.50$ | $323.701$ | $292.249$ | $286.890$
PAH2 | $11.30$ | $42.644$ | $41.541$ | $42.593$ | DUST2 | $28.50$ | $196.829$ | $190.268$ | $183.637$
PAH3 | $12.30$ | $48.847$ | $51.997$ | $53.343$ |  |  |  |  |

Note – The line and continuum band flux $I$ is normalized to the H$\beta$ flux
$I$(H$\beta$) = 100. The unit of the dust continuum flux density $F_{\nu}$ is
Jy. PAH1, PAH2, and PAH3 are the PAH band fluxes integrated in $8.30-9.40$ ,
$11.00-11.60$ , and $11.80-12.80$ , respectively. DUST1 and DUST2 are the dust
continuum fluxes integrated in $22.00-25.00$ and $27.00-30.00$ , respectively.

Table A4: The breakdown of the integrated nebular-line and continuum fluxes.  Line band () | Flux (erg s-1 cm-2)
---|---
$0.1~{}-~{}205$ | $3.41\times 10^{-12}$
Continuum band () | Flux (erg s-1 cm-2)
$0.1~{}-~{}5.2$ | $7.92\times 10^{-12}$
$5.2~{}-~{}36.0$ | $8.22\times 10^{-12}$
$36.0~{}-~{}600$ | $9.97\times 10^{-12}$
Total | $6.02\times 10^{-11}$