# Modeling Ionized Gas in the Small Magellanic Cloud: The Wolf-Rayet Nebula
N76

Elizabeth Tarantino Space Telescope Science Institute, 3700 San Martin Drive,
Baltimore, MD 21218, USA  Department of Astronomy, University of Maryland,
College Park, MD 20742, USA  Alberto D. Bolatto Department of Astronomy,
University of Maryland, College Park, MD 20742, USA  Rémy Indebetouw
Department of Astronomy, University of Virginia, Charlottesville, VA 22904,
USA  National Radio Astronomy Observatory, 520 Edgemont Road, Charlottesville,
VA 22903, USA  Mónica Rubio Departamento de Astronomía, Universidad de Chile,
Casilla 36-D, Santiago, Chile  Karin M. Sandstrom Center for Astrophysics &
Space Sciences, Department of Physics, University of California, San Diego,
9500 Gilman Drive, San Diego, CA 92093, USA  J.-D T. Smith Dept. of Physics &
Astronomy, University of Toledo, Toledo, OH 43606, USA  Daniel Stapleton Space
Telescope Science Institute, 3700 San Martin Drive, Baltimore, MD 21218, USA
Department of Astronomy, University of Maryland, College Park, MD 20742, USA
Mark Wolfire Department of Astronomy, University of Maryland, College Park, MD
20742, USA Elizabeth Tarantino<EMAIL_ADDRESS>
(Accepted April 10, 2024)

###### Abstract

We present Cloudy modeling of infrared emission lines in the Wolf-Rayet (WR)
nebula N76 caused by one of the most luminous and hottest WR stars in the low
metallicity Small Magellanic Cloud. We use spatially resolved mid-infrared
Spitzer/IRS and far-infrared Herschel/PACS spectroscopy to establish the
physical conditions of the ionized gas. The spatially resolved distribution of
the emission allows us to constrain properties much more accurately than using
spatially integrated quantities. We construct models with a range of constant
hydrogen densities between $\rm n_{H}=4-10\ \mbox{$\mathrm{cm^{-3}}$}$ and a
stellar wind-blown cavity of 10 pc which reproduces the intensity and shape of
most ionized gas emission lines, including the high ionization lines [O IV]
and [Ne V], as well as [S III], [S IV], [O III], and [Ne III]. Our models
suggest that the majority of [Si II] emission (91%) is produced at the edge of
the HII region around the transition between ionized and atomic gas while very
little of the [C II] emission ($<5\%$) is associated with the ionized gas. The
physical conditions of N76 are characterized by a hot HII region with a
maximum electron temperature of $T_{e}\sim 24,000\rm\,K$, electron densities
that range from $n_{e}\sim 4\ \rm to\ 12\rm\,cm^{-3}$, and high ionization
parameters of $\log(U)\sim-1.15\ \rm to\ -1.77$. By analyzing a low metallicty
WR nebula with a single ionization source, this work gives valuable insights
on the impact WR stars have to the galaxy-integrated ionized gas properties in
nearby dwarf galaxies.

Interstellar medium (847), H II regions (694), Small Magellanic Cloud (1468),
Wolf-Rayet stars (1806), Photoionization (2060)

††journal: ApJ††software: astropy (The Astropy Collaboration 2013, 2018,
2022), CUBISM (v1.8; Smith et al. 2007a), TinyTim (Krist et al. 2011), Cloudy
(Ferland et al. 2013), PAHFIT (Smith et al. 2007b)

##  1 Introduction

Massive stars inject energy into the surrounding gas, heating and influencing
the chemistry of their nearby interstellar medium (ISM). The effects massive
stars have on the ISM depend on a variety of factors, including the
metallicity of the surrounding material. Lower metallicity stars have harder
radiation fields and a decreased mass-loss rate (e.g., Hurley et al., 2000;
Vink et al., 2001). At the same time, the low metallicity ISM contains less
dust (e.g., Galliano et al., 2021), which consequently leads to a more porous
structure with a lower HII region effective covering fraction allowing
radiation fields to penetrate deeper into the ISM (e.g., Poglitsch et al.,
1995; Cormier et al., 2015, 2019). The chemistry and physical processes
governing the interaction between massive stars is thus highly affected by the
metallicity of the medium.

Investigation of ionized gas tracers are used to determine the properties of
the ionized gas across galaxies. In particular, the infrared (IR) lines such
as [Ne II], [S III], [O III], [Ne III], and [S IV], accessed by the Spitzer,
Herschel, and James Webb space telescopes trace a variety of conditions in the
ISM, including the ionizing radiation field strength and hardness, the
ionization parameter, the density and temperature of the ionized gas, and the
heating mechanism of the gas (e.g., Baldwin et al., 1981; Kaufman et al.,
2006; Cormier et al., 2015, 2019; Polles et al., 2019). The IR fine-structure
lines are less affected by extinction and dust attenuation than optical lines
and are therefore ideal for studying star-forming regions.

Previous studies of mostly spatially unresolved infrared spectroscopy in low
metallicity star-forming galaxies find harder radiation fields, extended,
bright [O III] emission, and an overall more porous structure (Hunt et al.,
2010; Cormier et al., 2015, 2019). The majority of spectra analyzed in these
studies, however, are unresolved and contain emission from the different
phases of the ISM and a variety of HII regions that vary in properties
(density, ionization parameter, stellar population, etc.) blended into one
beam or observation. Spatially resolved studies that focus on individual HII
regions, on the other hand, bring additional information and can be
considerably more powerful at constraining physical conditions. Of course,
they can only be carried out in nearby sources. The Small Magellanic Cloud
(SMC), at one-fifth solar metallicity and only 63 kpc (Dufour, 1984; Russell &
Dopita, 1992) away provides the ideal laboratory for studying the low
metallicity ionized gas in spatially-resolved detail.

This paper focuses on studying the ionized gas emission around the class of
the hottest and most luminous stars in the local universe, a Wolf-Rayet (WR)
star. After an O-type star loses its hydrogen-rich envelope (either from
binary stripping and accretion, or via stellar winds removing the outer layers
of the star), the inner, hot core is exposed and forms a WR star (Crowther &
Hadfield, 2006). WR stars are characterized by effective temperatures of
$T_{*}\sim$ 50 - 110 kK, luminosities of $\log L\sim 5-6.2\mbox{\,$\rm
L_{\odot}$}$, have strong stellar winds with velocities that range from $\rm
v_{\infty}\sim 1500-5000$ km s-1, and mass outflow rates of $\rm\dot{M}\sim
10^{-4}-10^{-5}\mbox{\,M${}_{\odot}$}\ yr^{-1}$ (Crowther & Hadfield, 2006).
These stellar winds compress the surrounding material, forming a stellar wind
blown bubble (Weaver et al., 1977). Simultaneously, the bright ultraviolet
(UV) flux from the WR star creates a highly ionized HII region, forming a
structure often called a WR nebula, a subclass of HII regions (e.g., Chu,
1981). Since WR stars inject massive amounts of energy into the surrounding
medium, they can have a profound influence on shaping the nearby ISM and
energy exchange in stellar clusters (e.g., Sokal et al., 2016).

This paper studies the WR nebula N76 that harbors one of the hottest, most
luminous WR stars in the SMC (Shenar et al., 2016). First classified by Henize
(1956) as an HII region, Garnett et al. (1991) identified broad He II emission
in N76, indicative of a WR nebula. N76 is powered by a well-studied binary
system consisting of a WN4 and a O6 I(f) star (Niemela et al., 2002; Shenar et
al., 2016). The goal of this work is to model the physical conditions of the
ionized gas in N76 using mid and far infrared spectroscopy and the Cloudy
photoionization code (Ferland et al., 2017) to provide a detailed, resolved
view of the impact the most energetic stars have on the metal-poor ISM. Many
low metallicity dwarf galaxies that have been studied in the past (Wu et al.,
2006; Hunt et al., 2010; Cormier et al., 2015, 2019) contain a large number of
WR stars and are undergoing a burst of active star formation, but these
galaxies are too far away to examine the specific impact of WR stars in
detail. Therefore N76 provides the ideal case study to isolate the impact WR
stars have on the low metallicity ISM due to its proximity and simple, well
characterized single ionization source.

Observations of Milky Way WR nebulae show that they vary in morphology, but
often appear as thin, bubble structures, or disrupted shells (Chu, 1981; Chu
et al., 1983; Toalá et al., 2015). In particular, the well-studied WR nebula
NGC 6888 shows a double shell model with a denser inner shell ($\rm n_{e}\sim
400\,\mbox{$\mathrm{cm^{-3}}$}$) and a thinner outer shell ($\rm n_{e}\sim
180\,\mbox{$\mathrm{cm^{-3}}$}$) (Fernández-Martín et al., 2012; Rubio et al.,
2020). Due to the rarity of Wolf-Rayet nebulae (e.g.,, Conti & Massey, 1989;
Hainich et al., 2014; Neugent et al., 2018), there have been few observations
investigating the properties of WR nebulae at low metallicities. It is
possible that the lower dust abundance and harder radiation fields at low
metallicities will allow high energy photons to penetrate deeper into the low
metallicity ISM, forming larger, lower density nebulae with high ionizing
photon escape fractions.

Further, understanding how WR stars directly affect their surrounding
environment and host galaxies is essential to deciphering the nature of
extreme emission line galaxies (EELGs). These are local analogues to galaxies
during the Epoch of Reionization, characterized by bright emission lines from
high ionization species (such as [O III], C III], C IV, He II, and [Ne V]),
low metallicities, and highly ionized gas (e.g., Atek et al., 2011; Maseda et
al., 2013; Rigby et al., 2015; Senchyna et al., 2017; Berg et al., 2019;
Olivier et al., 2022). The ionization source for the extreme emission found in
these galaxies is unclear, but ultraluminous X-ray sources and radiation from
WR stars are possible explanations. The observations of N76 presented in this
work will directly test whether low metallicity WR nebula can produce the high
ionization species ([O IV] and [Ne V]) associated with EELGs.

Below we outline the sections of this paper. In Section 2, we describe the
Spitzer and Herschel observations and the methods of producing line-integrated
spatially resolved emission line maps. Section 3 describes the Cloudy
photoionization model input conditions and N76 geometry. We then compare the
photoionization models to the spatially resolved IR data in Section 4. We
discuss the implications of these results, putting N76 in context with other
low metallicity galaxies and WR nebulae in Section 5. Lastly, our conclusions
are described in Section 6.

##  2 Observations

###  2.1 Spitzer IRS observations and spectral line images

The Spitzer data come from the InfRared Spectrograph (IRS) as part of the
Spitzer Spectroscopic Survey of the Small Magellanic Cloud (S4MC). The full
observations are described by Sandstrom et al. (2012), but we summarize them
briefly here. We downloaded the raw Basic Calibrated Data (BCD) files from the
Spitzer Heritage Archive from project GO 30491 (PI: A. Bolatto) with pipeline
version S18.18. We use the long-low (LL) and short-low (SL) modules of the IRS
to cover the 5.2 µm - 38.4 µm range with a resolving power ranging from $R\sim
60-120$ (coverage of the modules is shown in Figure 3 of Sandstrom et al.,
2012). The spatially resolved maps are produced by stepping the IRS slit
perpendicular and parallel to the source in steps of half a slit width for the
LL module (5″.08) and the full slit width for the SL module (3″.7). The LL
maps for N76 had $75\times 6$ steps ($376\arcsec\times$ 395″) with a 14 s
integration time per position. The SL maps are created similarly, except the
steps are using the full slit width in order to increase coverage and have
$120\times 5$ steps ($220\arcsec\times$ 208″).

Each set of LL and SL observations has an associated background observation
taken at R.A. $\rm 1^{h}9^{m}40^{s}$ Dec. $-73^{\circ}31\arcmin 30\arcsec$
(J2000), a region chosen from the Multiband Imaging Photometer for Spitzer
(MIPS) and Infrared Array Camera (IRAC) observations of the SMC to have
negligible emission from the SMC itself (Bolatto et al., 2007; Leroy et al.,
2007; Sandstrom et al., 2010). We subtract the background from the BCD files,
which removes emission from the zodiacal light of the Milky Way cirrus, and
reduces the number of “hot” pixels that contaminate the IRS data.

####  2.1.1 Constructing the data cubes

The BCD files are assembled into cubes using the software CUBISM (Smith et
al., 2007a, we use version 1.8). CUBISM processes the two-dimensional slit
spectral images created by scanning the IRS slit across the source into three-
dimensional spectral data cubes through polygon-clipping-based re-projection.
This algorithm is flux-conserving and is specially adapted for the IRS on
Spitzer. We apply a slit-loss correction for extended sources which assumes
uniform extended emission inside and outside the slits and corrects the native
tuning of the IRS pipeline that assumes point sources. The IRS is susceptible
to hot, “rogue”, or bad pixels that appear on the cube images as repeating
stripes, resulting from the bad pixel scanning across the image through the
grid of positions. These pixels are often caused by interactions between the
IRS CCD and solar wind particles. The positions and intensity of the rogue bad
pixels can vary on scales of hours to days, making them unpredictable. Some
bad pixels can be flagged automatically through CUBISM, but many must be
flagged by hand. We perform extensive flagging of the bad pixels in each data
cube ($\sim~{}7\%$ of data), with particular focus on the emission from the
spectral lines we plan to model, in order to produce the cleanest spectral
line maps.

Figure 1: Images of infrared line emission in N76 where the color scale
corresponds to the line brightness in $\rm 10^{-8}\ W\ m^{-2}\ sr^{-1}$. The
location of the WR-O star binary AB7, the ionization source of N76, is
illustrated by a gray star in the center of each image. The supernova remnant
(SNR) E0102-72 is located NE, on the upper left of N76, and is seen in some of
the emission lines ([O IV], [Ne V], [Ne III], [Ne II]). The full width at half
maximum (FWHM) of the point spread function (PSF) is shown in the upper right
corner of each panel, and the ionization potential of the ion is displayed in
the lower left corner together with the wavelength of each transition. These
images show how the structure of N76 changes depending on the ion.

####  2.1.2 Fitting spectral line maps

We use the program PAHFIT (Smith et al., 2007b) to create line maps for each
emission line in the IRS bandpass. PAHFIT uses a physically motivated model to
simultaneously fit multiple components in an IRS spectrum. Its model includes
dust continuum in fixed equilibrium temperature bins, starlight, bright
emission lines, individual and blended PAH features, and extinction from
silicate grains. In our fits we do not include extinction because the SMC has
very low levels of mid-IR dust absorption, due to its relatively low gas-to-
dust ratio (Lee et al., 2009; Roman-Duval et al., 2014). The mapping strategy
produces maps for both the SL and LL modules, but the placement of the
apertures corresponding to each spectral order creates an offset of the
mapping area between the different modules and orders, resulting in different
coverage for SL1, SL2, LL1, and LL2 (see Fig. 3 in Sandstrom et al., 2012, for
an illustration of the coverage on N 76). To maximize the area of the line
maps, we fit each IRS module and spectral order separately with PAHFIT.

In order to create the emission line map, we first extract the spectrum for
each pixel in the data cube. The native pixel size for the LL cubes is 5.07″,
which is smaller than the PSF at the lowest wavelength line ([S III] 18 µm ,
PSF FWHM of 4.15″) and highest wavelength ([Si II], PSF FWHM of 7.61″). The SL
native pixel scale is 1.85″, which samples the PSF of the lowest wavelength
line ([S IV], PSF FWHM of 2.40″) and the highest wavelength line ([Ne II], PSF
FWHM of 2.83″). We pass the spectrum into PAHFIT with a custom emission line
fitting list to ensure the fainter lines not included in the default fitting
list are incorporated. Then we extract the result of the line fit and
construct the resulting line map using the fits for the individual pixels. The
line maps of all emission lines (excluding PAHs and quadrupole H2 transitions)
found in N76 are presented in Figure 1. The integrated spectrum of N76 is
presented in Figure 2.

In addition to the N76 nebula, the IRS maps also include emission from a well-
studied nearby core-collapse supernova remnant SNR E0102-72 (Hayashi et al.,
1994; Blair et al., 2000; Chevalier, 2005). The focus of this work is on N76,
so we refer the reader to Sandstrom et al. (2009) which investigates the mid-
infrared properties of SNR E0102-72 in detail with the Spitzer IRS dataset
also used in this work. The SNR is masked out for all analysis in this paper.

###  2.2 Herschel PACS observations and spectral line images

In addition to the Spitzer IRS data, we also report far-infrared observations
from the Herschel Photodetector Array Camera and Spectrometer (PACS) project
1431 (PI: R. Indebetouw). We downloaded the data from the ESA Herschel Science
Archive and use the level 2 pipeline reduced data, which contains spectral
data cubes of the [C II] 158 µm, [O I] 63 µm, and [O III] 88 µm lines. To
produce the spectral line maps, we fit a polynomial of order one to the
baseline, remove the continuum, and integrate over the spectral line. This
process was repeated pixel-by-pixel over the full datacube to produce a line
intensity image for [C II], [O I], and [O III]. Lastly, we convert to the same
surface brightness units $\rm 10^{-8}\ W\ m^{-2}\ sr^{-1}$ used with the IRS.
We compared this reduction to an independent reduction using the procedure in
Cormier et al. (2015) and found excellent agreement within the uncertainties
(priv comm, D. Cormier).

###  2.3 Estimation of uncertainty in emission line surface brightness

For the Spitzer IRS data, CUBISM provides a uncertainty cube that corresponds
to the statistical uncertainty of the IRS measurement. These uncertainties are
then propagated into the PAHFIT code, where they are used for the chi-squared
minimization fitting. The resulting uncertainty from the PAHFIT output is
propagated into our emission line images. These estimates of the
uncertainties, however, are lower limits since they only include noise from
the detector readout and do not comprise any systematic errors such as
uncertainties in the absolute calibration or the effective PSF due to
perpendicular slit offsets in the map.

For the Herschel PACS data, we compute the RMS noise outside of the signal
integration per each pixel as a 1$\sigma$ statistical uncertainty for each
channel. We then propagate the channel RMS to calculate an uncertainty for the
integrated signal. This process is repeated per pixel to create a map of the
statistical uncertainties for [C II], [O I], and [O III].

In addition to the statistical uncertainties, we also calculate a “background”
level for each emission line. We concentrate on modeling the emission from the
HII region created by the Wolf-Rayet/O star binary ionization source AB7.
There are other ionization sources across the SMC that can produce diffuse
ionized gas throughout the galaxy (e.g., Haffner et al., 2009). In order to
quantify and remove this extraneous emission, we identify a blank part of the
map that is not contaminated by emission from the HII region (typically
$\sim$80 pc away from AB7). We calculate the median surface brightness in this
patch of sky and consider it the “background” for the given emission line. The
background is very low, typically $\lesssim 5\%$ of the emission line
intensity. The background levels are generated for the IRS maps only because
the PACS field of view only covers the N76 nebula and there are no regions in
the PACS maps that would not be contaminated by N76.

Figure 2: Left: Integrated spectrum of the N76 region showing the ionized gas
lines in N76 we model with Cloudy. Line labels are in blue and the orders of
the IRS instrument are presented through the background color. Right: spectrum
at a single line of sight towards the location of the WR binary AB7 that
highlights the [Ne V] and possible [S I] features.

###  2.4 PSF matching

The spatial resolution of the Spitzer IRS changes as a function of wavelength
as the telescope is diffraction limited, with a point spread function (PSF)
full width at half maximum (FWHM) ranging from 7.6″at 34.8 µm for [Si II] to
2.4″at 10.5 µm for [S IV]. In order to properly model the spectra, we need to
convolve the spectral line images to a common PSF (e.g., Smith et al., 2009;
Sandstrom et al., 2009). We follow the procedure described by Aniano et al.
(2011) to create custom convolution kernels that will transform an image with
a narrower PSF to a broader PSF. We refer the reader to Appendix A for a
detailed description of the PSF matching.

###  2.5 Radial profiles

The spatial information given by the resolved infrared emission lines helps
differentiate between photoionization models. Cloudy simulates HII regions
with a 1-dimensional (1D) radiative transfer code, but can accommodate a
variety of HII region geometries with the assumptions in a given model. Due to
the spherical nature of N76 (see Figure 1), we use the spherical geometry in
Cloudy to best match the observations. In order to properly compare with the
models, we average the spectral line images to create surface brightness
radial profiles. N76 is an ideal region for this method because it is
approximately spherically symmetric and mostly isolated from other bright
regions (see Figure 1). We use the location of AB7 as the center of N76 and
compute an azimuthal average with a radial bin size corresponding to one pixel
or $5\arcsec$. The radial profiles show how the line emission varies with
distance away from the ionization source, thereby preserving some of the
spatial information from the spectral line images but making it simple to
compare to properly integrated 1D Cloudy models. The radial profiles measured
for each line are presented in Figure 3, showing that higher ionization lines
generally peak closer to AB7. We also report the standard deviation in a given
azimuthal bin as the gray uncertainty bars in the profile.

###  2.6 Integrated Intensities

In addition to examining the resolved surface brightness of the infrared
emission lines through radial profiles ($\S$2.5), we also record the total
intensities by integrating across the nebula for the emission lines from
Spitzer IRS. We aim to compare the observed total intensity of N76 for each
line to the predictions of the total intensity from Cloudy as a secondary
check of the models. The lines from Herschel PACS ([C II], [O I], and [O
III]), however, only cover a strip of the N76 nebula and not correcting for
this missing area and flux leads to under-predictions of the total intensity
of N76. We identify the missing area for each PACS line by using the Spitzer
IRS emission line that has the most similar radial profile. For [C II] and [O
I] we use [Si II], and for [O III] we use [Ne III]. To correct for the missing
flux, we multiply the flux contained in the PACS footprint by the ratio of
total area of N76 for that emission line to the area of the PACS footprint. We
note that this modified version of the PACS maps are only used for the
integrated total intensities reported in Table 1 and are not used for the
radial profiles. The uncertainties in the integrated intensities are derived
from the statistical uncertainties propagated through the total intensity
calculation (see $\S$2.3). The integrated intensities for all the emission
lines are reported in Table 1.

###  2.7 [O IV] and [Fe II] Line Blend

The spectral resolution of R$\sim$70 in the Spitzer IRS LL module causes the
[O IV] 25.89 µm and [FeII] 25.99 µm lines to be blended. For this work, we
assume that the flux of the line observed at 25.9 µm is dominated by [O IV]
emission. The Cloudy models are tuned to the parameters of the AB7 system
determined from modeling the optical and UV spectra in Shenar et al. (2016)
(see §3) and find that [O IV] 25.89 µm emission is two orders of magnitude
brighter than [Fe II] 25.99µm line, strongly indicating that [Fe II] makes a
negligible contribution to the 25.9 µm. The Cloudy models are inconsistent
with a large contribution of [Fe II] to the 25.9 µm blend due to the nature of
hot WR star as the ionization source. Further, the resolved images of N76 show
that the emission at 25.9 µm is centrally concentrated and associated with the
high ionization zone of the HII region, making it much more likely to be a
line with an ionization potential of 55 eV ([O IV]) rather than 7.9 eV ([Fe
II]). There are also no detections of the iron lines at [Fe II] 17.94 µm or
[Fe III] 22.93 µm in N76. Baldovin-Saavedra et al. (2011) use the high
resolution Spitzer module (LH) which can resolve the 25.9 µm [O IV] and [Fe
II] blend. In their observations of the Taurus Molecular Cloud, they find that
the [Fe II] 25.99 µm line is a factor of three brighter than [Fe II] 17.94 µm,
consistent with theoretical predictions in Hollenbach & McKee (1989). Using
this relationship, we calculate an upper limit of the [Fe II] 25.99 µm in N76
of $1.3\times 10^{-10}\ \rm W\ m^{-2}\ sr^{-1}$, which would imply $<5\%$ of
the 25.9 µm line surface brightness comes from [Fe II] 25.99 µm emission.
Therefore, [Fe II] contributes a negligible component of the 25.9 µm line and
we will assume for this work that the line flux is dominated by the [O IV]
line. The [O IV] line dominating the 25.9 µm blend is also consistent with
observations of blue compact dwarf galaxies, which have similar metallicities
and ionization parameters as N76 (Hunt et al., 2010).

Table 1: Line properties and intensities

|  | Ionization |  |  |  | $\mathrm{I_{obs}\ (10^{-15}\ W\ m^{-2})}$ | $\mathrm{I_{cloudy}\ (10^{-15}\ W\ m^{-2})}$
---|---|---|---|---|---|---|---
Line | $\mathrm{\lambda\ (\mu m)}$ | Potential (eV) | $\mathrm{n_{crit}\ (cm^{-3})}$ | FWHM ($\arcsec$) | Instrument | Total | $\mathrm{Orig}$ | $\mathrm{n_{H}=4\ cm^{-3}}$ | $\mathrm{n_{H}=10\ cm^{-3}}$
$\rm\left[SiII\right]$ | 34.80 | 8.2 | 2 $\times\ 10^{3}$ | 7.61 | LL1 | 14.75 $\pm$ 0.01 | $\ldots{}$ | 17.16 | 14.37
$\rm\left[SIII\right]$ | 18.71 | 23.3 | 2 $\times\ 10^{4}$ | 4.15 | LL2 | 11.28 $\pm$ 0.73 | $\ldots{}$ | 13.00 | 12.49
$\rm\left[SIII\right]$ | 33.47 | 23.3 | 7 $\times\ 10^{3}$ | 7.17 | LL1 | 18.93 $\pm$ 0.00 | $\ldots{}$ | 20.57 | 19.52
$\rm\left[SIV\right]$ | 10.51 | 34.8 | 5 $\times\ 10^{4}$ | 2.40 | SL1 | 17.58 $\pm$ 2.35 | $\ldots{}$ | 35.05 | 41.65
$\rm\left[OIII\right]$ | 88.33 | 35.1 | 5 $\times\ 10^{2}$ | 5.35 | PACS | 111.68 $\pm$ 2.25 | 50.21 $\pm$ 1.01 | 138.40 | 143.00
$\rm\left[NeIII\right]$ | 15.55 | 41.0 | 3 $\times\ 10^{5}$ | 3.51 | LL2 | 22.05 $\pm$ 0.70 | $\ldots{}$ | 28.33 | 28.67
$\rm\left[SI\right]$ | 25.25 | 0.0 | 1 $\times\ 10^{5}$ | 5.44 | LL1 | 1.04 $\pm$ 0.14 | $\ldots{}$ | 0.00 | 0.00
$\rm\left[OI\right]$ | 63.17 | 0.0 | 9 $\times\ 10^{5}$ | 3.63 | PACS | 15.74 $\pm$ 3.37 | 4.49 $\pm$ 0.90 | 1.30 | 1.02
$\rm\left[CII\right]$ | 157.64 | 11.3 | 5 $\times\ 10^{1}$ | 9.43 | PACS | 60.35 $\pm$ 1.20 | 16.70 $\pm$ 0.33 | 1.95 | 1.31
$\rm\left[NeII\right]$ | 12.81 | 21.6 | 7 $\times\ 10^{5}$ | 2.83 | SL1 | 4.52 $\pm$ 1.70 | $\ldots{}$ | 0.92 | 0.71
$\rm\left[OIV\right]$ | 25.88 | 54.9 | 1 $\times\ 10^{4}$ | 5.66 | LL1 | 20.35 $\pm$ 0.04 | $\ldots{}$ | 28.99 | 29.38
$\rm\left[NeV\right]$ | 24.32 | 97.1 | 1 $\times\ 10^{5}$ | 5.22 | LL1 | 0.20 $\pm$ 0.06 | $\ldots{}$ | 0.49 | 0.37

Note. — Description of the spectral lines modeled, including the integrated
intensities measured and the predicted integrated intensities with two
densities from Cloudy. The Ionization Potential is the photon energy required
to produce the observed ionization state from the lower ionization state. The
critical density is for collisions with electrons. The lines observed with
Herschel/PACS, and given as $I_{obs}$ Orig, do not cover the full nebula (see
Figure 1); $I_{obs}$ Total contains the total intensity corrected for this
missing flux.

Figure 3: Radial profiles of observed emission line brightness and Cloudy
photoionization models. The black circles represent the average line
brightness of a given emission line as a function of distance away from AB7,
the gray region is the standard deviation in a given bin, and the dashed line
represents the background diffuse contribution to the emission line (see
$\S$2.3). The Cloudy models vary with hydrogen density and are represented by
the colored lines and symbols, with the average of all models presented in
purple crosses. The ions that originate from the photoionized gas (excluding
emission lines that can be produced through neutral material ([C II] and [O
I]) see section 4.4, or the diffuse ionized gas, $\S$4.5) are well predicted
by the Cloudy models that range in density from $\rm n_{H}\sim 4\ cm^{-3}-10\
cm^{-3}$. Both the surface brightness and the shape of the radial profiles are
well predicted in these cases.

##  3 Photoionization models

###  3.1 Cloudy model parameters

We use version C17.01 of Cloudy (Ferland et al., 2017) to model the properties
of the ionized gas in N76. Cloudy is a spectral synthesis code designed to
simulate physical conditions and resulting spectrum from gas and dust that is
exposed to an ionizing radiation field. In order to produce a model, a variety
of input parameters are required and listed below.

Ionizing source SED: The shape of the source spectrum as a function of photon
energy must be specified in Cloudy, which for N76 corresponds to the spectral
energy distribution (SED) of the central binary. We use SED models in Shenar
et al. (2016) that constrain the stellar parameters of AB7 binary through a
combination of UV spectroscopy from the International Ultraviolet Explorer
(IUE) and the Far Ultraviolet Spectroscopic Explorer (FUSE) as well as optical
spectroscopy taken from Foellmi et al. (2003). The stellar parameters of the
WR and its O star companion were identified by matching the observed spectra
to models from the non-LTE Potsdam Wolf-Rayet (PoWR) model atmosphere
code111http://www.astro.physik.uni-potsdam.de/PoWR.html (Gräfener et al.,
2002; Hamann & Gräfener, 2003; Sander et al., 2015; Hamann & Gräfener, 2004;
Hainich et al., 2019). These PoWR models are set to use the SMC metallicity
and derive a variety of parameters for the WR and O star. Shenar et al. (2016)
report that the WR is a WN4 spectral type with a temperature of
$T_{*}=105000\rm K$, wind speed of $v_{\infty}=1,700$ km s-1, surface gravity
of log $g_{*}=4.7\ \rm cm\ s^{-2}$, a low ($X_{\rm H}\approx 0.15$) hydrogen
fraction, a transformed radius of log $R_{t}=0.75\rm\ R_{\odot}$, and a mass
loss rate of log $\dot{M}=-5.0\rm\ M_{\odot}\ yr^{-1}$ . The O star companion
is classified as an O6 I(f) star and has a temperature of $T_{*}=36$ kK, wind
speed of $v_{\infty}=1,500$ km s-1, surface gravity of log $g=3.6$ $\rm cm\
s^{-2}$, and a mass loss rate of log $\dot{M}=-7.0\rm\ M_{\odot}\ yr^{-1}$.

The PoWR models are normalized by default to a luminosity of $\rm
L_{\odot}=5.3\mbox{\,$\rm L_{\odot}$}$. Therefore, we scale the SEDs fitted in
Shenar et al. (2016) by their luminosities of $\rm log\ L_{\odot}=6.1$ for the
WN4 and $\rm log\ L_{\odot}=5.5$ for the O6 components. We then add the
luminosity weighted WN and O6 SEDs together to produce a total SED for the AB7
binary system, shown in Figure 4. The feature at 40 eV in the spectrum is due
to a strong FeV emission line at 305.31 Å.

Figure 4: The SED of the AB7 system for Cloudy input. We use the Potsdam Wolf-
Rayet database models for the WN Wolf-Rayet (Hamann & Gräfener, 2004; Todt et
al., 2015) and O6 (Hainich et al., 2019) stars in the AB7 binary and combine
them into one input SED. The full SED for AB7 is in blue, the WR is green, and
the O star is orange.

Elemental abundances: Instead of scaling gas-phase Milky Way abundances to the
metallicity of SMC, we employ the best gas-phase abundances available for SMC
HII regions. Relative abundances of elements are different in the Magellanic
Clouds and the Milky Way due to different nucleosynthesis histories (Russell &
Dopita, 1992). For the majority of elements, we adopt the abundances given for
HII regions provided in Dufour (1984) and Russell & Dopita (1992), but we
slightly modify the abundances for sulphur and silicon to best match the
observations in N76 (see below). Our adopted abundances and the references
used are in Table 2.

We decrease the sulphur abundance by a factor of two from the
$\rm(S/H)=3.9\times 10^{-6}$ value reported in Russell & Dopita (1992).
Without this modification, the three sulphur lines observed in this work ([S
III] 18 µm, [S III] 33 µm and [S IV] at 10 µm) are consistently over-predicted
by a factor of two in our Cloudy models. This small variation in the sulphur
abundance is also consistent with the scatter presented in Russell & Dopita
(1990), the value of $\rm(S/H)=1.8\times 10^{-6}$ reported in Dennefeld &
Stasińska (1983), and the sulphur abundances found in Vermeij & van der Hulst
(2002) and Lebouteiller et al. (2008) that use infrared data.

Identifying an abundance for silicon is more complex because silicon is often
depleted in dust. While the dust abundance is often lower in HII regions than
HI regions or molecular clouds, dust can still be present and have an impact
on HII regions, especially towards the outer edges where the ionized hydrogen
fraction is lower and the [Si II] emission is brightest. Russell & Dopita
(1992) provide only a photospheric abundance for silicon that does not take
into account the depletion from dust. We therefore use elemental depletion
studies in the SMC to calculate the median gas phase abundance for silicon
(Tchernyshyov et al., 2015; Jenkins & Wallerstein, 2017). The median silicon
depletion from these studies is 0.4 which happens to accurately match the
observed flux of [Si II]. We note that this systematic uncertainty in the
silicon HII region abundance leads to an uncertainty in the prediction for the
[Si II] flux (see discussion in $\S$4.4).

Table 2: Adopted SMC Abundances

Element | Abundance (X/H) | Reference
---|---|---
He | 8.13 $\times\ 10^{-2}$ | 2
C | 1.45 $\times\ 10^{-5}$ | 1
N | 4.27 $\times\ 10^{-6}$ | 2
O | 1.07 $\times\ 10^{-4}$ | 2
Ne | 1.86 $\times\ 10^{-5}$ | 2
Al | 2.51 $\times\ 10^{-6}$ | 2
Si | 4.88 $\times\ 10^{-6}$ | 3,4
S | 1.95 $\times\ 10^{-6}$ | 2

Note. — References are (1) Dufour (1984) (2) Russell & Dopita (1992), (3)
Tchernyshyov et al. (2015), and (4) Jenkins & Wallerstein (2017).

Stopping criteria: We are interested in modeling the ionized gas properties of
N76 and therefore will run the Cloudy model to the end of the ionization zone
at the transition between the HII region and PDR. We therefore define the HII
region end to be the point at which the PDR begins and the gas is neutral. The
focus of this work is to model the ionized gas so continuing the model into
the PDR is unnecessary, but has the consequences of underestimating the
emission from the neutral dominated lines [Si II], [C II], and [O I] (see
$\S$4.4). After using a range of stopping conditions we find that an electron
temperature cutoff of 2000 K accurately probes the ionized gas to the end of
the ionization front.

Geometry: We assume a spherical geometry to approximate the geometry of N76.
The sphere command in Cloudy sets the geometry to be closed and assumes that
the gas fully covers the ionization source.

Inner radius: The inner radius is defined as the distance between the gas and
the ionization source. The massive, fast stellar winds produced from Wolf-
Rayet stars often form a stellar wind-blown bubble (SWB) structure in the
resulting nebula (e.g., Weaver et al., 1977). As we discuss in §4.3, the
region closest to AB7 has been shocked by the stellar winds and is occupied by
this tenuous very hot plasma. For the purposes of our modeling, this effect is
best reproduced by a central cavity in the gas. We vary the gas inner radius
and find a radius of 10 pc to best reproduce the observed emission lines
radial profiles.

Hydrogen density: The density input in Cloudy is the total hydrogen density.
We assume the density of the HII region is constant, a reasonable assumption
for classical, spherical HII regions, and the density remains constant up to
the ionization front of the region. We produce multiple constant densities
models that vary from $\rm n_{H}=1-50\ \mbox{$\mathrm{cm^{-3}}$}$ in order to
best constrain the density in N76 with the observed emission lines. The range
of hydrogen densities between $\rm n_{H}=4-10\ \mbox{$\mathrm{cm^{-3}}$}$ that
fit the data best are presented in this work.

Since the binary at the center of N76 is well-studied, most of the inputs in
the Cloudy model are using a priori knowledge of the ionization source,
abundances, and geometry.

###  3.2 Calculating predicted radial profiles

Figure 5: Schematic showing the geometry of calculating the projected
intensity from Cloudy models where $d$ is the projected distance, $r$ is the
physical distance, and $dr$ is the differential distance between shells. The
observer is peering from the right of the figure.

N76 is an approximately spherical shell, also known as a bubble nebula, with a
radius of about 40 pc. Because of the approximately symmetric nature of the
nebula, we compare the Cloudy models and the observations using radial
profiles. Cloudy models report the volume emissivity of each transition as a
function of depth into the nebula. In other words, they provide the radial
profile of emissivity in a sphere. To compare with the observations, we need
to appropriately project these profiles to compute the azimuthal average
surface brightness as a function of projected distance from AB7 for each of
the emission lines. A point that is at a distance $d$ from the source in
projection has contributions from shells that have radii between $d$ and
infinity in proportion to their emissivity and projected thickness. Figure 5
shows the geometry of the N76 Cloudy model and the relationship between
different shells, the projected distance $d$, and the physical distance $r$.
The projected surface brightness for the Cloudy model then is:

$I(d)=\frac{1}{2\pi}\int_{d}^{\infty}\frac{\epsilon(r)}{\sin\theta}dr\ ,$ (1)

where $r$ is the physical 3D distance from AB7, $\epsilon(r)$ is the radial
dependence of the line volume emissivity reported from Cloudy, $d$ is
projected distance from AB7, and $\theta$ is the angle between $d$ and $r$. We
convert Equation 1 to eliminate $\theta$ resulting in

$I(d)=\frac{1}{2\pi}\int_{d}^{\infty}\frac{\epsilon(r)}{\sqrt{1-d^{2}/r^{2}}}dr\
,$ (2)

and calculate the surface brightness predicted from the Cloudy model, $I(d)$,
as a function of the projected distance and volume emissivity reported from
Cloudy in order to compare to the observed data. We set the filling factor in
Cloudy to unity to reflect that the resolved emission in N76 is not clumpy.

The Cloudy models use AB7 as the only ionization source, but there may be
contributions of low ionization ionized gas from the warm ionized medium to
our observed emission lines (e.g., Haffner et al., 2009). To account for this,
we add the constant surface brightness background value calculated in $\S$2.3
to the radial profiles derived from the Cloudy models. This background is
typically very low representing $\lesssim 5\%$ of the total surface
brightness. Lastly, the final Cloudy profiles are convolved with a Gaussian
kernel with a FWHM of 12″ to make an equal comparison to the 12″ PSF matched
emission line images.

##  4 Results

###  4.1 Comparison of models to observations

Figure 3 shows the results of the constant density Cloudy model with hydrogen
densities $\rm n_{H}=4-10$ cm-3. We compare the models projected on the plane
of the sky (§3.2) with the measured emission line radial profiles (§2.5).
Therefore we compare not only the integrated intensity of the emission line,
but also how well its predicted distribution matches the data. The emission
lines we sample are at a variety of ionization energies and therefore trace
the different ionization zones in the nebula. For example, the [O IV] line
($\rm E_{ion}=55\ eV$) peaks at 10 pc while the [S III] line ($\rm E_{ion}=23\
eV$) peaks at 20 pc. A model that predicts an emission line well will have a
similar peak, spatial gradient, and ending position as the observed radial
profile

There is not a single constant density model that is able to predict the shape
and intensity of all of the emission lines, but most emission lines are
consistent with any given constant density model in the range of $\rm n_{H}=4$
cm-3 to $\rm n_{H}=10$ cm-3. In some cases, the spatial distribution of the
profile is better matched by a model that does not predict the intensity as
well. For example, in the [S IV] profile, the 8 cm-3 model predicts the shape
very well but is about 50% more luminous than the observations. These small
discrepancies are likely due to our constant density and spherical symmetry
simplifying assumptions. While there are small differences between the models
and observed emission line profiles, overall Cloudy predicts the emission from
ionized gas very well with simple models in a narrow density range, except for
a few transitions that we will discuss in detail.

In addition to the projected surface brightness profiles from the Cloudy
models, we also report the predicted total intensities in the Cloudy models in
Table 1 for the $4\ \mbox{$\mathrm{cm^{-3}}$}$ and $10\
\mbox{$\mathrm{cm^{-3}}$}$ models. Overall the radial profiles of the surface
brightness provide a much more precise way of differentiating between the
constant density Cloudy models than examining the integrated intensity. For
example, the integrated intensity of [Ne III] varies between the $4\
\mbox{$\mathrm{cm^{-3}}$}$ and $10\ \mbox{$\mathrm{cm^{-3}}$}$ models by 28.33
to 29.38 $\rm W\ m^{-2}$, a $<1\%$ difference. In contrast, the radial profile
of the $4\ \mbox{$\mathrm{cm^{-3}}$}$ model peaks 10 pc after the $10\
\mbox{$\mathrm{cm^{-3}}$}$ does, clearly showing it does not match the spatial
distribution of the [Ne III] emission. Thus a comparison between the spatial
resolved radial profiles provides much more information than using just
integrated intensities.

To produce both the Cloudy models and the measured radial profiles we assume
that N76 is spherically symmetric with AB7 in the center of the nebula. This
is approximately true, but not precisely correct. Figure 1 shows that AB7
appears offset in the eastern direction. Additionally, there is an overall
trend for the ionized gas lines to be much brighter on the eastern edge of the
nebula in comparison to some of the neutral gas lines ([CII] especially) that
have brighter emission on the western edge. We explore the effects of the
asymmetries on our results in Appendix B. Images of the H2 S(1) transition and
ancillary ALMA CO data (Figure 7) confirm that there is a molecular cloud
complex on the western side of N76. Therefore, there are likely internal
density gradients as well as a wall of dense gas on the western side of the
nebula breaking the symmetry assumption. For example, because AB7 is slightly
offset from the center, the radial averaging produces flatter slopes with
respect to the models. We have attempted to account for some of this offset by
showing the profiles as a strip instead of radially in Appendix B. These show
the structure in the different sides of the nebula and how a particular
density is a better predictor for one side compared to the other.

We also report the ionization structure predicted by Cloudy of each element in
Figure 6 to complement the projected profiles of surface brightness. The
panels show the ionization state of the relevant elements as a function of the
physical distance from AB7. In addition to the ions that are studied in this
paper, we also plot the ionization structure of helium. As the most abundant
element after hydrogen, helium can have a significant impact on the ionization
equilibrium of HII regions. This impact is most pronounced in HII regions
illuminated by hard radiation fields, due to the fact that their central
sources produce a lot of photons capable of ionizing and doubly ionizing
helium: the ionization potentials of $\rm He^{0}$ and $\rm He^{+}$ are 24.6 eV
and 54.4 eV respectively. Ionization cross-sections are sharply peaked at the
ionization energy: the cross-section of $\rm He^{0}$ to 24.6 eV photons is ten
times higher than that of $\rm H^{0}$, therefore these photons will
preferentially ionize helium over hydrogen. For sources with a hard enough
spectrum, characterized by the ratio of the rate of production of He-ionizing
photons ($Q_{1}$) to H-ionizing photons ($Q_{0}$), this can lead to the $\rm
He^{+}$ region extending slightly past the $\rm H^{+}$ region. In fact,
$Q_{1}/Q_{0}\sim 0.15$ is needed for the HII and HeII regions to be equivalent
sizes (Draine, 2011), while for N76 the $Q_{1}/Q_{0}\sim 0.55$ calculated from
the Cloudy models matches the slightly larger HeII emitting region seen in
Nazé et al. (2003). A similar effect is seen in the ionization structure of
neon, where the [Ne II] emission extends slightly past the hydrogen ionization
front.

Figure 6: The ionization structure of sulphur, oxygen, neon, silicon, carbon,
and helium from the 10 $\mathrm{cm^{-3}}$ constant density Cloudy model
beginning at the photoionized region (r = 10 pc). Ions that are presented in
this paper are bold, and we also show the structure of ionized hydrogen in
black. The various ionization zones can be seen throughout these elements in
physical distance, similar to the projected distance seen in Figure 1.

###  4.2 The high ionization lines: [Ne V] and [O IV]

The high ionization lines, [Ne V] and [O IV], with ionization potentials of 97
eV and 55 eV respectively, can often be difficult to reproduce through just
photoionization from O and B stars. Models of hot, low metallicity WR stars
around ages of 3-4 Myr like the WR in AB7, however, contain enough high energy
photons to ionize Ne+4 and O+3 and produce appreciable [Ne V] 24 µm and [O IV]
25 µm emission (Schaerer & Vacca, 1998; Schaerer & Stasińska, 1999; Crowther
et al., 1999). The N76 photoionization models presented here are able to
accurately predict both the flux of the [Ne V] and [O IV] lines (see Table 1)
as well as their spatial distribution (see Figure 3). The [Ne V] is slightly
over-predicted by the Cloudy models by a factor of two, which can be explained
by density variation (either through clumping or a lower density in the high
ionization zone) or uncertainty in the PoWR model (see below).

We note that the predicted flux of the [Ne V] and [O IV] lines are highly
dependent on small changes to the PoWR model that impact the high energy (55+
eV) tail of the SED. A Cloudy model using a PoWR SED with the same input
parameters but increasing the hydrogen fraction from $X_{H}\sim 0.15$ to
$X_{H}\sim 0.2$ produced no [Ne V] emission and a fraction ($<10\%$) of the
observed [O IV] emission. Additionally, the O star companion to the WR in AB7
contains a small X-ray contribution to the SED, implied by the presence of OVI
in its spectrum, but contributes a negligible ($<5\%$) amount to the [O IV]
and [Ne V] flux (see Shenar et al., 2016, for more details).

Further, Nazé et al. (2003) report narrow-band optical imaging of N76 and find
bright HeII $4686\ \rm\AA$ emission that is inconsistent with the hottest WR
models at the time. The Cloudy models presented here, however, are able to
reproduce the reported HeII emission within a factor of two. The spatial
distribution of the HeII emission (see Figure 4 in Nazé et al., 2003) is very
similar to the [O IV] emission, likely due to their similar ionization
potentials.

###  4.3 The wind-blown cavity

The Cloudy photoionization models require a 10 pc cavity to reproduce the
observed spatial distribution of the ionized gas lines. This origin of this
cavity is from the stellar winds of the WR and O star in AB7 interacting,
shocking, and heating the surrounding ISM, creating a stellar wind-blown
bubble (SWB). The resulting nebula consists of an onion like structure (see
Figure 1 in Weaver et al., 1977; Freyer et al., 2003) with four layers,
including (1) an innermost, free-flowing supersonic wind that travels into (2)
a hot ($\sim 10^{6}-10^{8}$ K), shocked gas region that will expand into (3)
the photoionized shell of swept up interstellar gas, and finally into (4) the
ambient interstellar medium. The hot, shocked material in shell 2 is typically
at much higher pressures than shell 3, and will expand as a function of time.
Assuming a uniform density spherical shell with equal pressure throughout and
that the thermal energy contained within shell 2 is much higher than the
kinetic energy, Weaver et al. (1977) derived a relationship between the radius
of hot, shocked gas bubble (shell 2) as a function of time, such that

$r_{s2}(t)=\left(\frac{125}{154\pi}\right)^{1/5}L_{w}^{1/5}\rho_{0}^{-1/5}t^{3/5}\
\rm{cm},$ (3)

where $\rho_{0}=n_{0}m_{H}$ is the density of the ambient medium in $\rm g\
cm^{-3}$, $L_{w}=\frac{1}{2}\dot{M_{w}}v_{w}^{2}$ is the stellar wind
luminosity in $\rm g\ cm^{2}\ s^{-3}$ with $\dot{M_{w}}$ as the mass loss rate
and $v_{w}$ as the terminal velocity of the stellar wind, and $t$ is the time
of expansion in seconds. We calculate the size of the bubble blown by the WR
star in AB7 through the values characterized by Shenar et al. (2016),
$\dot{M}=10^{-5}\ \mathrm{\mbox{\,M${}_{\odot}$}\ yr^{-1}}$, $v_{w}=1.7\times
10^{3}\ \mathrm{km\ s^{-1}}$, and assume $n_{0}=10\ \mathrm{cm^{-3}}$. Using
these parameters, a 10 pc radius bubble is created after 0.2 Myr of free
expansion. This is consistent with the evolutionary tracks studied in Shenar
et al. (2016), where they find an age of 3.4 Myr for the AB7 system and
estimate the WR star has been in the WR stage for $\sim$0.1-0.2 Myr.

Another contribution to the cavity in N76 may be through the colliding stellar
winds from the WR and O star. There is ample evidence of shocks originating
from colliding wind systems, such as the winds produced by the WN4 and O6 star
in AB7 (e.g., Stevens et al., 1992; Tuthill et al., 1999; Parkin & Pittard,
2008). However, the spatial scale of the colliding wind zone is comparably
insignificant to size of the photoionized cavity, corresponding to $\sim
10^{-6}$ pc for AB7. Further, Chandra observations of AB7 show X-ray emitting
gas attributed to shocks from the colliding winds, but the region is only
$\sim$3.5 pc in radius (Guerrero & Chu, 2008). Thus the SWB structure
described in Weaver et al. (1977) likely creates the implied cavity in the
presented photoionization models.

It is possible that shocks from either the colliding winds of the AB7 system
or the stellar wind blown bubble could contribute additional sources of [Ne V]
and [O IV] in N76. We investigate this by calculating the velocity and
postshock temperature needed to ionize Ne+3 into Ne+4. We estimate a
temperature of $\rm T\sim E_{ion}/k_{B}\sim 10^{6}$ K by taking the ionization
energy of [Ne V]. Using the relationship between the postshock temperature
(equation 36.28 of Draine, 2011) and the shock velocity, we find the minimum
velocity for the shock is $\rm v\approx 290$ km s-1. Toalá et al. (2017) study
the Wolf-Rayet Nebula NGC 3199 in the Milky Way and find diffuse X-ray
emission that corresponds to a plasma temperature of $\rm T\approx 1.2\times
10^{6}\ K$, very similar to the temperature required to produce the [Ne V]
emission observed in N76. Further, models of stellar wind blown bubble show
velocities of $\sim 200-400$ km s-1, enough to ionize Ne+3 into Ne+4 (Freyer
et al., 2003, 2006; Lancaster et al., 2021). Since the Cloudy photoionization
models already slightly overpredict the [Ne V] emission, it is unlikely there
is this additional contribution from shocks. The cavity may contain ions of
neon at much high energy states not detected by our current dataset.

Fast, radiative shocks with velocities $\sim$300-500 km s-1 can explain
optical observations of [Ne V] and HeII in low metallicity, blue compact dwarf
galaxies (e.g., Thuan & Izotov, 2005; Izotov et al., 2012, 2021). However, the
origin of these shocks is typically attributed to supernovae (SN), not stellar
winds. Indeed, the supernova remnant E0102-72 that serendipitously is in the
N76 field is a bright [Ne V] source (see Figure 1), suggesting shocks from SN
can produce appreciable [Ne V] emission. The SNR is too far away from N76 to
contribute to the [Ne V] emission associated with the nebula (see Figure 1).
Therefore, the [Ne V] within the N76 nebula is most likely produced through
photoionization. We note that if N76 was placed at the larger distances of the
blue compact dwarf galaxies studied in Izotov et al. (2021), it would be
difficult to differentiate between [Ne V] produced through the Wolf-Rayet
nebula and the SNR due to the lower physical resolutions.

###  4.4 The neutral gas dominated lines: [O I], [C II], and [Si II]

We define the neutral-dominated lines as ions or atoms that can come from
neutral gas, i.e. those ions that have an ionization potential less than that
of hydrogen. Most of these neutral-dominated emission lines are under-
predicted by the Cloudy models because the Cloudy models were designed to
model the ionized gas and do not probe farther than the ionization front after
the transition from ionized to neutral material. Further, Figure 6 shows the
ionization structure of carbon and oxygen, suggesting that the radiation field
of AB7 is too hard for $\rm O^{0}$ and $\rm C^{+}$ to exist in the interior of
the HII region where oxygen and carbon are in higher ionization states. Thus
the bulk of emission from [O I] and [C II] comes from the neutral phases;
future photo-dissociation region (PDR) modeling is needed to understand the
nature of these neutral lines, which is beyond the scope of this study.

The fine-structure transition of silicon, [Si II], however, is reasonably well
predicted by the Cloudy models. With an ionization potential of 8.15 eV for
$\rm Si^{0}$, [Si II] can come from the neutral gas. However, Cloudy modeling
of the HII region reproduces both the total intensity and the surface
brightness distribution of the [Si II] line (see Table 1 and Figure 3), unlike
for the other neutral gas tracers. Figure 6 shows that the modeled abundance
of Si+ grows faster toward the edge of the HII region than that of C+, which
is consistent with the observed radial profiles presented in Figure 3. Thus
the majority of the [Si II] emission is associated with the edge of the HII
region, around the transition between ionized and atomic gas. In subsection
5.1 we calculate the percentage [Si II] comes from each phase.

Why are there dominant PDR contributions for [C II] and [O I] but not [Si II]?
The [Si II] 34 µm fine-structure transition is at shorter wavelengths, and its
excitation requires higher temperatures than [O I] 63 µm or [C II] 158 µm.
Further, the [C II] and [O I] lines are also major coolants in the Warm and
Cold Neutral Medium (Wolfire et al., 1995, 2003) and contributions from that
can increase the total intensity of these lines by peering through the atomic
gas in projection towards the SMC. In contrast, models presented in Wolfire et
al. (1995, 2003) show the cooling from [Si II] in the atomic gas is two orders
of magnitude lower than that from the [C II] and [O I] lines, which leads to a
much fainter atomic gas contribution for the [Si II] emission. The intensity
of the [Si II] emission is also highly dependent on the silicon gas abundance.
Silicon is depleted heavily into dust, and in order to determine its gas phase
abundance we used depletion studies in the SMC (Tchernyshyov et al., 2015;
Jenkins & Wallerstein, 2017). The dust depletion of gas phase silicon with
respect to the photospheric abundances in the SMC (Russell & Dopita, 1992) has
a large range, factors from 1.25 to 0.05 (0.1 to -1.3 dex, our preferred
abundance corresponds to 0.4, or -0.4 dex). The assumed silicon abundance
matches well the observed radial profile in in Figure 3. A lower silicon
abundance leads to a [Si II] radial profile that does not match the overall
shape and would suggest a larger fraction of [Si II] comes from the PDR. While
we assume some silicon depletion in the ionized gas, our results suggest that
even more silicon is depleted in the neutral phase, especially when compared
to carbon or oxygen, which will decrease the available gas-phase silicon in
order to produce the [Si II] emission. These combined effects likely explain
the much smaller contribution of PDRs toward the observed [Si II] emission.

We present images of the other neutral gas tracers in N76 in Figure 7. The
Spitzer IRS covers the quadrupole transitions of molecular hydrogen, H2. Our
PAHFIT pipeline simultaneously fits these lines and we show the brightest H2
line, S(1) at 17 $\mu$m, in Figure 7. This line has $E/k\sim 1000$ K, traces
H2 at $T>300$ K, and its presence indicates molecular material that is heated
by a strong FUV field or through shocks from the central WR star. We compare
the H2 S(1) and [Si II] lines to the ACA CO(2-1) data presented in Tokuda et
al. (2021). Interestingly, these neutral gas species are not always spatially
coincident. As discussed above, the [Si II] emission originates mostly around
the ionization front of the HII region and therefore traces the the outer
shell of N76. The H2 and CO emission is brightest where AB7 illuminates the
molecular ridge region that is north and north-west of N76.

Figure 7: The neutral gas tracers in N76, with the H2 S(1) quadrupole
transition in the blue color scale, CO contours from Tokuda et al. (2021) in
black (in levels of 1, 6, 11, 16 Jy/beam), and the [Si II] emission in magenta
contours (in levels of 1, 1.5, and 3 $\rm\times\ 10^{-8}\ W\ m^{-2}\
sr^{-1}$). The beam/PSF size of each map is in the bottom left corner. The [Si
II] emission primarily traces the outer shell of N76, suggesting [Si II]
originates mostly from the transition between the ionized gas and the PDR. The
CO and H2 emission is brightest when AB7 illuminates the molecular ridge north
and north-west of N76.

###  4.5 Production of [Ne II] in N76

The Cloudy models predict a modest amount of [Ne II], about 20% of the
observed [Ne II]. At an ionization potential of 21.56 eV, [Ne II] traces lower
energy photons. Figure 6 shows how the ionization state of Neon transitions to
$\rm Ne^{+}$ at the very edge of N76 at a very small radius range compared to
the higher ionization states and rapidly transitions into neutral Neon.
Therefore, the total volume of the predicted [Ne II] emission from the Cloudy
models is quite small, leading to a very weak predicted [Ne II] line. The
excess of observed [Ne II] emission could be associated with other ionization
sources besides AB7, such as smaller HII regions. It is also possible that [Ne
II] traces the low density, diffuse ionized gas at densities $n_{e}\lesssim
0.1\ \rm cm^{-3}$. In the nearby dwarf galaxy IC 10, Polles et al. (2019) also
find that their models underpredict the [Ne II] observations by a factor
between three and four, suggesting that [Ne II] may be a tracer of the diffuse
ionized gas in low metallicity galaxies.

###  4.6 How changes in the model affect the inferred conditions

Besides the density, for which we produce a grid of models, a variety of other
Cloudy inputs can affect the final results of the modeling. The Cloudy models
are most sensitive to the input luminosity and SED.

For example, decreasing the luminosity of the WR by 50% leads to a 56%
reduction in the total intensity of the spectral lines at a constant density
of $n_{H}=10\ \mbox{$\mathrm{cm^{-3}}$}$. Naively, this reduction in the total
intensity provides a slightly better fit for some lines, but it is at the cost
of producing emission over a region with a smaller radius and a spectral line
radial profile that does not match the data. Thus the strength of comparing
the Cloudy models to the spatially resolved data shines through; using the a
priori values for the spectrum and luminosity provide an overall better fit
than if one were to directly match the total intensities.

We also experimented with using different stopping conditions, abundances, and
the presence of dust. All of these parameters have a negligible impact on the
overall results of the model, except for the abundances which we fixed in the
manner described in $\S$3.1 and Table 2.

###  4.7 Morphology of emission lines

Below we discuss some features that are apparent in our images of the N76
region shown in Figure 1.

Sharp and bright eastern edge. In most of the ionized gas emission lines,
there is a brightness gradient with a northeastern sharp edge compared to a
fainter and more diffuse southwestern edge. This is especially prominent in
the lower ionization lines such as [S III]. We test how the brightness
gradient impacts the Cloudy modeling in Appendix B. Modeling each individual
side suggests a density difference of order $\pm 2$ $\mathrm{cm^{-3}}$ from
the best density in the radially averaged models, or a typical density
difference between each side of approximately $20\%-30\%$.

Potential Spiral structure. In some of the high ionization lines, particularly
[O IV] and [S IV], there appears to be a “pinwheel” spiral structure emanating
from AB7. The spiral has a single arm beginning at the location of AB7 and
coils eastward then westward, ending about 10 pc away from AB7. This structure
is also present in narrow-band optical line imaging of HeII $4686\ \rm\AA$
from Nazé et al. (2003). Spirals in WR-O star binaries have been seen before
in dust emission (e.g., Tuthill et al., 1999; Lau et al., 2022) caused by the
interacting stellar winds and the orbits of the binary. Spiral features are
typically detected in the dust continuum of carbon-rich WC stars, but a spiral
structure was reported in a nitrogen-rich star similar to the Wolf-Rayet in
AB7 in radio continuum likely due to synchrotron emission (Rodríguez et al.,
2020). This suggests that although the presence of dust helps illuminate a
spiral pattern, it is intrinsic to binary systems and also present in other
tracers. The spiral features documented in the literature, however, are on the
scale of hundreds or thousands of AUs rather than the $\sim$10 pc feature we
see in N76. The wind-collision zone where a spiral structure is expected for
AB7 is $\sim$2 AU in size, and therefore unresolved in our data. Nevertheless,
it is possible that the high ionization lines of N76 show an imprint of the
spiral motions from past stellar wind collisions in the central region, which
then expand with the wind. Alternatively, the suggestive spiral structure may
be coincidental and due to density gradients throughout the photoionized gas.

Dense knot. In the southeastern edge of N76 there is a bright knot of emission
seen across the Spitzer IRS low ionization emission lines (and especially
prominent in the [Si II] line). Previous work studying N76 suggest that this
dense region is an independent compact HII region occasionally referred to as
a “high excitation blob” (Heydari-Malayeri et al., 2001; Nazé et al., 2003).
These are very dense, small regions that are associated with young massive
stars beginning to leave their natal molecular cloud. Figure 7 shows the
highly excited molecular gas in this blob, both in the H2 quadrupole
transitions as high as S(7) and a small clump of CO (2-1) (Tokuda et al.,
2021). There are two possible ionizing sources for this compact region,
including a Be star (Wisniewski & Bjorkman, 2006) and a candidate Young
Stellar Object (Bolatto et al., 2007). It is unclear whether either of these
stars can power the ionization seen in this clump. Alternatively, it is
possible that this feature is simply a dense blob of gas irradiated by AB7 and
shocked by the HII region expansion.

###  4.8 Tentative [S I] detection in the center of N76

One peculiar detection from Spitzer IRS is a spectral line at 25.25 µm in the
center of N76 that we tentatively identify as [S I]. The spectrum of this line
can be seen in Figure 2 and the image in Figure 1. Our fitting procedure from
PAHFIT identifies a spectral line between 25.2 and 25.3 µm throughout the
central region of N76. The neutral sulphur line is directly in the middle of
this range, at 25.245 µm. Using the NIST atomic spectral database, the nearest
possible spectral lines correspond to Mg II at 25.05 µm and Si I at 25.38 µm.
The redshift of these wavelengths due to the systemic velocity of the SMC is
$0.01$ µm. The average uncertainty PAHFIT reports on this spectral line
wherever there are detections is 0.017 µm, which corresponds to $\sim 10\%$ of
the wavelength step between spectral pixels in this region of the LL1 map
(0.178 µm).

The existence of neutral sulphur (or for that matter any line corresponding to
ions found in neutral gas, such as Mg II and Si I) in the center of such a
high ionization HII region is a mystery since neutral sulphur should be
ionized quickly nearby such a hot and luminous star. Indeed, Cloudy predicts a
total of $0.0005\ \times 10^{-15}\rm\ W\ m^{-2}$ for the total [S I]
intensity, which is almost four orders of magnitude fainter than the observed
total flux of $1.04\times\ 10^{-15}\rm\ W\ m^{-2}$. The radial profile Cloudy
predictions are also many orders of magnitude fainter than the observed flux,
predicting an average of $\sim 10^{-13}\ \rm W\ m^{-2}\ sr^{-1}$ compared to
the observed $\sim 10^{-9}\ \rm W\ m^{-2}\ sr^{-1}$ radial profile. If the
line is confirmed as [S I], then it is the first detection of such a feature
in the center of the highly irradiated environment of a WR nebula. Previous
observations show [S I] resides in dense, neutral gas and can be produced
through fast, dissociative J-shocks (Haas et al., 1991; Rosenthal et al.,
2000). These observations, however, consistently show the presence of other
neutral lines, such as H2 and [O I]. In N76, there is no corresponding neutral
emission where we observe [S I], making the production of this line a mystery.
We explore possible explanations for this emission in $\S$5.2.

##  5 Discussion

###  5.1 Ionized gas contributions to the neutral lines: [Si II], [C II], and
[O I]

The ions and atoms that can exist in the neutral gas, $\rm Si^{+}$ and $\rm
C^{+}$, have ionization potentials less than that of hydrogen, but their next
ionization stage is above or close to 13.6 eV. The emission from [Si II] and
[C II] can therefore be present in both the neutral and ionized gas phases.
For [O I], the ionization potential of oxygen (13.618 eV) is very close to
that of hydrogen and facilitates a very efficient charge-exchange reaction,
ensuring that there is very little [O I] emission associated with the ionized
gas. Although this work is focused on modeling only the ionized gas,
understanding the proportion the ionized gas contributes to the overall
intensity of the neutral-dominated lines is of interest for using these lines
as neutral gas tracers. Our Cloudy simulations were designed to end at the
ionization front, the position where the ionized gas transitions into fully
neutral gas. To match the method in Cormier et al. (2019), we compute the
contribution the neutral gas has on [Si II] and [C II] by defining a more
conservative “end” of the ionized gas region as the point in which the
electron fraction and proton fraction are equal to 0.5. We use this value to
define the end of the ionized gas region.

As described in $\S$4.4, the Cloudy models underpredict the emission from [O
I] and [C II] that is observed in the direction of the nebula, because the
neutral gas regions that surround the ionized gas but are not included in our
Cloudy models contribute the bulk of the emission in these lines. In contrast,
the [Si II] intensity is well matched to the predicted emission from ionized
gas in the Cloudy models, suggesting a substantial amount of [Si II] is
produced at the edge of the HII region itself and therefore not outside the
HII region in the PDR. We quantify this by calculating the fraction of the
intensity originating from the ionized gas for these lines in our models by
defining the end of the ionized gas region as the point where the electron and
neutral fraction are equal to 0.5. We find that for N76 91% of the [Si II]
emission originates from the ionized gas while the [C II] has an ionized gas
contribution of 3%. Cormier et al. (2019) use their Cloudy models of the Dwarf
Galaxy Survey (DGS) to perform a similar calculation and report that
$\lesssim$10% of the [C II] and $\lesssim$40% of [Si II] emission originates
from the HII regions. Although we see a much higher HII region contribution to
the [Si II] emission in N76 than in the galaxy-wide integrated measurements in
the DGS, keep in mind that N76 is a Wolf-Rayet nebula and thus not
representative of most HII regions.

###  5.2 Production of [S I] in the center of N76

We discussed in $\S$4.8 the surprising identification of a line in the central
regions of N76 as neutral sulphur. Below we go over some possible ideas about
the origin of this transition.

Velocity shifted [O IV] emission: The nearest emission line in wavelength to
our possible [S I] emission that is bright and clearly identified in N76 is [O
IV] at 25.93 µm. We considered whether the line we identify as [S I] could be
high velocity [O IV] emission. The velocity required to shift [O IV] to the
observed wavelength of 25.25 µm is 7,600 km s-1. While the WR and O stars in
AB7 produce stellar winds, the speed of these winds is below 2,000 km s-1
(Shenar et al., 2016). The X-ray emission observed in N76, resulting from the
colliding WR and O star winds, also implies a wind speed of 2,000 km s-1
(Guerrero & Chu, 2008). Thus it is extremely unlikely that the emission is [O
IV] shifted to 25.25 µm. Further, the wind would have to be asymmetric, since
we do not see a corresponding redshifted component. Moreover, instead of
discrete components we would expect a symmetric, broadened [O IV] profile
because the emission should be optically thin.

Dust destruction through shocks from stellar winds: The main mystery of the
origin of [S I] at the center of this nebula is how sulphur can remain in a
neutral state where ionized species caused by $>$55 eV photons are also
observed. It is possible that the winds and radiation from AB7 can shock the
surrounding dust, destroy dust grains, and move their constituent atoms,
including neutral sulphur, into the gas phase. We then observe the [S I]
emission before neutral sulphur is ionized into higher states. Most previous
detections of [S I] are associated with shocks (e.g., Haas et al., 1991;
Rosenthal et al., 2000; Neufeld et al., 2009; Simpson, 2018) and shocks
provide a mechanism to both release neutral sulphur material from the dust
grains into the gas phase and to excite [S I] emission.

We evaluate whether this is possible by performing a back-of-the-envelope
calculation. First, we estimate the timescale it would take to ionize sulphur
in the hard radiation field of AB7. We integrate the ionizing spectrum over
the ionization potential of sulphur (10.36 eV) to find the rate of sulphur-
ionizing photons at a given distance of the [S I] emitting region. We
calculate a timescale of 1.7 years to ionize sulphur after assuming a cross-
section of $\rm 3\times 10^{-18}\ cm^{-3}$ (Draine, 2011) and a radius of 10
pc. We next calculate the mass of neutral sulphur that must be replenished
within this time. From the observed [S I] intensity we find the column density
by assuming collisional excitation and use the equation given in Crawford et
al. (1985), a critical density for collisions with $\rm e^{-}$ of $\rm
n=1.55\times 10^{5}\ cm^{-3}$ (Draine, 2011), a spontaneous decay rate of $\rm
A=1.40\times 10^{-3}\ s^{-1}$ (Froese Fischer et al., 2006), a temperature of
$\rm T=1.4\times 10^{4}\ K$, an abundance of $3.89\times 10^{-6}$ (Russell &
Dopita, 1992), and a density of $\rm n=40\ cm^{-3}$. The density is chosen by
assuming a strong shock wave will have four times the density of the pre-
shocked medium of 10 $\mathrm{cm^{-3}}$(e.g., Draine, 2011). We calculate a
column density of $\rm N_{S^{0}}=4.7\times 10^{14}\ cm^{-2}$ and consequently
a mass of $\rm M_{S^{0}}=3.7\times 10^{-2}\ \mbox{\,M${}_{\odot}$}$ over a 10
pc radius. Thus our estimate suggests that $\rm 3.7\times 10^{-2}\
\mbox{\,M${}_{\odot}$}$ of sulphur will need to be lifted off of destroyed
dust grains every 1.7 years. Given the low abundance of sulphur relative to
other elements that make the bulk of the dust mass (C, Si, O), this represents
a very large dust mass destroyed every year, and it seems impossible this
reservoir would exist over the life of the nebula (or be in some way
replenished by the central source mass loss rate). Further, there is no
signature of other neutral species besides [S I] near the center of the nebula
to support the scenario that dust destruction from shocks could produce the
neutral sulphur.

Further, we would also expect that the other neutral species present in dust
are also at the central part of nebula, which is not the case.

We note that if the [S I] emission is produced through dust destruction due to
shocks from stellar winds, there would likely be a substantial amount of [Fe
II] emission produced from the same mechanism since iron is heavily depleted
in dust. The 25.9 µm blend of [O IV] 25.89 µm emission and [Fe II] 25.99µm
described in §2.7 may therefore have a substantial contribution from [Fe II]
25.99 µm. However, since the [Fe II] 17.93 µm line is not detected and we
predict $<5\%$ of the 25.9 µm line flux is from [Fe II], we conclude that dust
destruction scenario is unlikely.

Artifact in Spitzer IRS Data: The feature at 25.25 µm has multiple detections
through the spectral slit scans and is not a one time, random anomaly with the
data. We instead explore two detector effects associated with Spitzer IRS data
that can manifest as spurious spectral features and explain how they are
unlikely to cause the possible [S I] feature at 25.25 µm. The IRS Handbook
reports that spectral “ghost” features have appeared near extremely bright
spectral lines in the short-high (SH) module near the 12.8 µm [Ne II] line
(Simpson et al., 2007; IRS Instrument Team & Science User Support Team, 2012).
However, it is extremely unlikely that the 25.25 µm feature is a ghost of the
25.9 µm [O IV] line as the IRS Handbook specifically reports that no similar
features had been found for the long-low (LL) module that covers the 25.25 µm
line. The ghost features manifest in the SH module due to order mixing in the
cross-disperser and there are no equivalent order-sorting filters for the LL
module. Further, the ghost features that appear in the SH module are symmetric
across the [Ne II] line due to smearing of the flux on the plane of the slit
and are also 0.5-1% of the total flux of the line. The [S I] 25.25 µm feature
observed in N76 is completely asymmetric, does not appear as a smeared feature
on the slit plane, and is 5% of the total [O IV] flux, therefore strengthening
the position that it is unlikely for the [S I] line to be a ghost feature of
[O IV]. The other possible artifact we explore is spectral fringing, reported
in the IRS handbook to be associated with the LL module where the 25.25 µm
line is found. Spectral fringing manifests as residual non-sinusoidal wiggles
in the IRS spectrum and begins at $\sim$20 µm (IRS Instrument Team & Science
User Support Team, 2012). This fringing was more prevalent earlier in the
Spitzer mission and has been seen in objects such as M82 and Cassiopeia A
(Beirão et al., 2008; Smith et al., 2009). In these objects, the fringes are
in a continuous pattern from 20 - 30 µm, are 1-2% of the flux, and are
coincident with the continuum and correlated with its brightness. The 25.25 [S
I] µm feature in N76 is therefore unlike these spectral fringes because only a
single feature appears in the spectrum, not the continuous pattern seen in all
other cases. It is also a factor of 2 or more brighter than the fringe
features and the [S I] line is not spatially correlated with the continuum.
Therefore, we conclude that the 25.25 µm line is likely not due to a Spitzer
IRS artifact, but future observations (such as with JWST) will be need to
confirm and characterize this curious line.

Unidentified spectral line: We used the NIST atomic spectral database to
search for lines near the 20.0-20.5 µm range and did not identify any possible
candidates besides [S I] (see discussion in $\S$4.8). It is possible that the
25.25 µm emission originates from a molecular compound, but a molecule would
have similar difficulties surviving the intense radiation field of AB7.
Instead of [S I], however, this line could be an unidentified spectral line
that corresponds to a very high ionization state not in the NIST catalog.

###  5.3 Physical conditions in N76

Figure 8 presents how the electron temperature and electron density vary for
each constant hydrogen density model as a function of physical distance from
AB7. The electron density is slightly higher at the center of N76, due to the
contribution from the $\rm He^{++}$ region that forms around very hot WR stars
(see also sixth panel in Figure 6). In our constant hydrogen density modeling
$\rm n_{e}$ stays constant up to the ionization front, where the gas
transitions from ionized to neutral hydrogen. Overall, the electron density
ranges from $n_{e}\sim 4-12\ \rm cm^{-3}$ The temperature profile is smooth
distribution, with a central peak near the edge of the photoionized region
reaching 22,000 to 24,000 K, which then stays approximately constant at around
14,800 K until reaching the hydrogen ionization front. Note that the inner
regions close to the star are simply a cavity of radius 10 pc in the
photoionization models, so the precise details of the transition between
shock-ionized and photoionized gas are not accurately represented in these
models.

The ratio of [S III] 18 µm/ [S III] 33 µm also provides an estimate for
$n_{e}$ of an HII region. However, this ratio is only sensitive to gas with
densities $n_{e}>100\ \mbox{$\mathrm{cm^{-3}}$}$ (Dudik et al., 2007), thus it
does not probe the less dense ionized gas in N76. Our photoionization modeling
of N76, however, is able to explore a wider region of parameter space than
line ratios alone. Figure 8 shows how dramatically the size of the HII region
changes for small changes in hydrogen density, illustrating the usefulness of
spatially resolved data for characterizing HII regions.

We also calculate the ionization parameter, $U$, defined as the dimensionless
ratio of the hydrogen-ionizing photons to the total hydrogen density. This
parameter is often used to characterize HII regions, and is defined as

$U=\frac{Q(H)}{4\pi r_{o}^{2}n(H)c}\ ,$ (4)

where $r_{o}$ is the ionized region radius in units of cm, $n(H)$ is the total
hydrogen number density in units of $\mathrm{cm^{-3}}$, $c$ is the speed of
light, and $Q(H)$ is the number of hydrogen-ionizing photons emitted by the
central source per second. The ionization parameter we measure for N76 ranges
from $U=1.7\times 10^{-2}$ to $U=7.1\times 10^{-2}$, which is quite high and
comparable to the massive HII region 30 Doradus (see Figure 9). The high
radiation fields of AB7 and relatively low density of N76 indicate that N76
has a high ionization parameter, even with a single ionization source.

Figure 8: The electron temperature and density for the constant density Cloudy
models as a function of physical distance away from AB7. Colors follow the
same scheme as in Figure 3.

We also note that the Cloudy models presented in this work are radiation-
bounded, which are defined as HII regions wherein Cloudy models end at the
ionization front and do not allow appreciable ionizing photons to escape. The
Cloudy models presented here are an idealized case and represent the average
of N76, but it is likely that some parts of the HII region are matter-bounded
while others are radiation-bounded. The south-western edge of N76, for
example, is much more diffuse than the rest of the nebula, which can allow
ionizing photons to escape. A similar morphology is found in N76’s HII region
neighbor, the supergiant HII region N66 (Geist et al., 2022).

###  5.4 N76 in context

The N76 nebula is unique for its very luminous central ionizing source, low
density ($\sim$10 $\mathrm{cm^{-3}}$), and consequent large ($\sim$40 pc)
size. In order to compare it to other HII regions, we examine our results
through two different lenses: observations of the ionized gas in low
metallicity dwarf galaxies, and in-depth studies of WR nebulae.

####  5.4.1 Comparison with other low metallicity systems

The largest study focused on understanding the properties of the mid-infrared
transitions from gas in low metallicity systems is the Dwarf Galaxy Survey
(DGS), which modeled IR observations of 48 nearby galaxies using Cloudy
(Cormier et al., 2015, 2019). The results of DGS made considerable progress in
the understanding of ionized gas properties in low metallicity systems, but
the majority of the sample consists of modeling an unresolved spectrum of each
galaxy using a single or two component Cloudy model. While recent, follow up
studies use a topological model to infer the statistical properties of the
multiphase ISM from the unresolved data (Lebouteiller & Ramambason, 2022;
Ramambason et al., 2022, 2023), this work instead focuses on the nearby HII
region N76 as it is fully resolved and contains a single, well characterized
ionization source. The models of N76 presented here provide an ideal case
study to understand the impact Wolf-Rayet stars have in low metallicity
environments, which are often characterized by active star formation and can
contain large number of WR stars (e.g., Izotov et al., 1997).

We compare the ionized gas properties found in the large surveys of dwarf
galaxies below. Cormier et al. (2019) reports electron densities that range
$\rm n_{e}=10^{0.5}-10^{3.0}\ \mbox{$\mathrm{cm^{-3}}$}$ and ionization
parameters of $\log U=-3.0\ \rm{to}\ -0.3$ with a correlation such that the
highest $U$ values are found in the lowest metallicity galaxies. This trend is
also seen in the blue compact dwarf (BCD) galaxies analyzed in Hunt et al.
(2010), where they report elevated [Ne III]/[Ne II], [S IV]/[S III], and [O
IV]/[Si II] ratios when compared to higher metallicity counterparts,
indicative of harder radiation fields at low metallicities. Additionally, Hunt
et al. (2010) report electron densities ranging from $\rm n_{e}=30-600\
\mbox{$\mathrm{cm^{-3}}$}$ through the [S III] 18 µm/ [S III] 33 µm ratio.
These two programs focused on analyzing spatially unresolved infrared
spectroscopy of galaxy-wide measurements in order to form conclusions about
the global ionized gas properties in dwarf galaxies. With a metallicity of 12
+ log(O/H) $\sim$ 8.0, N76, as part of the SMC, falls at the median of the
metallicity ranges investigated in these studies. The resolved measurements of
N76 are well in the range of the globally integrated properties determined by
these surveys, although the electron density in N76 is at the lower boundary
found in these studies. This suggests that the brightest HII regions can
dominate the excitation properties of the ISM in the ionized gas of these
galaxies.

We also compare N76 to photoionization modeling studies of nearby dwarf
galaxies that are more spatially resolved than those just described. To be
consistent, all of these studies are done using similar IR spectroscopic
observations and modeled with the Cloudy photoionization code. Indebetouw et
al. (2009) model 30 Doradus in the Large Magellanic Cloud (LMC), the largest
supergiant HII region in the nearby universe. They find a relatively high
electron density, $n_{e}=10^{2.4}-10^{2.7}\ \mbox{$\mathrm{cm^{-3}}$}$ and a
hard radiation field with $\log U=-1.7\ \rm{to}-1.3$. Similarly, Polles et al.
(2019) investigates the ionized gas in five individual star-forming clumps in
IC10, and their Cloudy modeling indicates $n_{e}=10^{2.0}-10^{2.6}\
\mbox{$\mathrm{cm^{-3}}$}$ and $\log U=-3.8\ \rm{to}-1.0$. Lastly, Dimaratos
et al. (2015) focus their study on NGC 4214 and identify a separate central
and southern star forming region, calculating lower densities of $n_{e}=440$
and $n_{e}=170$, respectively, and ionization parameters of $U=-2.3$ and
$U=-2.7$. Taken as a whole, N76 has a lower electron density than found in
those objects, with somewhat larger ionization parameters.

In order to summarize the infrared Cloudy modeling results across these
different surveys and studies, we plot the [Ne III] 15 µm/ [Ne II] 12 µm and
the Cloudy model ionization parameter, $U$, in Figure 9. This ratio traces the
intensity of the radiation field and, in principle for a single main sequence
star hardness should correlate well with $U$. We also plot a summary of the
[Ne III] 15 µm/ [Ne II] 12 µm ratio from a sample of high metallicity galaxies
(Dale et al., 2009; Inami et al., 2013). While these analyses did not produce
an ionization parameter, we do show their range in [Ne III] 15 µm/ [Ne II] 12
µm line ratio. We confirm that the higher metallicity galaxies have a much
lower [Ne III] 15 µm/ [Ne II] 12 µm line ratio than the low metallicity
systems. N76 and 30 Doradus are at the upper third of the observed line ratio,
with [Ne III] 15 µm/ [Ne II] 12 µm $\sim$ 10-15, indicating hard radiation
fields in these HII regions. Many of the galaxies in the DGS sample have a
similar ionization parameter and line ratio as these highly energetic HII
regions, suggesting that conditions in these galaxies may be very similar to
the resolved studies of N76 and 30Dor.

When examining scales of $\sim$200 pc, Cormier et al. (2019) find evidence of
additional, low ionization components of ionized gas from their Cloudy models.
Similarly, Polles et al. (2019) and Dimaratos et al. (2015) find excess
emission of [Ne II] that they attribute to a diffuse ionized gas component. In
N76, our Cloudy models also underpredict [Ne II], providing more evidence for
an additional low ionization source that fuels the diffuse ionized gas in
these galaxies.

Figure 9: Comparing the [Ne III] 15 µm/ [Ne II] 12 µm ratio to the ionization
parameter of a variety of low metallicity systems, including the Dwarf Galaxy
Survey (Cormier et al., 2015, 2019), IC10 (Polles et al., 2019), NGC 4214
(Dimaratos et al., 2015), 30 Doradus (Indebetouw et al., 2009), and N76 (this
work). These studies include a mixture of integrated data from entire galaxies
(DGS), to individual measurements in divided regions of these galaxies (IC10
and NGC 4214), to single HII regions (N76 and 30 Dor). Overall, the [Ne III]
15 µm/ [Ne II] 12 µm ratio shows that the low metallicity systems plotted have
much harder radiation fields than the averages shown for the higher
metallicity galaxies from Dale et al. (2009) and Inami et al. (2013).

####  5.4.2 Comparison with other Wolf-Rayet nebulae

Most well-studied WR nebula are defined as bubble nebulae, characterized by
thin ionized shell covering a large cavity, excavated by the strong stellar
winds of the central WR star. Optical IFU measurements of the Milky Way WR
nebula NGC 6888 show a three-shell structure centered around a 400 pc
evacuated cavity: an elliptical inner broken shell formed by shocks from
interactions between the WR and supergiant winds, an outer spherical shell
that represents the wind-blown bubble, and a skin of ISM material surrounding
the nebula (Fernández-Martín et al., 2012). Through photoionization modeling,
Fernández-Martín et al. (2012) calculate an electron density of $n_{e}=400\
\mbox{$\mathrm{cm^{-3}}$}$ for the outer wind-blown bubble shell. Infrared
spectroscopy described in Rubio et al. (2020) characterize NGC 6888 further,
corroborating the density found for the outer shell in Fernández-Martín et al.
(2012) and identifying the density of the ISM skin as $n_{e}=180\
\mbox{$\mathrm{cm^{-3}}$}$.

The overall morphological picture of NGC 6888 is not dissimilar to that of
N76: a wind-shocked inner cavity caused by the powerful stellar winds of the
WR star that transitions into a photoionized gas region and the surrounding
ISM. However, N76 has a much smaller cavity (10 pc in our models compared to
400 pc), a lower density photoionized region with $n_{e}=4-10\
\mbox{$\mathrm{cm^{-3}}$}$, and this photoionized shell is rather thick
compared to the cavity size (35 pc). These differences may be related to the
age of N76. Shenar et al. (2016) estimates that AB7 is 3.4 Myr old, 1.3 Myr
younger than WR 136 that powers NGC 6888 (Moore et al., 2000). The size of a
stellar wind-blown bubble directly correlates with its age (Weaver et al.,
1977; Freyer et al., 2003, 2006), suggesting that the younger N76 nebula will
be smaller than NGC 6888. Observations of WR nebulae in the LMC show the
traditional bubble structure, i.e. a large cavity with a thin, photoionized
shell of material in the majority of nebulae (Hung et al., 2021).

Since the regions around WR nebulae contain material blown off by the powerful
stellar winds of the star, their abundances can be enhanced, yielding elevated
nitrogen abundances in those with a central WN star (e.g., Fernández-Martín et
al., 2012; Stock & Barlow, 2014). AB7 contains a WN4 type star, but we do not
have a nitrogen spectral line to probe the nitrogen abundance. The abundances
described in Table 2 do a satisfactory job of predicting the emission line
surface brightness so we do not have evidence of elemental enhancement in N76.
Optical spectroscopy of N76 by Nazé et al. (2003), however, show a 40%
elevated N/O ratio compared to global SMC values, suggesting that the chemical
make up of N76 is influenced by AB7.

##  6 Summary and Conclusions

We use Spitzer IRS spectroscopic maps of mid-infrared transitions and Cloudy
to model the photoionization in the HII region N76 in the Small Magellanic
Cloud. The ionization source is AB7, a binary composed of a WN4 nitrogen-rich
Wolf-Rayet and an O6 I supergiant star, and we use PoWR models (Todt et al.,
2015) to estimate its SED and luminosity (Shenar et al., 2016). AB7 contains
one of the hottest WR stars in the SMC. We assume the N76 nebula is
spherically shaped and project the predicted Cloudy line surface brightness as
a function of distance to the source on the plane of the sky. We only consider
emission from the ionized gas. The results from the photoionization models
match the spatially resolved images of most ionized emission lines measured
from Spitzer IRS and Herschel PACS observations, with just density as the free
parameter. Our conclusions are as follows:

  1. 1.

Comparing the models to the observed spatial distribution of the measured
emission line brightness allows us to determine much more precise densities
than just using the total line intensity (Figure 3 and Table 1).

  2. 2.

Constant density Cloudy photoionization models reproduce most of the ionized
gas emission line intensity distributions ([S III], [Ne III], [S IV], [Ne V],
[O IV], and [O III]) for a narrow hydrogen density between $\rm n_{H}=4\
cm^{-3}-10\ cm^{-3}$ (Figure 1).

  3. 3.

The Cloudy photoionization models require a 10 pc cavity of hot, shocked
plasma originating from a stellar wind-blown nebula in order to predict the
spatial distribution of the ionized gas emission lines (Weaver et al., 1977;
Freyer et al., 2003, 2006).

  4. 4.

The high ionization lines, [O IV] and [Ne V] ($\rm E_{ion}>50\ eV$), are well
predicted through the Cloudy models using the SED modeled in Shenar et al.
(2016). The flux of these lines is highly sensitive to changes in the high
energy tail of the WR SED.

  5. 5.

The neutral-dominated emission lines (lines coming from ionic species with an
ionization potential less than hydrogen, 13.6 eV), can originate from both the
ionized gas and the neutral gas. Our modeling finds that very little of the
observed [C II] emission comes from the ionized gas ($<5\%$), while most (91%)
of the [Si II] emission can be produced in the transition from HII region to
PDR.

  6. 6.

The [Ne II] line is under-predicted in our Cloudy models, suggesting the
excess emission may arise from the diffuse ionized gas not associated with N76
(e.g., Dimaratos et al., 2015; Cormier et al., 2019; Polles et al., 2019).

  7. 7.

We observe a feature at $25.25$ µm that we tentatively identify as [S I] in
the high ionization zone of the N76 nebula (see Figure 1). We explore possible
causes for [S I] and alternative possibilities for the line in §5.2, including
the possibility of producing neutral sulphur from dust destruction. None of
these satisfactorily explains the line.

  8. 8.

The physical conditions of N76 derived from the Cloudy models characterize a
hot HII region with a maximum electron temperature range from $T_{e}\sim
22,000-24,000\rm K$, a relatively small density range of $n_{e}\sim 4-12\ \rm
cm^{-3}$, and high ionization parameters of $log(U)=-1.15--1.77$ (§5.3).

  9. 9.

Low metallicity galaxies have harder radiation fields than their higher
metallicity counterparts. The N76 region has a similar ionization parameter to
samples of low metallicity galaxies, suggesting that the bright HII regions
dominate the spectrum of individual galaxies (see Figure 9).

Our work shows that just photoionization from low metallicity WR stars can
produce appreciable [O IV] and [Ne V] emission in the surrounding ISM. Thus
N76 is a vital case study in determining the role of WRs in producing large
fluxes of high ionization ions ($\rm E_{ion}>50\ ev$) found in low metallicity
Extreme Emission Line Galaxies (EELG) and blue compact dwarf galaxies (e.g.,
Thuan & Izotov, 2005; Izotov et al., 2012, 2021; Kehrig et al., 2015, 2018;
Leitherer et al., 2018; Olivier et al., 2022). Due to the rarity of WRs in
population synthesis modeling (e.g., Leitherer et al., 2014), WR stars are
often ruled out as a contributor to the high ionization lines. Our work
suggests that the choice of WR model can make a significant impact on the high
energy tail of the SED and consequently the amount of [Ne V] and [O IV]
produced, suggesting that special care should be taken when modeling WR stars
in unresolved systems.

We would like to thank T. Shenar for providing the SED of AB7 from Shenar et
al. (2016). We would also like to thank the anonymous referee for providing
feedback that greatly improved the manuscript. ET and ADB acknowledge support
from NSF-AST2108140 and NASA ADAP 80NSSC19K1015. RI acknowledges support from
NSF-AST2009624. M.R. wishes to acknowledge support from ANID(CHILE) through
FONDECYT grant No1190684 and partial support from ANIDA Basal FB210003.

## References

  * Aniano et al. (2011) Aniano, G., Draine, B. T., Gordon, K. D., & Sandstrom, K. 2011, PASP, 123, 1218, doi: 10.1086/662219
  * Atek et al. (2011) Atek, H., Siana, B., Scarlata, C., et al. 2011, ApJ, 743, 121, doi: 10.1088/0004-637X/743/2/121
  * Baldovin-Saavedra et al. (2011) Baldovin-Saavedra, C., Audard, M., Güdel, M., et al. 2011, A&A, 528, A22, doi: 10.1051/0004-6361/201015622
  * Baldwin et al. (1981) Baldwin, J. A., Phillips, M. M., & Terlevich, R. 1981, PASP, 93, 5, doi: 10.1086/130766
  * Beirão et al. (2008) Beirão, P., Brandl, B. R., Appleton, P. N., et al. 2008, ApJ, 676, 304, doi: 10.1086/527343
  * Berg et al. (2019) Berg, D. A., Chisholm, J., Erb, D. K., et al. 2019, ApJ, 878, L3, doi: 10.3847/2041-8213/ab21dc
  * Blair et al. (2000) Blair, W. P., Morse, J. A., Raymond, J. C., et al. 2000, ApJ, 537, 667, doi: 10.1086/309077
  * Bolatto et al. (2007) Bolatto, A. D., Simon, J. D., Stanimirović, S., et al. 2007, ApJ, 655, 212, doi: 10.1086/509104
  * Chevalier (2005) Chevalier, R. A. 2005, ApJ, 619, 839, doi: 10.1086/426584
  * Chu (1981) Chu, Y. H. 1981, ApJ, 249, 195, doi: 10.1086/159275
  * Chu et al. (1983) Chu, Y. H., Treffers, R. R., & Kwitter, K. B. 1983, ApJS, 53, 937, doi: 10.1086/190914
  * Conti & Massey (1989) Conti, P. S., & Massey, P. 1989, ApJ, 337, 251, doi: 10.1086/167101
  * Cormier et al. (2015) Cormier, D., Madden, S. C., Lebouteiller, V., et al. 2015, A&A, 578, A53, doi: 10.1051/0004-6361/201425207
  * Cormier et al. (2019) Cormier, D., Abel, N. P., Hony, S., et al. 2019, A&A, 626, A23, doi: 10.1051/0004-6361/201834457
  * Crawford et al. (1985) Crawford, M. K., Genzel, R., Townes, C. H., & Watson, D. M. 1985, ApJ, 291, 755, doi: 10.1086/163113
  * Crowther et al. (1999) Crowther, P. A., Beck, S. C., Willis, A. J., et al. 1999, MNRAS, 304, 654, doi: 10.1046/j.1365-8711.1999.02380.x
  * Crowther & Hadfield (2006) Crowther, P. A., & Hadfield, L. J. 2006, A&A, 449, 711, doi: 10.1051/0004-6361:20054298
  * Dale et al. (2009) Dale, D. A., Cohen, S. A., Johnson, L. C., et al. 2009, ApJ, 703, 517, doi: 10.1088/0004-637X/703/1/517
  * Dennefeld & Stasińska (1983) Dennefeld, M., & Stasińska, G. 1983, A&A, 118, 234
  * Dimaratos et al. (2015) Dimaratos, A., Cormier, D., Bigiel, F., & Madden, S. C. 2015, A&A, 580, A135, doi: 10.1051/0004-6361/201526447
  * Draine (2011) Draine, B. T. 2011, Physics of the Interstellar and Intergalactic Medium (Princeton University Press)
  * Dudik et al. (2007) Dudik, R. P., Weingartner, J. C., Satyapal, S., et al. 2007, ApJ, 664, 71, doi: 10.1086/518685
  * Dufour (1984) Dufour, R. J. 1984, in Structure and Evolution of the Magellanic Clouds, ed. S. van den Bergh & K. S. D. de Boer, Vol. 108, 353–361, doi: 10.1017/S0074180900040481
  * Ferland et al. (2017) Ferland, G. J., Chatzikos, M., Guzmán, F., et al. 2017, Rev. Mexicana Astron. Astrofis., 53, 385\.  https://arxiv.org/abs/1705.10877
  * Fernández-Martín et al. (2012) Fernández-Martín, A., Martín-Gordón, D., Vílchez, J. M., et al. 2012, A&A, 541, A119, doi: 10.1051/0004-6361/201118421
  * Foellmi et al. (2003) Foellmi, C., Moffat, A. F. J., & Guerrero, M. A. 2003, MNRAS, 338, 360, doi: 10.1046/j.1365-8711.2003.06052.x
  * Freyer et al. (2003) Freyer, T., Hensler, G., & Yorke, H. W. 2003, ApJ, 594, 888, doi: 10.1086/376937
  * Freyer et al. (2006) —. 2006, ApJ, 638, 262, doi: 10.1086/498734
  * Froese Fischer et al. (2006) Froese Fischer, C., Tachiev, G., & Irimia, A. 2006, Atomic Data and Nuclear Data Tables, 92, 607, doi: 10.1016/j.adt.2006.03.001
  * Galliano et al. (2021) Galliano, F., Nersesian, A., Bianchi, S., et al. 2021, A&A, 649, A18, doi: 10.1051/0004-6361/202039701
  * Garnett et al. (1991) Garnett, D. R., Kennicutt, Robert C., J., Chu, Y.-H., & Skillman, E. D. 1991, ApJ, 373, 458, doi: 10.1086/170065
  * Geist et al. (2022) Geist, E., Gallagher, J. S., Kotulla, R., et al. 2022, PASP, 134, 064301, doi: 10.1088/1538-3873/ac697b
  * Gräfener et al. (2002) Gräfener, G., Koesterke, L., & Hamann, W. R. 2002, A&A, 387, 244, doi: 10.1051/0004-6361:20020269
  * Guerrero & Chu (2008) Guerrero, M. A., & Chu, Y.-H. 2008, ApJS, 177, 216, doi: 10.1086/587059
  * Haas et al. (1991) Haas, M. R., Hollenbach, D., & Erickson, E. F. 1991, ApJ, 374, 555, doi: 10.1086/170144
  * Haffner et al. (2009) Haffner, L. M., Dettmar, R. J., Beckman, J. E., et al. 2009, Reviews of Modern Physics, 81, 969, doi: 10.1103/RevModPhys.81.969
  * Hainich et al. (2019) Hainich, R., Ramachandran, V., Shenar, T., et al. 2019, A&A, 621, A85, doi: 10.1051/0004-6361/201833787
  * Hainich et al. (2014) Hainich, R., Rühling, U., Todt, H., et al. 2014, A&A, 565, A27, doi: 10.1051/0004-6361/201322696
  * Hamann & Gräfener (2003) Hamann, W. R., & Gräfener, G. 2003, A&A, 410, 993, doi: 10.1051/0004-6361:20031308
  * Hamann & Gräfener (2004) —. 2004, A&A, 427, 697, doi: 10.1051/0004-6361:20040506
  * Hayashi et al. (1994) Hayashi, I., Koyama, K., Ozaki, M., et al. 1994, PASJ, 46, L121
  * Henize (1956) Henize, K. G. 1956, ApJS, 2, 315, doi: 10.1086/190025
  * Heydari-Malayeri et al. (2001) Heydari-Malayeri, M., Charmandaris, V., Deharveng, L., et al. 2001, A&A, 372, 495, doi: 10.1051/0004-6361:20010269
  * Hollenbach & McKee (1989) Hollenbach, D., & McKee, C. F. 1989, ApJ, 342, 306, doi: 10.1086/167595
  * Hung et al. (2021) Hung, C. S., Ou, P.-S., Chu, Y.-H., Gruendl, R. A., & Li, C.-J. 2021, ApJS, 252, 21, doi: 10.3847/1538-4365/abcc00
  * Hunt et al. (2010) Hunt, L. K., Thuan, T. X., Izotov, Y. I., & Sauvage, M. 2010, ApJ, 712, 164, doi: 10.1088/0004-637X/712/1/164
  * Hurley et al. (2000) Hurley, J. R., Pols, O. R., & Tout, C. A. 2000, MNRAS, 315, 543, doi: 10.1046/j.1365-8711.2000.03426.x
  * Inami et al. (2013) Inami, H., Armus, L., Charmandaris, V., et al. 2013, ApJ, 777, 156, doi: 10.1088/0004-637X/777/2/156
  * Indebetouw et al. (2009) Indebetouw, R., de Messières, G. E., Madden, S., et al. 2009, ApJ, 694, 84, doi: 10.1088/0004-637X/694/1/84
  * IRS Instrument Team & Science User Support Team (2012) IRS Instrument Team & Science User Support Team. 2012, IRS Instrument Handbook, IPAC, doi: 10.26131/IRSA487
  * Izotov et al. (1997) Izotov, Y. I., Foltz, C. B., Green, R. F., Guseva, N. G., & Thuan, T. X. 1997, ApJ, 487, L37, doi: 10.1086/310872
  * Izotov et al. (2021) Izotov, Y. I., Thuan, T. X., & Guseva, N. G. 2021, MNRAS, 508, 2556, doi: 10.1093/mnras/stab2798
  * Izotov et al. (2012) Izotov, Y. I., Thuan, T. X., & Privon, G. 2012, MNRAS, 427, 1229, doi: 10.1111/j.1365-2966.2012.22051.x
  * Jenkins & Wallerstein (2017) Jenkins, E. B., & Wallerstein, G. 2017, ApJ, 838, 85, doi: 10.3847/1538-4357/aa64d4
  * Kaufman et al. (2006) Kaufman, M. J., Wolfire, M. G., & Hollenbach, D. J. 2006, ApJ, 644, 283, doi: 10.1086/503596
  * Kehrig et al. (2018) Kehrig, C., Vílchez, J. M., Guerrero, M. A., et al. 2018, MNRAS, 480, 1081, doi: 10.1093/mnras/sty1920
  * Kehrig et al. (2015) Kehrig, C., Vílchez, J. M., Pérez-Montero, E., et al. 2015, ApJ, 801, L28, doi: 10.1088/2041-8205/801/2/L28
  * Lancaster et al. (2021) Lancaster, L., Ostriker, E. C., Kim, J.-G., & Kim, C.-G. 2021, ApJ, 914, 90, doi: 10.3847/1538-4357/abf8ac
  * Lau et al. (2022) Lau, R. M., Hankins, M. J., Han, Y., et al. 2022, Nature Astronomy, 6, 1308, doi: 10.1038/s41550-022-01812-x
  * Lebouteiller et al. (2008) Lebouteiller, V., Bernard-Salas, J., Brandl, B., et al. 2008, ApJ, 680, 398, doi: 10.1086/587503
  * Lebouteiller & Ramambason (2022) Lebouteiller, V., & Ramambason, L. 2022, A&A, 667, A34, doi: 10.1051/0004-6361/202243865
  * Lee et al. (2009) Lee, M.-Y., Stanimirović, S., Ott, J., et al. 2009, AJ, 138, 1101, doi: 10.1088/0004-6256/138/4/1101
  * Leitherer et al. (2018) Leitherer, C., Byler, N., Lee, J. C., & Levesque, E. M. 2018, ApJ, 865, 55, doi: 10.3847/1538-4357/aada84
  * Leitherer et al. (2014) Leitherer, C., Ekström, S., Meynet, G., et al. 2014, ApJS, 212, 14, doi: 10.1088/0067-0049/212/1/14
  * Leroy et al. (2007) Leroy, A., Bolatto, A., Stanimirovic, S., et al. 2007, ApJ, 658, 1027, doi: 10.1086/511150
  * Maseda et al. (2013) Maseda, M. V., van der Wel, A., da Cunha, E., et al. 2013, ApJ, 778, L22, doi: 10.1088/2041-8205/778/1/L22
  * Moore et al. (2000) Moore, B. D., Hester, J. J., & Scowen, P. A. 2000, AJ, 119, 2991, doi: 10.1086/301389
  * Nazé et al. (2003) Nazé, Y., Rauw, G., Manfroid, J., Chu, Y. H., & Vreux, J. M. 2003, A&A, 408, 171, doi: 10.1051/0004-6361:20030847
  * Neufeld et al. (2009) Neufeld, D. A., Nisini, B., Giannini, T., et al. 2009, ApJ, 706, 170, doi: 10.1088/0004-637X/706/1/170
  * Neugent et al. (2018) Neugent, K. F., Massey, P., & Morrell, N. 2018, ApJ, 863, 181, doi: 10.3847/1538-4357/aad17d
  * Niemela et al. (2002) Niemela, V. S., Massey, P., Testor, G., & Giménez Benítez, S. 2002, MNRAS, 333, 347, doi: 10.1046/j.1365-8711.2002.05395.x
  * Olivier et al. (2022) Olivier, G. M., Berg, D. A., Chisholm, J., et al. 2022, ApJ, 938, 16, doi: 10.3847/1538-4357/ac8f2c
  * Parkin & Pittard (2008) Parkin, E. R., & Pittard, J. M. 2008, MNRAS, 388, 1047, doi: 10.1111/j.1365-2966.2008.13511.x
  * Pereira-Santaella et al. (2010) Pereira-Santaella, M., Diamond-Stanic, A. M., Alonso-Herrero, A., & Rieke, G. H. 2010, ApJ, 725, 2270, doi: 10.1088/0004-637X/725/2/2270
  * Poglitsch et al. (1995) Poglitsch, A., Krabbe, A., Madden, S. C., et al. 1995, ApJ, 454, 293, doi: 10.1086/176482
  * Polles et al. (2019) Polles, F. L., Madden, S. C., Lebouteiller, V., et al. 2019, A&A, 622, A119, doi: 10.1051/0004-6361/201833776
  * Ramambason et al. (2022) Ramambason, L., Lebouteiller, V., Bik, A., et al. 2022, arXiv e-prints, arXiv:2207.06146.  https://arxiv.org/abs/2207.06146
  * Ramambason et al. (2023) Ramambason, L., Lebouteiller, V., Madden, S. C., et al. 2023, arXiv e-prints, arXiv:2306.14881, doi: 10.48550/arXiv.2306.14881
  * Rigby et al. (2015) Rigby, J. R., Bayliss, M. B., Gladders, M. D., et al. 2015, ApJ, 814, L6, doi: 10.1088/2041-8205/814/1/L6
  * Rodríguez et al. (2020) Rodríguez, L. F., Arthur, J., Montes, G., Carrasco-González, C., & Toalá, J. A. 2020, ApJ, 900, L3, doi: 10.3847/2041-8213/abad9d
  * Roman-Duval et al. (2014) Roman-Duval, J., Gordon, K. D., Meixner, M., et al. 2014, ApJ, 797, 86, doi: 10.1088/0004-637X/797/2/86
  * Rosenthal et al. (2000) Rosenthal, D., Bertoldi, F., & Drapatz, S. 2000, A&A, 356, 705.  https://arxiv.org/abs/astro-ph/0002456
  * Rubio et al. (2020) Rubio, G., Toalá, J. A., Jiménez-Hernández, P., et al. 2020, MNRAS, 499, 415, doi: 10.1093/mnras/staa2837
  * Russell & Dopita (1990) Russell, S. C., & Dopita, M. A. 1990, ApJS, 74, 93, doi: 10.1086/191494
  * Russell & Dopita (1992) —. 1992, ApJ, 384, 508, doi: 10.1086/170893
  * Sander et al. (2015) Sander, A., Shenar, T., Hainich, R., et al. 2015, A&A, 577, A13, doi: 10.1051/0004-6361/201425356
  * Sandstrom et al. (2010) Sandstrom, K. M., Bolatto, A. D., Draine, B. T., Bot, C., & Stanimirović, S. 2010, ApJ, 715, 701, doi: 10.1088/0004-637X/715/2/701
  * Sandstrom et al. (2009) Sandstrom, K. M., Bolatto, A. D., Stanimirović, S., van Loon, J. T., & Smith, J. D. T. 2009, ApJ, 696, 2138, doi: 10.1088/0004-637X/696/2/2138
  * Sandstrom et al. (2012) Sandstrom, K. M., Bolatto, A. D., Bot, C., et al. 2012, ApJ, 744, 20, doi: 10.1088/0004-637X/744/1/20
  * Schaerer & Stasińska (1999) Schaerer, D., & Stasińska, G. 1999, A&A, 345, L17.  https://arxiv.org/abs/astro-ph/9903430
  * Schaerer & Vacca (1998) Schaerer, D., & Vacca, W. D. 1998, ApJ, 497, 618, doi: 10.1086/305487
  * Senchyna et al. (2017) Senchyna, P., Stark, D. P., Vidal-García, A., et al. 2017, MNRAS, 472, 2608, doi: 10.1093/mnras/stx2059
  * Shenar et al. (2016) Shenar, T., Hainich, R., Todt, H., et al. 2016, A&A, 591, A22, doi: 10.1051/0004-6361/201527916
  * Simpson (2018) Simpson, J. P. 2018, ApJ, 857, 59, doi: 10.3847/1538-4357/aab55b
  * Simpson et al. (2007) Simpson, J. P., Colgan, S. W. J., Cotera, A. S., et al. 2007, ApJ, 670, 1115, doi: 10.1086/522295
  * Smith et al. (2009) Smith, J. D. T., Rudnick, L., Delaney, T., et al. 2009, ApJ, 693, 713, doi: 10.1088/0004-637X/693/1/713
  * Smith et al. (2007a) Smith, J. D. T., Armus, L., Dale, D. A., et al. 2007a, PASP, 119, 1133, doi: 10.1086/522634
  * Smith et al. (2007b) Smith, J. D. T., Draine, B. T., Dale, D. A., et al. 2007b, ApJ, 656, 770, doi: 10.1086/510549
  * Sokal et al. (2016) Sokal, K. R., Johnson, K. E., Indebetouw, R., & Massey, P. 2016, ApJ, 826, 194, doi: 10.3847/0004-637X/826/2/194
  * Stevens et al. (1992) Stevens, I. R., Blondin, J. M., & Pollock, A. M. T. 1992, ApJ, 386, 265, doi: 10.1086/171013
  * Stock & Barlow (2014) Stock, D. J., & Barlow, M. J. 2014, MNRAS, 441, 3065, doi: 10.1093/mnras/stu724
  * Tchernyshyov et al. (2015) Tchernyshyov, K., Meixner, M., Seale, J., et al. 2015, ApJ, 811, 78, doi: 10.1088/0004-637X/811/2/78
  * Thuan & Izotov (2005) Thuan, T. X., & Izotov, Y. I. 2005, ApJS, 161, 240, doi: 10.1086/491657
  * Toalá et al. (2015) Toalá, J. A., Guerrero, M. A., Ramos-Larios, G., & Guzmán, V. 2015, A&A, 578, A66, doi: 10.1051/0004-6361/201525706
  * Toalá et al. (2017) Toalá, J. A., Marston, A. P., Guerrero, M. A., Chu, Y. H., & Gruendl, R. A. 2017, ApJ, 846, 76, doi: 10.3847/1538-4357/aa8554
  * Todt et al. (2015) Todt, H., Sander, A., Hainich, R., et al. 2015, A&A, 579, A75, doi: 10.1051/0004-6361/201526253
  * Tokuda et al. (2021) Tokuda, K., Kondo, H., Ohno, T., et al. 2021, ApJ, 922, 171, doi: 10.3847/1538-4357/ac1ff4
  * Tuthill et al. (1999) Tuthill, P. G., Monnier, J. D., & Danchi, W. C. 1999, Nature, 398, 487, doi: 10.1038/19033
  * Vermeij & van der Hulst (2002) Vermeij, R., & van der Hulst, J. M. 2002, A&A, 391, 1081, doi: 10.1051/0004-6361:20020864
  * Vink et al. (2001) Vink, J. S., de Koter, A., & Lamers, H. J. G. L. M. 2001, A&A, 369, 574, doi: 10.1051/0004-6361:20010127
  * Weaver et al. (1977) Weaver, R., McCray, R., Castor, J., Shapiro, P., & Moore, R. 1977, ApJ, 218, 377, doi: 10.1086/155692
  * Wisniewski & Bjorkman (2006) Wisniewski, J. P., & Bjorkman, K. S. 2006, ApJ, 652, 458, doi: 10.1086/507260
  * Wolfire et al. (1995) Wolfire, M. G., Hollenbach, D., McKee, C. F., Tielens, A. G. G. M., & Bakes, E. L. O. 1995, ApJ, 443, 152, doi: 10.1086/175510
  * Wolfire et al. (2003) Wolfire, M. G., McKee, C. F., Hollenbach, D., & Tielens, A. G. G. M. 2003, ApJ, 587, 278, doi: 10.1086/368016
  * Wu et al. (2006) Wu, Y., Charmandaris, V., Hao, L., et al. 2006, ApJ, 639, 157, doi: 10.1086/499226

##  Appendix A PSF Matching Procedure

We follow the procedure described by Aniano et al. (2011) to create custom
convolution kernels that will transform an image with a narrower PSF to a
broader PSF.

Mathematically, we generate a kernel, $K$, that when convolved transforms
$\Psi_{B}$, the narrower PSF, to $\Psi_{A}$, the broader PSF:

$\Psi_{A}=\Psi_{B}\ast K(B\rightarrow A).$ (A1)

To calculate the convolution kernel, $K$, we take the Fourier Transform (FT)
and solve:

$K(B\rightarrow
A)=\rm{FT^{-1}}\bigg{(}\rm{FT}(\Psi_{B})\times\frac{1}{\rm{FT}(\Psi_{A})}\bigg{)}.$
(A2)

Computing this kernel numerically while ensuring its accuracy and stability
requires several steps outlined in Aniano et al. (2011) and are summarized
below for completeness.

####  A.0.1 Preparing the PSFs

The IRS PSFs are from version 2.0 of the Spitzer Tiny Tim
software222https://irsa.ipac.caltech.edu/data/SPITZER/docs/dataanalysistools/tools/contributed/general/stinytim/.
We create model PSFs for the full wavelength range of IRS, 5 µm - 38 µm. We
note that these are simulated PSFs and there are some differences between the
predicted PSF for IRS and the measured PSF from calibration stars, but these
differences are negligible for our application (Pereira-Santaella et al.,
2010). The PSFs range in pixel scale and image size, depending on the module.
In order to calculate the convolution kernel, $K$, we regrid all PSFs to a
common pixel scale of $0.2\arcsec$ and an image size of $1859\ \times\ 1859$
pixels. The PSF matching algorithm outlined in Aniano et al. (2011) works best
when the PSFs are rotationally symmetric. We circularize the PSFs by rotating
the PSF 14 times and averaging after each rotation in order to make the final
PSF invariant for any angle that is a multiple of $360^{o}/2^{14}=0.022^{o}$.

####  A.0.2 Creating the convolution kernel

We compute the Fourier Transform (FT) of the narrower PSF, $\Psi_{B}$, and the
broader PSF, $\Psi_{A}$. We use the PYTHON package FFTPACK for all FT
calculations. The FT of the target PSF is in the denominator (see Equation A4)
and will amplify any small high-frequency components of the FT when creating
the convolution kernel, $K$. We account for this by introducing a low-pass
filter, $f$, in kernel construction:

$f(k)=\begin{cases}1&\mathrm{for}\ k\leq k_{L}\\\
\frac{1}{2}\times\bigg{[}1+\cos\big{(}\pi\times\frac{k-k_{L}}{k_{H}-k_{L}}\big{)}\bigg{]}&\mathrm{for}\
k_{L}\leq k\leq k_{H}\\\ 0&\mathrm{for}\ k_{H}\leq k\end{cases}$ (A3)

where k is the spatial frequency in the Fourier domain, $k_{H}$ is the high
frequency cutoff for the filter, and $k_{L}=0.7\times k_{H}$ is the low-
frequency cutoff. We create the convolution kernel by modifying Equation A4
with the low-pass filter:

$K(B\rightarrow
A)=\rm{{FT^{-1}}\bigg{(}\rm{FT}(\Psi_{B})\times\frac{1}{\rm{FT}(\Psi_{A})}}\times
f(k)\bigg{)}.$ (A4)

repeating the process by looping over the narrower PSFs for $\Psi_{B}$.
Lastly, we normalize the kernel to unity to ensure flux conservation.

####  A.0.3 Testing the kernel

We test the convolution kernel for accuracy and stability through the metrics
defined by Aniano et al. (2011). The convolution of the narrower PSF,
$\Psi_{B}$, and the kernel, $K(B\rightarrow A)$, should reproduce the broader
PSF $\Psi_{A}$ (see Equation A1). Further, kernels with many negative values,
which redistribute flux, can be unstable. Aniano et al. (2011) define two
parameters, $D$ (their equation 20), the integral of the difference between
the convolved PSF and the target PSF, and $W^{-}$, a sum of the negative
values in the kernel, to evaluate the kernels. We calculate these metrics and
find very low values for $D$ ($<$ 0.01 for all kernels) and satisfactory
values for $W^{-}$ (ranges between 0.2 - 1.0). Aniano et al. (2011) report
that $W^{-}<1.2$ are stable. Further, the value chosen for the filtering
function $k_{H}$ affects these metrics, where an increase in $k$ leads to a
lower $D$ but a higher $W^{-}$. When constructing kernels for each wavelength,
we use a range of $k_{H}$ values and select a $k_{H}$ with the minimum of both
$W^{-}$ and $D$. The $k_{H}$ filtering value ranges from $0.7-1.0$ for the
final convolution kernels.

####  A.0.4 The final PSF matched images

We produce two sets of PSF matched line images: one that provides maximum
uniform resolution by matching to the PSF of the [Si II] transition (our
longest wavelength line), and one that matches everything to a 12″ Gaussian
for ease of comparison with other data. In our analysis we use the 12″
Gaussian PSF. The [Si II] 34.8 µm line has the broadest PSF out of the
spectral lines detected with Spitzer IRS with a FWHM of 7.6″. We match each
narrower IRS PSF to the PSF of the [Si II] line. Then, we create and apply a
convolution kernel to go from the [Si II] PSF to a 2D Gaussian PSF. We find
that a 2D Gaussian with a FWHM of 12″ best matches the [Si II] PSF while
balancing kernel stability.

Each emission line image is convolved with its custom kernel using the ASTROPY
package CONVOLUTION. We then regrid each image to the coordinates of the [Si
II] image to match its pixel scale and to simplify the analysis. For the
analysis we also mask non-physical bright artifacts and the emission from the
SNR E0102-72. We use 12″ Gaussian PSF-matched final images for the analysis in
this work. Note that Figure 1 presents the original resolution images to show
the maximum detail in the morphology of N76.

In addition to creating the PSF matched line maps that are described in
$\S$2.1, we also produce a PSF matched cube for the full IRS spectrum,
including all orders: SL1, SL2, LL2, LL1. We follow the same procedure on each
frame of the original data cubes. A spectrum of the central position of the
N76 nebula from the PSF-matched cube is presented in Figure 2.

We PSF match the IRS Spitzer data, but not the Herschel PACS images (see
$\S$2.2). The emission line images from PACS have larger PSFs that would
degrade too much the resolution of the IRS maps if we were to match each image
to the largest Herschel PSF.

##  Appendix B Effect of asymmetries in N76

In order to compare the Cloudy models directly to various emission lines in
N76, we need to assume spherical symmetry. However, it is clear that the
images of the emission lines show spatial asymmetries in Figure 1. The eastern
side of N76 is brighter in almost every IR line and the western edge is less
well-defined, with more diffuse, extended features. To explore how these
asymmetries affect the Cloudy modeling, we show here the results for a narrow
strip that intersects with the AB7 and designate it the central point. We then
sample the surface brightness to the left (eastern edge) and right (western
edge) of the strip and report it in Figure 10.

The asymmetric profiles show that the western edge prefers a lower density
model in almost all cases (except for [O I], [C II], and [Ne II]). The
idealized symmetric model presented in the radial profiles in Figure 3 is
simply an average of these two sides. The difference between these densities
is relatively low, however, with the higher density side preferring a model
only $4\ \mbox{$\mathrm{cm^{-3}}$}$ larger. This is typically a $\pm 20-30\%$
change in the result of the spherically symmetric modeling.

Figure 10: The surface brightness distribution of emission lines taken as a
strip across the nebula with AB7 in the center and the eastern and western
edges on the left and right hand of the plot, respectively. This shows the
asymmetry of N76, where the eastern edge is brighter than the western side.
The Cloudy models prefer a slightly lower density for the western side.