# Investigating the Drivers of Electron Temperature Variations in HII Regions
with Keck-KCWI and VLT-MUSE

Ryan J. Rickards Vaught Center for Astrophysics and Space Sciences, Department
of Physics, University of California, San Diego, 9500 Gilman Dr., La Jolla, CA
92093, USA Ryan J. Rickards Vaught<EMAIL_ADDRESS>Karin M. Sandstrom
Department of Astronomy & Astrophysics, University of California, San Diego,
9500 Gilman Dr., La Jolla, CA 92093, USA  Francesco Belfiore INAF –
Osservatorio Astrofisico di Arcetri, Largo E. Fermi 5, I-50157 Firenze, Italy
Kathryn Kreckel Astronomisches Rechen-Institut, Zentrum für Astronomie der
Universität Heidelberg, Mönchhofstr. 12-14, D-69120 Heidelberg, Germany  J.
Eduardo Méndez-Delgado Astronomisches Rechen-Institut, Zentrum für Astronomie
der Universität Heidelberg, Mönchhofstr. 12-14, D-69120 Heidelberg, Germany
Eric Emsellem European Southern Observatory, Karl-Schwarzschild Straße 2,
D-85748 Garching bei München, Germany  Univ Lyon, Univ Lyon1, ENS de Lyon,
CNRS, Centre de Recherche Astrophysique de Lyon UMR5574, F-69230 Saint-Genis-
Laval France  Brent Groves International Centre for Radio Astronomy Research,
University of Western Australia, 7 Fairway, Crawley, 6009 WA, Australia
Guillermo A. Blanc Observatories of the Carnegie Institution for Science, 813
Santa Barbara Street, Pasadena, CA 91101, USA  Departamento de Astronomía,
Universidad de Chile, Camino del Observatorio 1515, Las Condes, Santiago,
Chile  Daniel A. Dale Department of Physics & Astronomy, University of
Wyoming, Laramie, WY 82071, USA  Oleg V. Egorov Astronomisches Rechen-
Institut, Zentrum für Astronomie der Universität Heidelberg, Mönchhofstr.
12-14, D-69120 Heidelberg, Germany  Simon C. O. Glover Universität Heidelberg,
Zentrum für Astronomie, Institut für Theoretische Astrophysik, Albert-Ueberle-
Str. 2, 69120, Heidelberg, Germany.  Kathryn Grasha ARC DECRA Fellow  Research
School of Astronomy and Astrophysics, Australian National University,
Canberra, ACT 2611, Australia  ARC Centre of Excellence for All Sky
Astrophysics in 3 Dimensions (ASTRO 3D), Australia  Ralf S. Klessen
Universität Heidelberg, Zentrum für Astronomie, Institut für Theoretische
Astrophysik, Albert-Ueberle-Str. 2, 69120, Heidelberg, Germany.  Universität
Heidelberg, Interdisziplinäres Zentrum für Wissenschaftliches Rechnen, Im
Neuenheimer Feld 205, D-69120 Heidelberg, Germany.  Justus Neumann Max-Planck-
Institut für Astronomie, Königstuhl 17, D-69117 Heidelberg, Germany  Thomas G.
Williams Sub-department of Astrophysics, Department of Physics, University of
Oxford, Keble Road, Oxford OX1 3RH, UK

(Received August 29, 2023)

###### Abstract

H ii region electron temperatures are a critical ingredient in metallicity
determinations and recent observations reveal systematic variations in the
temperatures measured using different ions. We present electron temperatures
($T_{e}$) measured using the optical auroral lines ([N II]$\lambda 5756$, [O
II]$\lambda\lambda 7320,7330$, [S II]$\lambda\lambda 4069,4076$, [O
III]$\lambda 4363$, and [S III]$\lambda 6312$) for a sample of H ii regions in
seven nearby galaxies. We use observations from the Physics at High Angular
resolution in Nearby Galaxies survey (PHANGS) obtained with integral field
spectrographs on Keck (Keck Cosmic Web Imager; KCWI) and the Very Large
Telescope (Multi-Unit Spectroscopic Explorer; MUSE). We compare the different
$T_{e}$ measurements with H ii region and interstellar medium environmental
properties such as electron density, ionization parameter, molecular gas
velocity dispersion, and stellar association/cluster mass and age obtained
from PHANGS. We find that the temperatures from [O II] and [S II] are likely
over-estimated due to the presence of electron density inhomogeneities in H ii
regions. We observe that differences between [N II] and [S III] temperatures
are weakly correlated with stellar association mass and molecular gas velocity
dispersion. We measure high [O III] temperatures in a subset of regions with
high molecular gas velocity dispersion and low ionization parameter, which may
be explained by the presence of low-velocity shocks. Of all the temperatures
considered, the $T_{\rm{e}}$–$T_{\rm{e}}$ between [N II] and [S III]
temperatures have the lowest observed scatter and generally follow predictions
from photoionization modeling, which suggests that these tracers reflect H ii
region temperatures across the various ionization zones better than [O II], [S
II], and [O III].

##  1 Introduction

The characterization of abundance variations within galaxies provides insight
into the physical processes that drive galaxy and chemical evolution. A
galaxy’s gas-phase metal abundance (i.e., metallicity) reflects the history of
chemical enrichment from stars and the net balance of gas flows (mixing,
outflows, inflows of pristine material, etc.). In addition, the metallicity of
ISM gas directly controls its cooling and other important ISM physics (Draine,
2011; Peimbert et al., 2017).

The distribution of gas-phase metals in a galaxy is commonly traced by the
abundance of oxygen, nitrogen, sulfur, and other metals using the emission
from ionized gas located inside H ii regions (e.g. Kennicutt & Garnett, 1996;
Bresolin et al., 2012; Hernandez et al., 2013; Ho et al., 2017; Kreckel et
al., 2019, 2020; van Loon et al., 2021; Grasha et al., 2022). There are
several indirect methods calibrated using strong optical emission lines to
derive an estimate of the H ii region metallicity (e.g. Kewley & Ellison,
2008; Blanc et al., 2015). A “direct” measure of an H ii region metallicity
requires knowledge of the electron temperature ($T_{\rm{e}}$) of the gas. Due
to their exponential dependence on $T_{\rm{e}}$, one of the ways to infer
electron temperature is through the auroral-to-nebular line ratios of
collisionally excited lines (CEL; Peimbert, 1967; Osterbrock & Ferland, 2006;
Peimbert et al., 2017). Nebular and auroral lines originate from different
excited states of ions. Auroral lines are from higher energy levels, but are
still accessible for collisional excitation in a $T\sim 10^{4}$ K gas. If the
density of the gas is below the auroral and nebular line critical densities
(i.e. when collisional de-excitation is negligible), then the auroral-to-
nebular line ratio is sensitive to the electron temperature (e.g., Osterbrock
& Ferland, 2006). Given that the excitations to the auroral level are only
accessible to electrons of higher energy, auroral line emission can be $>100$
times weaker than nebular lines (Kennicutt et al., 2003; Esteban et al., 2004;
Berg et al., 2020). One alternative way to to measure $T_{\rm{e}}$ include the
ratio between recombination line (RL) emission from H and other species, and
which exhibit linear sensitivity to $T_{\rm{e}}$ (Peimbert, 1967; Osterbrock &
Ferland, 2006; Peimbert et al., 2017). However, optical RLs useful for use as
temperature diagnostics are typically reserved for deep, high S/N spectra, as
RL emission is typically much fainter than the emission from auroral lines.

For ions with optical auroral lines studied in this work, we can measure the
temperatures for each ion using the following line ratios:

$T_{\rm{e},\text{[{O}{III}]}}$ $\rightarrow$ [O III]$\lambda
4363$/$\lambda\lambda 4959,5007$,

$T_{\rm{e},\text{[{O}{II}]}}$ $\rightarrow$ [O II]$\lambda\lambda
7320,7330$/$\lambda\lambda 3726,3729$111[O II]$\lambda\lambda 7320,7330$ is an
unresolved quadruplet with transitions at $\lambda 7319$, $\lambda 7320$,
$\lambda 7330$ and $\lambda 7331$ Å.,

$T_{\rm{e},\text{[{S}{III}]}}$ $\rightarrow$ [S III]$\lambda
6312$/$\lambda\lambda 9069,9532$,

$T_{\rm{e},\text{[{S}{II}]}}$ $\rightarrow$ [S II]$\lambda\lambda
4069,4076$/$\lambda\lambda 6716,6731$,

$T_{\rm{e},\text{[{N}{II}]}}$ $\rightarrow$ [N II]$\lambda
5756$/$\lambda\lambda 6548,6584$.

The O+, N+, S+ ions require energies of $13.6$ eV, $14.5$ eV, and $10.3$ eV to
be produced while S++ and O++ require energies 23 eV and 35 eV, respectively.

Several effects play competing roles in determining the ionization and
temperature structure of H ii regions. These include a radially decreasing
intensity and hardening of the radiation field (photons closest to 13.6 eV are
absorbed first) as well as a change in the ions which dominate gas cooling,
and therefore the cooling efficiency (Stasińska, 1980; Garnett, 1992). Because
of the varying degree of ionization within an H ii region, a model of three
ionization zones—low-, intermediate- and high—is often used to describe them.
Because each ionization zone could theoretically have different temperatures,
this further stresses the importance of observing multiple auroral lines and
developing temperature priorities for use in accurately determining abundances
(e.g. Berg et al., 2015, 2020).

Observing the full set of optical auroral lines in an H ii region can be
challenging. In addition to the large wavelength range needed, $\sim$
3700–10000 Å, some auroral lines are weaker than others depending on the
metallicity and temperature of the gas. Because of these challenges, it is
very important to establish temperature relationships that allow us to infer
the conditions of a certain ionization zone from the others. Photoionization
modeling (e.g. Garnett, 1992; Vale Asari et al., 2016) has been used to derive
temperature relationships, but standard models consider only simple geometries
and homogeneous physical and ionization conditions, that might not be suitable
for more complex regions potentially affected by shocks, stellar feedback or
other mechanism that produce density or temperature inhomogeneities (Peimbert,
1967; Peimbert et al., 1991; Binette et al., 2012; Berg et al., 2015, 2020;
Arellano-Córdova & Rodríguez, 2020; Nicholls et al., 2020; Méndez-Delgado et
al., 2023b, a).

In the presence of temperature fluctuations, the exponential dependence of CEL
strengths on temperature will bias auroral-to-nebular temperatures towards
higher values than the true average (Peimbert, 1967; Peimbert & Costero,
1969). Such inhomogeneities may be related to the presence of turbulence,
density structure, and shocks associated with either stellar winds or
radiation–pressure driven expansion. If the sources of temperature
inhomogeneities are confined to the central part of the nebula, the effects
that these phenomena have on temperature may primarily affect only the high
ionization zone. This has been suggested by Méndez-Delgado et al. (2023a) who
presented evidence for temperature inhomogeneities affecting only the highly
ionized gas traced by [O III]. In a sample of Galactic and extra-galactic H ii
regions, they observed that differences between [O III] and [N II]
temperatures correlated with the degree of deviation from the average
temperature measured using faint O II recombination line emission.
Furthermore, a strong correlation between the O II recombination and [N II]
temperatures observed by Méndez-Delgado et al. (2023a) implies that
temperatures inferred from the [N II] auroral line accurately measures the
average $T_{\mathrm{e}}$ of the low-ionization zone (Méndez-Delgado et al.,
2023a, b).

Due to the importance of obtaining accurate temperatures for precise
abundances, significant effort has been devoted to advancing our understanding
of the temperatures of different H ii region ionization zones. For example,
previous works have found that the scatter between temperatures of different
ionization zones may be correlated with other properties of the gas such as
the ionization parameter and metallicity (Berg et al., 2015; Arellano-Córdova
& Rodríguez, 2020; Berg et al., 2020; Yates et al., 2020).

A number of important questions about electron temperatures in H ii regions
remain. Past studies (e.g. Kennicutt et al., 2003; Rosolowsky & Simon, 2008;
Berg et al., 2015, 2020; Wei et al., 2020) have typically targeted the
brightest H ii regions in order to maximize the probability of detection of
auroral lines. This procedure, however, may have introduced selection biases
towards the brightest regions or those with particular physical conditions. In
addition, due to observational challenges (e.g., wide wavelength range needed
to cover the [S III] lines) it remains uncertain how temperature
inhomogeneities may impact the intermediate ionization zone temperatures
traced by [S III] (Méndez-Delgado et al., 2023b).

To explore these questions, we use deep 3600–9500 Å IFU mapping to measure the
set of optical auroral lines and nebular lines for a sample of H ii regions.
We use observations obtained from the Keck Cosmic Web Imager (KCWI, Morrissey
et al., 2018) and Multi-Unit Spectroscopic Explorer (MUSE, Bacon et al., 2010)
to measure the electron temperature from all 3 ionization zones in H ii
regions in nearby galaxies. In Section 2 we present our sample galaxies as
well as primary and supplemental observations. In Section 3 we discuss the
reduction of the KCWI data. We assess the quality of the KCWI mosaics in
Section 1. We construct our H ii region sample in Section 4. We present the
auroral line measurements in Section 5. We derive H ii region properties from
nebular diagnostics and from ALMA and HST data in Section 6. The results and
discussion are presented in Section 7 and 8.

##  2 Observations

The analysis presented here makes joint use of multi-wavelength observations
of seven galaxies obtained with Keck-KCWI, VLT-MUSE, Atacama Large
Millimeter/submillimeter Array (ALMA), and the Hubble Space Telescope (HST).

###  2.1 Sample Selection

The seven galaxies in this analysis are drawn from the PHANGS-MUSE sample
(Emsellem et al., 2022). To date, 90 galaxies make up the full PHANGS
sample222http://phangs.org/ (Leroy et al., 2021a), and 19 have been observed
by MUSE. In order to be observed with KCWI in the northern hemisphere, we
selected the seven target galaxies from a subset of PHANGS-MUSE galaxies with
declination, $\delta>-30^{\circ}$. Table 1 presents general properties of
these galaxies, including distances, masses, sizes, and the angular resolution
of the MUSE data.

###  2.2 Keck Cosmic Web Imager

We observed each galaxy using KCWI on the Keck II telescope with multiple
pointings taken over several nights between the years 2017 and 2021. Clear
conditions were present for the majority of observations, except for the
nights of October 16 and 17, 2018, which suffered from variable cloud
coverage. These poor conditions primarily affect the observations of NGC 628.
The instrument was configured with the “Large” slicer and BL grating centered
at 4600 Å. The usable spectral range afforded by this configuration is
3650–5550 Å with a spectral resolution $R\sim 900$, corresponding to a full
width at half maximum (FWHM) $\sim 5.1$ Å (or $\sim$ 300 km s-1) at the
central wavelength. The Large slicer has an angular slice width of 1.35″. The
field of view (FoV) using the Large slicer and BL grating is 33″ perpendicular
and 20.4″ parallel to the slicer.

Because the FoV is small compared with the large angular size of each galaxy,
we observed each galaxy over multiple fields. Most fields were observed two
times using 1200 s (i.e., 20 min) integration times. The only exceptions were:
all fields in NGC 3627, which were observed five times each with 120 s (2 min)
integration times; field 17 in NGC 628 which was observed 3 times using 1200
s; and field 2 in NGC 5068 and field 5 in NGC 1385, both having only a single
observation of 1200 s. A half slice width, or 0.675″, dither was applied
between each exposure. Immediately following the second exposure, we observed
an off-galaxy region, selected to be free of extended emission and/or bright
sources, in order to measure a sky spectrum close in time to the observations.
These sky frames were used for sky subtraction during data reduction. We
summarize the number of fields, exposure times, and dates in Table 7 of
Appendix A. The full data reduction is outlined in Section 3.

###  2.3 Multi-Unit Spectroscopic Explorer

MUSE observations of these galaxies come from the PHANGS-MUSE survey (Emsellem
et al., 2022). MUSE covers wavelengths between 4800–9500 Å. Taken in
combination KCWI and MUSE span the full optical spectrum. The full details of
the MUSE data reduction and data products are presented in Emsellem et al.
(2022), and we provide a brief overview here. The PHANGS-MUSE program observed
19 galaxies using 168 individual 1${}^{\prime}\times 1^{\prime}$ pointings.
The median spatial resolution across all pointings is $\sim 50$ pc (or $\sim$
0.80″) with a typical spectral FWHM of $\sim 2.5$ Å (but varying with
wavelength). The data were reduced using the
pymusepipe333https://pypi.org/project/pymusepipe/ and spectral fitting and
analysis was performed using the Data Analysis Pipeline
444https://gitlab.com/francbelf/ifu-pipeline packages described in Emsellem et
al. (2022). The individual MUSE pointings were homogenized to a uniform
Gaussian point-spread-function (PSF) with FWHM set to the largest FWHM
measured for each target, resulting in “convolved and optimized” (COPT)
mosaics. The PSFs of the COPT mosaics are listed in Table 1. We use the COPT
mosaics in the following work.

###  2.4 ALMA

Our analysis makes use of Atacama Large Millimeter/submillimeter Array (ALMA)
data obtained as part of PHANGS–ALMA (Leroy et al., 2021a). PHANGS-ALMA
observed the $J=2-1$ rotational transition of 12CO, hereafter CO, for $90$
galaxies. The details of the data reduction are described in Leroy et al.
(2021b). We make use of the ALMA datacubes with combined CO measurements from
the 12m and 7m arrays plus Total Power (12m+7m+TP). The nominal angular
resolution of 12m+7m+TP observations is $\sim 1.3\arcsec$, similar to the
resolution of both KCWI and MUSE. The velocity resolution is $2.5$ km s-1.

###  2.5 HST

The PHANGS-HST survey (Lee et al., 2022) observed555Lee et al. (2022) use of
previous NGC 628 HST imaging obtained as part of the Legacy ExtraGalactic
Ultraviolet Survey (Legus; Calzetti et al., 2015). our target galaxies using 5
HST filters: F275W (NUV), F336W (U), F438W (B), F555W (V), F814W (I). Of the
high-level data products produced from this dataset, we make use of compact
star cluster catalogs (Thilker et al., 2022) and multi-scale stellar
association catalogs (Larson et al., 2023). The association catalog identifies
sources using both the V and NUV filters and has been convolved to several
physical resolutions (8, 16, 32 and 64 pc, respectively). Following
Scheuermann et al. (2023) we used the NUV selected, 32 pc catalogs.

Table 1: Properties of the PHANGS-MUSE galaxies observed with KCWI.

Name | DistanceaaFrom the compilation of Anand et al. (2021) | log${}_{10}(M_{*})$ bbDerived by Leroy et al. (2021a), using GALEX UV and WISE IR photometry. |  $R_{25}$ccFrom LEDA Makarov et al. (2014). | PSFMUSEddThe FWHM of the Gaussian PSF for the homogenized COPT mosaic from PHANGS-MUSE (Emsellem et al., 2022). | PSFKCWIeeThe average FWHM of the KCWI PSF for the set of a galaxy’s observed pointings.
---|---|---|---|---|---
| [Mpc] | [$M_{\odot}$] | [Arcmin] | [Arcsec] | [Arcsec]
NGC 628 | 9.8$\pm$0 | 10.34$\pm$0.1 | 4.9 | 0.92 | 2.0 $\pm$ 0 |  |
NGC 1087 | 15.9$\pm$2 | 9.93$\pm$0.1 | 1.5 | 0.92 | 1.2 $\pm$ 0 |  |
NGC 1300 | 19.0$\pm$2 | 10.62$\pm$0.1 | 3.0 | 0.89* | 1.3 $\pm$ 0 |  |
NGC 1385 | 17.2$\pm$2 | 9.98$\pm$0.1 | 1.7 | 0.77* | 1.3 $\pm$ 0 |  |
NGC 2835 | 12.2$\pm$0 | 10.00$\pm$0.1 | 3.2 | 1.15 | 1.4 $\pm$ 0 |  |
NGC 3627 | 11.3$\pm$0 | 10.83$\pm$0.1 | 5.1 | 1.05 | 1.1 $\pm$ 0 |  |
NGC 5068 | 5.2$\pm$0 | 9.40$\pm$0.1 | 3.7 | 1.04 | 1.5 $\pm$ 0 |  |

**footnotetext: Denotes galaxies observed with MUSE using ground based
adaptive optics.

##  3 KCWI Data Reduction

The KCWI observations were reduced using the Version 1.0.1 Python
implementation of the KCWI Data Extraction and Reduction Pipeline
(KDRP)666KCWI DRP-Python. It was built using the Keck Data Reduction Pipeline
Framework package777KeckDRPFramework and is a port of the initial IDL
reduction pipeline888KCWI DRP-IDL (Morrissey et al., 2018). The pipeline
performs basic CCD reduction including bias and over-scan subtraction, gain
correction, cosmic ray removal, dark and scattered light subtraction as well
as a flat field correction.

Following these basic reductions, the KDRP used the continuum bar and
Thorium/Argon arc lamp images to generate geometric and wavelength solutions
to convert each 2D science image into a spectral datacube. The accuracy of the
wavelength solutions were similar across all the observation nights. The
average RMS for the derived wavelength solutions was 0.1 Å.

We derived an inverse sensitivity curve to flux calibrate each datacube from
standard star observations. The measured standard deviation between all of the
derived inverse sensitivity curves was $\sim 9\%$ at $\lambda=4600$Å. The
maximum standard deviation within the wavelength range containing the lines
used in this analysis is $\sim 11\%$. Details of each standard star
observation can be found in Table 8 of Appendix A. After flux calibration,
each datacube was corrected for differential atmospheric refraction.

Because the instrument FoV is much smaller than each galaxy, our images
contained no sky pixels. To perform sky subtraction, we observed dedicated sky
positions interspersed between science observations. We assigned the sky frame
closest in time to each science observation to be used for sky subtraction
according to the instructions in the KDRP documentation999https://kcwi-
drp.readthedocs.io/en/latest/sky_subtraction.html. The KDRP then performed sky
subtraction using our preferred frames. The sky in all pixels was averaged
together and scaled by the ratio of science-to-sky exposure time to estimate
the sky observed in the “on” position. The final data products output by the
pipeline include flux calibrated science and $1\sigma$ uncertainty cubes, as
well as a bad-pixel mask cube.

###  3.1 Image Re-Projection

Next we constructed mosaics from the individual KCWI pointing datacubes. The
steps involved included image registration, matching the flux calibration to
MUSE, and image co-addition. We also compared the absolute flux calibration of
the final KCWI mosaics to MUSE and SDSS

The datacubes output by the KDRP have rectangular pixels with pixel-scale
1.35″ $\times$ 0.29″. We reprojected the cubes onto a square pixel grid using
the astronomical mosaic software Montage101010See
http://montage.ipac.caltech.edu/. Prior to running Montage, we converted the
KCWI data to surface brightness units (erg s-1 cm-2 Å-1 sr-1) by dividing the
flux per pixel by the pixel solid angle in steradians. The reprojected images
have a final square pixel grid with a uniform pixel-scale of 0.29″$\times$
0.29″. We validated the flux-conservation in our data before and after
reprojection.

###  3.2 Image Registration

To co-add each galaxy’s set of cubes into a spectral mosaic, we placed each
cube onto a common world coordinate system (WCS). It is typical to perform
image registration by matching the location of foreground stars or background
galaxies to known positions found in catalogs. However, in our case most
individual fields did not contain a sufficient number of bright point sources.
Instead, we performed image registration by maximizing the cross-correlation
between individual KCWI fields and overlapping MUSE data in order to match the
KCWI pointing astrometry to MUSE. The astrometry of the MUSE galaxy mosaics
were validated against wide-field broadband imaging and stellar positions from
the Gaia DR1 as described in Emsellem et al. (2022). The MUSE astrometry, when
compared to broadband imaging, exhibited better than 100 milli-arcseconds RMS.
in both R.A. and Dec.

To cross-correlate KCWI and MUSE, we created synthetic photometry ($P_{S}$)
images from spectral regions where the wavelength coverage of KCWI and MUSE
overlap. Because there is some saturation in H$\beta$ and [O III] at the
brightest locations in the KCWI data, we masked out those lines in both cubes
to avoid any issues with the comparison between the two data cubes.

In order to determine the optimal astrometric shifts to apply to each KCWI
science frame, we utilized a two-step grid search operation, first shifting in
1 pixel (or 0.29″) steps $\pm$ 17 pixels (or $5\arcsec$) from the center of
the KCWI pointing in order to find the optimal R.A. and Dec. offsets which
maximize the correlation of the KCWI and corresponding MUSE $P_{S}$ images.
After locating 1st-pass shifts, we performed a secondary grid search using
finer 1/2 pixel increments over a smaller range ($\pm 1\arcsec$ from the image
center). The 0.5 pixel sampling corresponds to 0.145″ which is less than the
MUSE pixel scale of 0.20″ but also corresponds to 1/10 the typical KCWI FWHM
which is equal to 1.2″. Across the galaxy sample, the final offsets correspond
to correlation coefficients $>$ 0.9 between KCWI and MUSE.

###  3.3 Matching the MUSE Flux Calibration

In order to correct for any additive and/or multiplicative offsets between the
MUSE and KCWI flux calibrations, we compared the surface brightness (SB)
calculated in apertures in overlapping position and wavelength. To do this we
made use of the $P_{S}$ images, described in Section 3.2, and measured the
surface brightness inside a number of 3″ radius apertures located at randomly
drawn positions inside the combined KCWI and MUSE coverage. The aperture size
was chosen to be large enough to minimize effects arising from the different
PSFs of KCWI and MUSE. We determined the best fit line to the measured
$SB_{\mathrm{KCWI}}$ vs. $SB_{\mathrm{MUSE}}$ relationship, where the slope,
$m$, and $y$-intercept, $b$, reflect the multiplicative and additive offset
between the KCWI and MUSE flux calibration. We applied the correction by
multiplying the KCWI datacubes by $m$ and adding $b$ to the full spectrum in
each pixel. The average multiplicative and additive offsets were $m=1.03\pm
0.02$ and $b=-7.6\pm 1.7\times 10^{-20}\ \mathrm{erg}\ \mathrm{s}^{-1}\
\mathrm{cm}^{-2}\ \text{\AA}^{-1}$. This $\sim 3\%$ deviation from a 1–1 slope
and low level of additive offset show that the calibrations were already in
good agreement.

###  3.4 PSF of Individual KCWI Fields

Because the MUSE mosaics have been homogenized to a uniform Gaussian PSF, we
have an image with a known (parameterized) PSF. This is advantageous, as we
can directly measure the KCWI PSF per pointing using cross-convolution,
following the steps outlined in Emsellem et al. (2022). We briefly summarize
the procedure here. 1) We produced $P_{S}$ images of both the MUSE mosaic and
the individual KCWI pointings. 2) We re-projected the MUSE cutout onto the
KCWI 0.29″$\times$ 0.29″ pixel grid. 3) We convolved the KCWI pointing with a
2D Gaussian Kernel with PSF equal to the reference MUSE PSF. 4) In an
iterative manner, the reference MUSE image is then convolved with a 2D
Gaussian Kernel, PSFk, where the PSFk represents the KCWI pointing PSF which
is a free-parameter. We then varied the FWHM of this kernel until we maximized
the correlation between the KCWI pointing and the reference MUSE image. For
the set of images observed for each galaxy, we present the average and
standard deviation of the measured PSFs in Table 1. The average PSF across the
galaxy sample is $1.4^{\prime\prime}\pm 0.2^{\prime\prime}$ which is
consistent with the PSFs measured using the Standard Star observations
presented in Appendix A. We chose not to homogenize the PSF of the KCWI data.
To do so, would mean convolving the KCWI data to the largest observed PSF. In
turn, this would increase the mismatch in resolution between KCWI and MUSE as
well as lead to larger H ii region boundaries. The larger regions have higher
susceptibility of contamination from the diffuse ionized gas.

###  3.5 Image Co-Addition

After KCWI datacubes had been aligned and flux calibrated to match the MUSE
mosaics, the KCWI datacubes were co-added to create KCWI galaxy mosaics. The
image co-addition was performed with Montage. The mAddCube call to Montage
initiates the co-addition. The co-addition is performed by stacking the
aligned images, according to the an output WCS determined by Montage, and then
taking the average value between any overlapping pixels. Pseudo $g$-band
images for the final mosaics of each galaxy are shown in Figure 1.

Figure 1: Pseudo $g$-band images of the KCWI mosaics. A 10″ scale bar is shown
in the bottom of each image.

###  3.6 Absolute calibration of KCWI Compared to MUSE

We assesed the absolute calibration of the co-added KCWI mosaics by comparing
the SB’s between the KCWI and MUSE mosaics. Comparisons of the MUSE mosaics
with SDSS imaging in the $r$-band described in Emsellem et al. (2022) showed
that the MUSE absolute flux calibration is consistent with SDSS calibration
within the instrument uncertainty of $0.06$ mag (Padmanabhan et al., 2008). We
calculated the SB inside $r=3\arcsec$ apertures, randomly placed in the KCWI
coverage. The SB offsets between the KCWI and MUSE mosaics are shown in Figure
2. The resulting offsets, summarized in Table 2, reveal acceptable agreement
between the MUSE and KCWI absolute calibration. The average percent SB offset
with respect to the SDSS imaging,
$\Delta\mathrm{SB}/\mathrm{SB}_{\mathrm{SDSS}}$, is between $-1.1$% and $0.7$%
with a median value across all galaxies of $-0.1$% $\pm$ 4%.

Figure 2: Surface Brightness comparison measured in $3$″ diameter apertures
between overlapping KCWI and MUSE mosaics. In each panel we show the
fractional SB differences between KCWI and MUSE versus the SB of MUSE. The
median fractional offset (black-dashed) is shown relative to the zero line
(black-solid). Across all galaxies, the median offset is $\sim-0.1\%$. Figure
3: Comparison of $g$-band surface brightness measured in $r=3$″ apertures from
synthetic KCWI $g$-band mosaics and SDSS imaging for the 4 KCWI galaxies with
available SDSS data. The surface brightness offset between KCWI and SDSS is
shown versus the SDSS $g$-band surface brightness. The median offset (black-
dashed) is shown relative to the zero line (solid-black). Across the sample,
the offset between the KCWI and SDSS surface brightness is $\sim 1$%.

###  3.7 Absolute calibration of KCWI Compared to SDSS

We have shown agreement between KCWI and MUSE, but this comparison is only an
assessment of the flux calibration in the overlapping wavelength range of KCWI
and MUSE. To assess the absolute flux calibration across a wider wavelength
range, we compared synthetic g-band images of the KCWI mosaics to Sloan
Digital Sky Survey (SDSS, Abazajian et al., 2003) images of the same galaxies.
Only four galaxies with KCWI mosaics have SDSS imaging: NGC 628, NGC 1087, NGC
3627 and NGC 5068. We constructed synthetic g-band images of these galaxies by
convolving the spectrum in each pixel, $F_{\lambda}(x,y)$, with the $g$-band
transmission curve, $T_{g}(\lambda)$, according to the following equation:

$F_{g}(x,y)=\frac{\int F_{\lambda}(x,y)T_{g}(\lambda)d\lambda}{\int
T_{g}(\lambda)d\lambda}.$ (1)

The SDSS imaging is presented in units of nanomaggies or $f_{\nu}=$
3.631$\times 10^{-6}$ Jy. To compare with the KCWI data, with native units of
flux density, $f_{\lambda}$, we converted to flux density with the following
expression (Tokunaga & Vacca, 2005),

$f_{\lambda}=\frac{c}{\lambda_{p}^{2}}f_{\nu},$ (2)

where $\lambda_{p}$ is the pivot wavelength of the band-pass. The pivot
wavelength for the g-band filter is $\lambda_{p}=4702$ Å. The surface
brightness for both KCWI and SDSS are measured inside $3\arcsec$ radius
apertures. The mean and 1$\sigma$ scatter of the measured surface brightness
are shown in Figure 3. We found overall agreement of the absolute calibration
between KCWI and SDSS for the four galaxies. Summarized in Table 2, the median
offset across the sample galaxies,
$\Delta\mathrm{SB}/\mathrm{SB}_{\mathrm{SDSS}}$ range between $-$1.0% and 3.0%
and exhibits scatter between 3% and 10%. The largest scatter appears to be
limited to a single field of NGC 628 that only contains signal from a single,
isolated region, and will have negligible impact on the remaining analysis.
The median offset across all galaxies is $1\pm 7$% which suggests that the
KCWI calibration is in good agreement with SDSS.

Table 2: Flux Calibration Comparisons between KCWI, SDSS and MUSE. Name |  $\mu(\Delta$SB)aaThe median fractional surface brightness offset between KCWI and MUSE in %. |  $\sigma(\Delta$SB)bbThe standard deviation of the fractional surface brightness offset between KCWI and MUSE in %. |  $\mu(\Delta$SBg)ccThe median $g$-band surface brightness offset between KCWI and SDSS. |  $\sigma(\Delta$SBg)ddThe standard deviation of $g$-band surface brightness offset between KCWI and SDSS.
---|---|---|---|---
| [%] | [%] | [%] | [%]
NGC 628 | 0.3 | 6.0 | 0.3 | 10.0 |  |
NGC 1087 | 0.5 | 3.6 | -1.0 | 3.4 |  |
NGC 1300 | 1.5 | 4.0 | \-  | - |  |
NGC 1385 | -0.6 | 4.3 | \-  | - |  |
NGC 2835 | -0.4 | 3.6 | \-  | - |  |
NGC 3627 | 1.7 | 1.6 | -0.05 | 2.6 |  |
NGC 5068 | 0.3 | 6.7 | 2.9 | 8.4 |  |

##  4 H II Region Catalog

In order to assess the emission line properties of H ii regions, we determine
the H ii region location and boundaries using H$\beta$ maps constructed from
the KCWI spectral datacubes and the image segmentation software HIIphot
(Thilker et al., 2000). Although H ii region masks have previously been
constructed from the MUSE H$\alpha$ maps for these galaxies (Kreckel et al.,
2019; Santoro et al., 2022; Groves et al., 2023; Congiu et al., 2023), the
MUSE angular resolution is higher than that of KCWI. Because of this, the H ii
region boundaries derived from MUSE may not fully encapsulate the spatial
extent of the H ii regions observed using KCWI. Furthermore, simply convolving
or reprojecting the MUSE H ii masks to the KCWI resolution or grid would
introduce uncertainty on the boundaries for tightly spaced H ii regions. We
therefore perform H ii region identification directly on the KCWI data.

###  4.1 Construction of H$\beta$ Maps

H$\beta$ is the brightest H I recombination line observed by KCWI, and maps of
this emission for the galaxies will be used to define our H ii regions. The
continuum near and underlying H$\beta$ emission must be removed in order to
accurately map its emission. To remove the continuum we used LZIFU, an
emission line fitting code designed specifically for use with IFUs (Ho et al.,
2016). LZIFU implements and streamlines the penalized pixel-fitting software
(PPXF, Cappellari & Emsellem, 2004; Cappellari, 2017) for using PPXF on IFU
emission line maps. To fit the continuum of the input spectrum LZIFU matches a
series of input single-metallicity, $-1.31<$ [Z/H] $<0.22$, stellar population
models (Vazdekis et al., 2010, MILES) that have been redshifted, and convolved
to match the input spectrum PSF. LZIFU fits Gaussian models at the location of
emission lines. In the output H$\beta$ map, pixels with weak or no H$\beta$
can contain ‘NaN’ values which can be problematic for HIIphot. To avoid these
artifacts, we subtract the stellar continuum from each pixel’s spectrum and
construct the final H$\beta$ maps by integrating the continuum-subtracted
spectra between 4856–4876 Å. The final maps are suitable for H ii region
identification using HIIphot.

###  4.2 H II Region Identification

HIIphot was designed to identify H ii regions and complexes (unresolved or
blended H ii regions) while also minimizing the inclusion of surrounding
diffuse ionized gas (or DIG). HIIphot works by first defining “seeds” at the
location of peak emission in H$\beta$ (or H$\alpha$), then iteratively grows
each “seed” and terminates only when the gradient of the H$\beta$ (or
H$\alpha$) surface brightness distribution matches a termination value, in
mandatory units of emission measure (EM), set by the user. The gradient of the
surface brightness distribution is more robust method of stopping uncontrolled
growth at lower S/N compared to using only the average local background level.
For each galaxy we apply the the same termination gradient, $\Delta=$ 5 EM
pc-1 or 2.43 $\times\ \rm{erg}\ \rm{s}^{-1}\ \rm{arcsec}^{-2}\ \rm{pc}^{-1}$,
as the recent PHANGS-MUSE work by Santoro et al. (2022).

HIIphot uses the PSF to convolve the input H$\beta$ map to different spatial
scales to identify seeds. Using a constant PSF for all galaxies can
potentially miss valid regions or generate non-physical regions. The PSF of
the input H$\beta$ emission map is required by HIIphot. For each galaxy
mosaic, we used the average PSF from the its KCWI pointings (see Table 7 in
Appendix A) as the input for HIIphot. The resulting 2D mask returned by
HIIphot contains H ii regions with smooth and reasonable boundaries, judged by
the distinction between clearly separated HII regions, the minimization of
spurious small and pixelated regions or runaway growth. In total, HIIphot
identifies $\sim$ 688 H ii regions or complexes across all of the KCWI
mosaics. This number is smaller than the 2169 potential H ii regions
identified for the Nebular catalog (Kreckel et al., 2019; Santoro et al.,
2022; Groves et al., 2023), and 2124 potential H ii regions from Congiu et al.
(2023), inside the same KCWI footprints. This is largely due to the
differences in angular resolution between KCWI and MUSE, and the three-fold
decrease in strength of H$\beta$ emission relative to H$\alpha$. In Figure 20
of Appendix C we show histograms of the H$\beta$ luminosity and radii for KCWI
and Nebular catalog regions as well a comparisons of the spatial masks in
Figures 21–27 of Appendix C. Additionally, we also present in Table 9 of
Appendix C the number of regions detected per galaxy.

Figure 4: Region boundaries returned by HIIPhot. Pixels within the (red)
boundaries are identified as corresponding to a potential H ii region. The
boundaries for the remaining sample galaxies are shown in Appendix C. The R.A.
and Dec. offset are centered on the R.A. and Dec. coordinates 02hr 46.0m
25.53s, -00∘ 29′ 38.8″.

###  4.3 Generation of Integrated H II Region Spectra

The KCWI H ii region masks, produced by HIIphot, are used to isolate and sum
the spectra in pixels belonging to each H ii region, resulting in an
integrated H ii region spectrum. To produce a matching MUSE H ii region mask
we transformed the H ii regions coordinates/boundaries from the KCWI pixel
grid onto the MUSE pixel grid. These are then used to construct MUSE
integrated spectra for each H ii region. The KCWI and MUSE H ii region spectra
for a single H ii region, with the full set of auroral lines highlighted, is
shown in Figure 5.

Figure 5: A KCWI, shown in the top panel, and MUSE, shown in the bottom panel,
H ii region spectrum. Example spectrum for an H ii region in NGC 5068. The
full wavelength range afforded by combining both KCWI and MUSE capture the
full set of optical auroral lines: [S II]$\lambda\lambda 4068,4078$, [O
III]$\lambda 4363$, [N II]$\lambda 5756$, [S III]$\lambda 6312$ and [O
II]$\lambda\lambda 7320,7330$, which are identified with zoomed insets.

We also produce integrated variance spectra for each H ii region. The variance
spectra for both KCWI and MUSE H ii regions are constructed by propagating the
pipeline-produced variance datacubes for pixels contained within each H ii
region boundary. In the case of MUSE, the datacubes have undergone an
additional convolution process in order to generate mosaics with uniform PSF
which introduces a correlation between neighboring pixels (Emsellem et al.,
2022). Because of the additional convolution, we generated the MUSE variance
spectra, assuming fully correlated conditions, by adding the pixel variance
spectra linearly (Taylor, 1997).

We verify this choice by comparing the median standard deviation of the
propagated MUSE variance spectra, $\sigma_{\mathrm{propagated}}$, to the
median standard deviation of the H ii region spectrum,
$\sigma_{\mathrm{measured}}$, in the emission line free wavelength range
$5400-5450$ Å. We measured an average ratio between the propagated and
measured error of $\sigma_{\mathrm{measured}}/\sigma_{\mathrm{propagated}}=$
$1.1\pm 0.2$, implying that we are appropriately propagating the error. We
also perform a similar comparison for KCWI, and find agreement,
$\sigma_{\mathrm{measured}}/\sigma_{\mathrm{propagated}}=1.9\pm 0.5$, between
the measured and propagated error using uncorrelated pixel error propagation.
To generate the appropriate variance we propagated the error using
uncorrelated, for KCWI, and fully correlated error, for MUSE, propagation
methods.

###  4.4 H II Region Stellar and Emission Line Fitting

We modeled the stellar continuum and emission lines of the integrated H ii
region spectra using the general PPXF toolkit. Although the LZIFU
implementation of PPXF allowed for the streamlined, full datacube, fitting of
the KCWI H$\beta$ emission line map, the general PPXF toolkit offers more
flexibility in the input fitting parameters. For example, we can input a
wavelength dependent LSF function as well as fix the kinematics between
emission lines of doublets and lines with similar levels of ionization. We
followed Emsellem et al. (2022) and fit the emission lines simultaneously with
the stellar continuum. This particular fitting recipe was chosen to mirror the
philosophy of the Mapping Nearby Galaxies at APO Data Analysis Pipeline (MaNGA
DAP, Law et al., 2016; Emsellem et al., 2022) and is suggested to mitigate the
biases on emission line fluxes introduced from the masking of stellar
absorption features around affected lines (Sarzi et al., 2006; Oh et al.,
2011; Belfiore et al., 2019). PPXF robustly fits the stellar continuum by
matching a set of templates to the observed H ii region spectrum. The
templates originate from the E-MILES library of simple stellar population
(SSP) models (Vazdekis et al., 2016). The SSP ages were between $0.15-13.5$
Gyr. Each age bin contained SSPs with the following metallicities,
$\mathrm{[Z/H]}=[-1.49,-0.96,-0.35,0.06,0.26,0.4]$. Typically, PPXF convolves
the SSP templates with a Gaussian model accounting for the spectral resolution
of the input spectrum and the stellar velocity dispersion. However, because
the KCWI line profile deviates significantly from a Gaussian (see Appendix B)
we convolved the PPXF templates with a 4-moment Gauss-Hermite function while
fitting KCWI emission lines. We constrained the fits of $h_{4}$ and $h_{3}$ in
the Gauss-Hermite functions to values listed in Appendix B. We performed the
PPXF fitting of the KCWI H ii region spectra independently from the MUSE H ii
region spectra. The fractional difference between the H ii region H$\beta$
flux for KCWI and MUSE is $-2.3\pm 7.5\%$.

We obtained errors on the emission line fluxes from the output of PPXF. The
output errors are considered reliable if the PPXF derived reduced
$\chi^{2}\approx 1$. Together the resulting fits for both the KCWI and MUSE
have an average reduced $\chi_{\mathrm{reduced}}^{2}\approx 2.0$, indicative
the input variance spectra are under-estimated. We obtained a better estimate
of the errors by re-scaling the returned errors, for each fit to KCWI and MUSE
spectra, by a factor of $\sqrt{\chi_{\mathrm{reduced}}^{2}}$ (Cappellari &
Emsellem, 2004; Cappellari, 2017; Emsellem et al., 2022).

###  4.5 Dust Correction

We derive the V-band extinction ($A_{V}$) for each H ii region using the
Balmer decrement. To evaluate this decrement while also taking into account
the errors on the measured H$\alpha$ and H$\beta$ emission, we construct a
distribution of the MUSE H$\alpha$/H$\beta$ ratios by sampling the error for
each line. Next, using PyNeb, we calculate the extinction by comparing the
average of the H$\alpha$/H$\beta$ distribution to the theoretical value
assuming the Case B recombination conditions $n_{e}=$103 cm-3 and
$T_{e}=10^{4}$ K (Storey & Hummer, 1995). We explore how changing the assumed
$T_{e}$ could affect the derived $A_{V}$ by sampling a range of temperatures
between 5000 K and 1.5$\times 10^{4}$ K and find that the standard deviation
of $A_{V}$ for a fixed Balmer decrement is $\sim 0.06$ mag. We apply the
wavelength dependent extinction correction assuming a O’Donnell (1994)
extinction curve. We present a histogram of the derived E($B$-$V$) in Figure
6. The average E($B$-$V$) for the regions is $0.30$ mag and corresponds to an
$A_{V}\sim 0.9$ mag. We find a negligible difference when using KCWI H$\beta$
in place of the MUSE H$\beta$ flux. We add that a recent investigation shows
that correcting, or not correcting, the Balmer lines for DIG contamination can
impact the measured E($B$-$V$) (Congiu et al., 2023). For the range of
E($B$-$V$) observed, a DIG-corrected E($B$-$V$) may return values of $A_{V}$
0.05–0.1 mag lower than presented here. We discuss the effects that
under/over-estimated extinction may have on the measured electron temperature
in Section 7.1.1.

Figure 6: Histogram of derived extinction, E(B-V). We measured E($B-V$) using
the O’Donnell (1994) extinction law for the regions identified by HIIPhot. A
(black dashed) line is located at the mean E($B-V$), 0.30 $\pm 0.14$ mag.

###  4.6 Diffuse Ionized Gas

Emission from Balmer transitions, [S II]$\lambda\lambda 6716,6731$, [N
II]$\lambda 6584$, and other lines originating from DIG surrounding the H ii
regions can contaminate the H ii region emission line fluxes of the same
transitions. Measuring the DIG contribution to H ii region line flux is only
beginning to being explored by IFU studies (see Belfiore et al., 2022; Congiu
et al., 2023). In order to removes regions with large degree of DIG
contamination we calculate the contrast between H ii/DIG emission following
the scheme outlined in Kreckel et al. (2022). First, we mask each H ii region
in the MUSE emission line maps for: H$\alpha$, H$\beta$, [S II], [N II] and [O
I]. For each emission line, we place a $10\arcsec\times 10\arcsec$ aperture
around each H ii region and measure the median DIG flux in pixels with
$\rm{S/N}>3$ and H$\alpha$ surface brightness
$\rm{Log}_{10}(SB_{\rm{H}\alpha}/[\mathrm{erg}\ \mathrm{s}^{-1}\
\mathrm{kpc}^{-2}])<38$ (see Belfiore et al., 2022). Finally, we calculate the
“integrated” DIG emission by multiplying the median DIG flux by the H ii
region size. For each DIG emission line, we calculate the percent contrast
between the measured H ii region flux and DIG flux. If the contrast between
the H ii region flux and DIG flux of any low-ionization emission line is
$<50$%, we exclude it from the sample.

###  4.7 Quality Assessment and Classification of Regions

The regions identified by HIIphot are potentially a mix of H ii regions,
planetary nebulae, and supernova remnants. Some also have low S/N. We perform
a set of cuts to reject non-H ii regions and/or low S/N spectra from the
catalog of H ii regions.

  * •

We exclude any H ii region whose centroid coordinates are within 2″ from the
edge of the mosaic. This step removes 86 H ii regions

  * •

We require the strong lines used for temperature determinations: H$\beta$,
H$\alpha$, [O III]$\lambda\lambda 4959,5007$ [O II]$\lambda 3727$, [N
II]$\lambda\lambda 6548,6584$, [S III]$\lambda 9069$ and [S II]$\lambda\lambda
6716,6731$ to be detected above a threshold of S/N $>$ 5\. This step cuts 72 H
ii regions from the sample.

  * •

Using the lines of H$\beta$, H$\alpha$, [O III] and [N II] we construct a
Baldwin-Phillips-Terlevich (BPT, Baldwin et al., 1981) diagram. We require H
ii regions to be consistent with photoionization by massive stars. Therefore,
we require them to fall below the empirical [O III]/H$\beta$ vs. [N
II]/H$\alpha$ (Kauffmann et al., 2003) and [O III]/H$\beta$ vs. [S
II]/H$\alpha$ (Kewley et al., 2001) lines. The BPT diagram showing the
location of each H ii region is shown in Figure 7. Out of the sample, 11 H ii
regions are above the empirical and theoretical line cut-offs and are removed.

  * •

We exclude H ii regions that fail our DIG contrast check in Section 4.6. This
step removes 98 H ii regions.

The constraints together remove 267 out of the 688 detected H ii regions
leaving 421 H ii regions remaining for use in our analysis. For the final
sample, comparison of line ratios in Figure 7 to model classifications in
Figure 3 of Congiu et al. (2023) suggest that final sample is consistent with
their H ii region classification.

Figure 7: BPT diagrams showing [O III]/H$\beta$ vs. [N II]/H$\alpha$, in the
left panel, and [O III]/H$\beta$ vs. [S II]/H$\alpha$ ratios, in the right
panel, for each region identified by HIIphot. Regions with emission line
ratios consistent with photoionization by stars are expected to populate the
parameter space below the theoretical (solid-black, Kewley et al., 2001) and
empirical (dashed-black, Kauffmann et al., 2003) classification lines. The H
ii regions are marked by filled markers. The regions above either of the
classification lines are labeled with non-filled markers and are rejected from
the catalog. Regions marked by ‘+’ symbols have been rejected by our other
considered constraints.

##  5 Measurement of Auroral Line Emission

In order to robustly measure the flux and uncertainty for the faint,
temperature sensitive auroral lines from [N II]$\lambda 5756$, [O
II]$\lambda\lambda 7320,7330$, [S II]$\lambda\lambda 4069,4076$, [O
III]$\lambda 4363$, and [S III]$\lambda 6312$ we implement a specific auroral
line fitting scheme in place of PPXF. This is necessary because any under/over
subtraction of the continuum at the location of the faint auroral lines can
bias the measured auroral line flux.

The framework of the auroral line fit is as follows. We first subtract the
stellar continuum spectrum fitted by PPXF. This results in an H ii region
spectrum that contains only emission lines and residuals from the continuum
subtraction. Next, we measure the standard deviation of the residuals,
$\sigma_{cont}$, in a region near the auroral line but also free of emission.
We then perform a large number of fits to the auroral line with random noise
added to each wavelength bin, drawn from a normal distribution with standard
deviation $\sigma_{cont}$. In each trial, we fit a single Gaussian (or double
depending on the auroral line) plus constant continuum level to the the
spectrum. The constant term is needed to correct for any remaining continuum
present in the continuum-subtracted spectrum. After completing the $N$ trials,
we calculate the average integrated line flux, $F_{avg}$, and the standard
deviation of the measured fluxes, $\sigma_{avg}$. If S/N $>3$ we consider the
auroral line detected. In Figure 29 of Appendix D we show the model fits and
residuals for a region in NGC 5068.

Although this general process is performed for all the auroral lines, the
auroral line from [OIII] is subject to additional constraints because emission
from the [O III]$\lambda 4363$ can be blended with [Fe II]$\lambda 4360$. This
has been observed in both stacked galaxy and individual H ii region spectra
(Curti et al., 2017; Berg et al., 2020; Arellano-Córdova & Rodríguez, 2020).
The strength of the [Fe II]$\lambda 4360$ has been observed to increase with
the metallicity of the gas (Curti et al., 2017), although continuum pumping
fluorescence contributes strongly to its emissivity (Rodríguez, 1999). [Fe
II]$\lambda 4360$ and [Fe II]$\lambda 4288$ arise from the same atomic upper
level and their relative intensities are independent of the physical
conditions of the gas and depend solely on the atomic transition
probabilities. Since $I(\lambda 4360)/I(\lambda 4288)=0.73$ (Mendoza et al.,
2023), if [Fe II]$\lambda 4360$ is detected, [Fe II]$\lambda 4288$ should be
present too.

To remove blending of the [O III]$\lambda 4363$ line by [Fe II]$\lambda 4360$,
we use the fixed ratio, $I(\lambda 4360)/I(\lambda 4288)=0.73$, to estimate
the degree of contamination by measuring the strength of the brighter [Fe
II]$\lambda 4288$ emission. We first fit the [Fe II]$\lambda 4288$ line using
a single Gaussian plus a constant continuum level model. The initial guesses
for the line center and width are taken from the parameters of H$\gamma$
returned by the PPXF fits. In the case that the [Fe II]$\lambda 4288$ is not
detected above a S/N threshold of 3, we instead use the 3$\sigma$ upper limit.
We next fit [O III]$\lambda 4363$ using a single Gaussian plus constant offset
model. The initial guesses for the kinematics of [O III]$\lambda 4363$ are
taken from the PPXF fit of [O III]$\lambda 4958$. After adding random noise,
we generate a model for the [Fe II]$\lambda$4360 using the best fit parameters
and errors derived from the fit to [Fe II]$\lambda 4288$, using the fixed
$I(\lambda 4360)/I(\lambda 4288)$ ratio. After generating the [Fe
II]$\lambda$4360 model, we subtract it from the trial spectrum and proceed to
then fit for the emission line flux originating from [O III]$\lambda 4363$. We
detect emission from [Fe II]$\lambda 4288$ in 30 H ii regions, two of these
are in regions with measurable [O III]$\lambda 4363$. From these, we
determined that the 3$\sigma$ upper-limit on the [Fe II]$\lambda 4360$ flux is
3.5 $\times 10^{-17}\rm{erg}\ \rm{s}^{-1}\ \rm{cm}^{-2}\ \rm{pc}^{-1}$ which
is $\sim 20\%$ of the average [O III]$\lambda 4363$ flux. We show in Figure 30
Appendix D an example [O III]$\lambda 4363$ fit which has had significant [Fe
II]$\lambda 4360$ contribution removed.

##  6 H II Region Nebular, Environmental, and Stellar Properties

We assess the ionized gas physical conditions—$n_{\rm{e}}$, $T_{e}$, and
ionization parameter $U$—of each H ii region using a subset of the dust-
corrected emission line fluxes. We also measure a number of characteristics of
the H ii region’s environment and local stellar population, as described
below.

###  6.1 Electron Density

Using PyNeb we calculate the electron density $n_{\rm{e}}$ for each H ii
region using the [S II]$\lambda\lambda 6716,6731$ doublet. Based on our
constraints discussed in Section 4.7, each H ii region is guaranteed to have
measured emission from this doublet at S/N $>5$. The [O II]$\lambda\lambda
3726,3729$ doublet is also commonly used to estimate $n_{\rm{e}}$. The atomic
levels responsible for [S II]$\lambda\lambda 6716,6731$ and [O
II]$\lambda\lambda 3726,3729$ both have critical densities, when collisional
and radiative de-excitation are occurring at equal rates, that are of order
$10^{3}$ cm-3. The critical densities, as well as other references for the
atomic data used, are listed in Table 3. Both the [S II]$\lambda\lambda
6716,6731$ and [O II]$\lambda\lambda 3726,3729$ doublets are sensitive to
densities $10^{2}\ \mathrm{cm}^{-3}<n_{\rm{e}}<10^{3.5}\ \mathrm{cm}^{-3}$.
However, in the KCWI measurements, the [O II]$\lambda\lambda 3726,3729$
doublet is unresolved, and therefore we do not use it to measure $n_{\rm{e}}$.
With the exception of a handful of regions, the measured electron densities
are in the low-density limit, $n_{\rm{e}}\ <\ 100$ $\mathrm{cm}^{-3}$.

Table 3: Transitions Probabilities, Collision Strengths, and Critical Densities for the relevant emission lines used in the PyNeb Temperature Determinations. Ion | Transition Probabilities | Collision Strengths | $n_{\mathrm{crit,nebular}}$
---|---|---|---
|  |  | ($10^{3}$ cm-3)
$[$O II$]$ | Zeippen (1982) | Kisielius et al. (2009) | 2aaAverage critical density for [O II]$\lambda\lambda 3726,3729$ .
$[$O III$]$ | Froese Fischer & Tachiev (2004) | Storey et al. (2014) | 691
$[$N II$]$ | Froese Fischer & Tachiev (2004) | Tayal (2011) | 88
$[$S II$]$ | Rynkun et al. (2019) | Tayal & Zatsarinny (2010) | 3bbAverage critical density for [S II]$\lambda\lambda 6716,6731$ .
$[$S III$]$ | Froese Fischer et al. (2006) | Tayal & Gupta (1999) | 543ccAverage critical density for [S III]$\lambda\lambda 9069,9532$ .

###  6.2 H II Region Electron Temperatures

The auroral lines allow the determination of the electron temperatures for the
O+, O2+, N+, S+ and S2+ ions, or $T_{\rm{e},\text{[{O}{II}]}}$,
$T_{\rm{e},\text{[{O}{III}]}}$, $T_{\rm{e},\text{[{N}{II}]}}$,
$T_{\rm{e},\text{[{S}{II}]}}$, and $T_{\rm{e},\text{[{S}{III}]}}$,
respectively. The temperatures are calculated via PyNeb using the collision
strengths and transition probabilities, listed in Table 3, as well the
measured upper-limits of the electron density, $n_{\rm{e}}$, to convert an
auroral-to-nebular ratio to temperature. The critical densities of the
requisite lines, see Table 3, used to estimate $T_{\rm{e},\text{[{O}{III}]}}$,
$T_{\rm{e},\text{[{N}{II}]}}$, and $T_{\rm{e},\text{[{S}{III}]}}$ are high
enough such that the auroral-to-nebular lines ratios are insensitive to choice
of $n_{\rm{e}}<10^{4}$ $\mathrm{cm}^{-3}$. The density sensitivity in the
lines used for $T_{\rm{e},\text{[{O}{II}]}}$ and $T_{\rm{e},\text{[{S}{II}]}}$
begins at smaller densities $n_{\rm{e}}\approx 10^{3}$ $\mathrm{cm}^{-3}$ and
is further discussed in Section 7.2. The uncertainty for each temperature
measurement is the standard deviation of the distribution of temperatures
constructed by Monte Carlo sampling of the error for each auroral and nebular
line included in the temperature determination. We summarize the number of
detections and median temperature for each ion in Table 4.

Table 4: Median temperature for each ion. Te,Ion | MedianaaThe 50th $\pm$ the 16th – 84th percentile of the measured temperatures. |  $N_{\mathrm{regions}}$bbNumber of H ii regions with temperature measured from the particular ion.
---|---|---
| [104 K] |
$T_{\rm{e},\text{[{O}{II}]}}$ | 0.95${}^{+0.22}_{-0.12}$ | 156
$T_{\rm{e},\text{[{N}{II}]}}$ | 0.81${}^{+0.12}_{-0.09}$ | 245
$T_{\rm{e},\text{[{S}{II}]}}$ | 0.94${}^{+0.19}_{-0.14}$ | 305
$T_{\rm{e},\text{[{S}{III}]}}$ | 0.89${}^{+0.30}_{-0.10}$ | 143
$T_{\rm{e},\text{[{O}{III}]}}$ | 1.20${}^{+0.35}_{-0.16}$ | 26

###  6.3 H II Region Ionization Parameters

The ionization parameter is an indicator of the strength of the ionizing
radiation field. The Strömgren sphere descriptions of H ii regions define the
ionization parameter as $U=Q_{0}/(4\pi R^{2}n_{\mathrm{H}}c)$, where $Q_{0}$
is the emission rate of photons capable of ionizing hydrogen (i.e., with
energy $>13.6$ eV), $R$ is the radius of the ionized region, $n_{\mathrm{H}}$
and $c$ are the hydrogen density and speed of light. However, calculating the
ionization parameter using this definition is difficult (Kreckel et al., 2019,
2022) as resolved studies show that H ii regions exhibit a range of non-
spherical morphologies and non-uniform densities (Wood & Churchwell, 1989).

Instead, we trace the ionization parameter using the [S III]/[S II] and [O
III]/[O II] emission line ratios. Photoionization modeling has shown that both
the [S III]$\lambda\lambda 9069,9532$/[S II]$\lambda\lambda 6716,6731$ and [O
III]$\lambda\lambda 4959,5007$/[O II]$\lambda\lambda 3626,3729$ correlate with
the ionization parameter, $U$ (Dors et al., 2011). Although it has a positive
correlation with the ionization parameter, the [O III]/[O II] also has a
secondary dependence on metallicity, increasing with decreasing metallicity.
[S III]/[S II] is not as sensitive to metallicity, making the ratio a more
reliable tracer of the ionization parameter (Kewley & Dopita, 2002). While [S
III]$\lambda 9532$ is not observed with MUSE, we measure the [S III]$\lambda
9069$ and assume the fixed theoretical line ratio of [S III]$\lambda
9532/\lambda 9069=2.5$ (Froese Fischer et al., 2006) in all calculations.

###  6.4 ALMA-CO: Intensity, Peak Temperature, and Velocity Dispersion

Using the ALMA 12m+7m+TP datacubes, we calculate moment 0, and 2 (integrated
intensity and velocity dispersion) for molecular gas near each H ii region.
Because molecular and ionized gas are not entirely co-spatial at our
resolution, it is necessary to make a selection to capture gas near the H ii
region. One possibility is to match H ii regions to molecular clouds via a
nearest-neighbor algorithm (Grasha et al., 2019; Zhang et al., 2021;
Zakardjian et al., 2023) We instead choose to integrate the ALMA spectrum
contained in the H ii region boundaries in order to measure the properties of
the molecular gas closest in projection to the the ionized gas (i.e. in front,
behind, or blended due to the resolution of KCWI) and likely affected by the
radiative feedback. To measure the CO spectra, we reproject the H ii region
masks onto the grid of the ALMA datacubes and integrate the ALMA spectra for
pixels located inside the footprint of each H ii region.

From the ALMA spectrum, we then calculate the integrated intensity, $I_{\rm
CO}$, peak temperature, $T_{\rm peak}$ and the velocity dispersion,
$\sigma_{v,\mathrm{CO}}$. In order to accurately measure these CO moments in
the presence of noise, we construct a signal mask following the basic approach
from Leroy et al. (2021b). To do this, we locate the velocity channel which
contains the peak emission and construct an integration window around this
channel by including velocity channels with signal above the $1\sigma$ noise.

As a check on our analysis of molecular gas, we compare the calculated
velocity dispersions to those from a sample of nearest-neighbor matched H ii
regions-GMC’s constructed by Zakardjian et al. (2023). We find that our
$\sigma_{v,{\rm CO}}$ span a similar range, up to Log10($\sigma_{v,{\rm
CO}}/$[km s-1])=1.5, with an average and standard deviations of
Log${}_{10}(\sigma_{v,{\rm CO}}/$[km s-1])=0.88 $\pm$ 0.25; in line with the
average from Zakardjian et al. (2023) which suggests that this method of
extracting and measuring the CO properties is reasonable.

###  6.5 H ii Region Compact Clusters and Associations

In order to test for correlations between the young stellar populations that
power H ii regions and their electron temperatures, we compile the stellar
mass and age of compact clusters and multi-scale stellar associations within
our KCWI H ii regions using results from HST observations.

We match the HST clusters to individual H ii regions by simply selecting all
the clusters whose on-sky coordinates fall inside any H ii region’s spatial
footprint. For the associations, we match the NUV selected, 32 pc scale,
stellar associations to the individual H ii regions in the same manner as
Scheuermann et al. (2023). The associations catalog comes with spatial masks
identifying the footprint of all the detected associations. Because the
association masks have a finer pixel scale than the KCWI H ii region masks, we
reproject the KCWI mask onto the association mask pixel grid. We find for
cluster matches that only 99 or 23.5% of our H ii regions are matched to a
single cluster. 274 or 65% of our H ii regions have zero matches. For
association matches, only 149 or 35.4% of our H ii regions are matched to a
single association 117 or 27.8% of our H ii regions have zero matches.

For the remaining regions with more than one clusters or association match,
because we expect that the youngest and most massive clusters or associations
contribute most to the overall ionization of the H ii region, we assign to
each H ii region the age (mass) of the youngest (most massive) available
cluster/association.

##  7 Results

We present electron temperatures derived using the auroral-to-nebular line
ratios from [O II], [N II], [S II], [S III], and [O III]. We construct
$T_{\rm{e}}$–$T_{\rm{e}}$ diagrams to compare any multi-ionization zone
$T_{\rm{e}}$ relationships to recently measured and/or modeled trends. We then
compare the temperatures to properties of the ionized gas such as electron
density, $n_{\rm{e}}$, and ionization parameter, $U$. We also relate the
temperatures to properties of the molecular gas and stellar populations.

###  7.1 Temperature-Temperature Relations

####  7.1.1 Low-Ionization Zone

We show $T_{\rm{e}}$–$T_{\rm{e}}$ relations for our sample H ii regions in
Figures 8, 9 and 10. The left panel displays the individual H ii region
measurements. To better see $T_{\rm{e}}$–$T_{\rm{e}}$ relationships in
individual galaxies we show in the right panel the $T_{\rm{e}}$–$T_{\rm{e}}$
relations for temperatures binned in steps of 2000 K in the x-axis, with
minimum 2 H ii regions per bin, for each galaxy. Alongside the binned data we
show recent $T_{\rm{e}}$–$T_{\rm{e}}$ relations from Berg et al. (2020,
hereafter CHAOS), Zurita et al. (2021, hereafter Z21), Garnett (1992,
hereafter G92), Vale Asari et al. (2016, hereafter BOND), and Méndez-Delgado
et al. (2023b, hereafter MD23).

The low-ionization zone $T_{\rm{e}}$–$T_{\rm{e}}$ relations are shown in
Figure 8. The top figure shows the
$T_{\rm{e},\text{[{N}{II}]}}$–$T_{\rm{e},\text{[{O}{II}]}}$ comparison. We
observe that $T_{\rm{e},\text{[{O}{II}]}}$ gives higher values than
$T_{\rm{e},\text{[{N}{II}]}}$ by 1000 K on average. The magnitude of this
offset is largest when $T_{\rm{e},\text{[{O}{II}]}}$$>1.0\times 10^{4}$ K. The
$T_{\rm{e},\text{[{O}{II}]}}$$>$$T_{\rm{e},\text{[{N}{II}]}}$ inequality is
also reflected in the relations from Z21, CHAOS, and MD23. With the exception
of the CHAOS trend, the slope of the relations also indicate a steeper rise in
$T_{\rm{e},\text{[{O}{II}]}}$ than $T_{\rm{e},\text{[{N}{II}]}}$. Although the
data lie below the 1–1 line, they also do not clearly prefer either CHAOS,
MD23, and Z21 relations.

The comparison between
$T_{\rm{e},\text{[{S}{II}]}}$–$T_{\rm{e},\text{[{N}{II}]}}$ is shown in the
middle panel of Figure 8. We observe that there is an offset in between
$T_{\rm{e},\text{[{S}{II}]}}$–$T_{\rm{e},\text{[{N}{II}]}}$ in favor of hotter
$T_{\rm{e},\text{[{S}{II}]}}$. Similar to
$T_{\rm{e},\text{[{O}{II}]}}$–$T_{\rm{e},\text{[{N}{II}]}}$, the scatter in
the temperatures is larger in $T_{\rm{e},\text{[{S}{II}]}}$ rather than
$T_{\rm{e},\text{[{N}{II}]}}$. A majority of the binned data points congregate
above the literature relations. Only a few exceptions agree with the 1–1 line.

We show the temperatures for 129 regions with both
$T_{\rm{e},\text{[{S}{II}]}}$-$T_{\rm{e},\text{[{O}{II}]}}$ in the bottom-left
panel. We observe large scatter for both $T_{\rm{e},\text{[{S}{II}]}}$ and
$T_{\rm{e},\text{[{O}{II}]}}$. Again, the binned data do not strongly prefer
any of the selected relations.

We observe in the $T_{\rm{e}}$–$T_{\rm{e}}$ relations that
$T_{\rm{e},\text{[{O}{II}]}}$ and $T_{\rm{e},\text{[{S}{II}]}}$ are
systematically hotter than $T_{\rm{e},\text{[{N}{II}]}}$. Based on
photoionization models (e.g., Campbell et al., 1986; Garnett, 1992), the low-
ionization zone temperatures are expected to be equal. There are systematic
effects that could increase $T_{\rm{e},\text{[{O}{II}]}}$ and
$T_{\rm{e},\text{[{S}{II}]}}$ temperatures such that
$T_{\rm{e},\text{[{O}{II}]}}$$\sim$$T_{\rm{e},\text{[{S}{II}]}}$$>$$T_{\rm{e},\text{[{N}{II}]}}$.
The wide wavelength range between [O II]$\lambda\lambda 7320,7330$ and [O
II]$\lambda\lambda 3726,3729$ as well as [S II]$\lambda\lambda 4068,4076$ and
[S II]$\lambda\lambda 6716,6731$, make these ratios sensitive to the applied
reddening correction. An overestimate in the extinction would lead to an
underestimate of $T_{\rm{e},\text{[{O}{II}]}}$ and underestimate of
$T_{\rm{e},\text{[{S}{II}]}}$. At the same time this would leave
$T_{\rm{e},\text{[{N}{II}]}}$ relatively unchanged due to the closeness in
wavelength of the requisite lines. We see that both
$T_{\rm{e},\text{[{O}{II}]}}$ and $T_{\rm{e},\text{[{S}{II}]}}$ are greater
than $T_{\rm{e},\text{[{N}{II}]}}$ which is not compatible with the effects of
overestimated extinction correction. Despite the above, our temperatures
hierarchy could potentially be produced by under-estimated extinction as this
would lead to over-estimates of $T_{\rm{e},\text{[{O}{II}]}}$ and
$T_{\rm{e},\text{[{S}{II}]}}$ while again leaving
$T_{\rm{e},\text{[{N}{II}]}}$ unchanged. We tested two different extinction
prescriptions (O’Donnell, 1994; Fitzpatrick, 1999) and found no change in our
results. Due to the possible DIG contribution to the measured extinction, it
is more likely that we could be overestimating the extinction than
underestimating it. Telluric contamination to the line emission at [O
II]$\lambda\lambda 7320,7330$ could be an additional systematic error, but is
unlikely to be significant in our data because the [O II]$\lambda
7320$/$\lambda 7330$ ratio measured for our sample is $1.27\pm 0.3$, in
agreement with values predicted by the transition probabilities and the
collisional strengths (Zeippen, 1982; Kisielius et al., 2009) and observed in
nearby H ii regions (Seaton & Osterbrock, 1957; Kaler et al., 1976; Yates et
al., 2020; Méndez-Delgado et al., 2022). The [S II]$\lambda\lambda 4068,4076$
doublet may contain $\sim 10\%$ contamination due to O II recombination
emission near the location of [S II]$\lambda\lambda 4068,4076$ (Méndez-Delgado
et al., 2023b) that could bias $T_{\rm{e},\text{[{S}{II}]}}$ upwards by $\sim
0.2\times 10^{4}$ K for $T_{\rm{e},\text{[{S}{II}]}}$ $<2.0\times 10^{4}$ K.
The KCWI spectral resolution is too low to separate out any contamination in
the measurements of [S II]$\lambda\lambda 4068,4076$. Lastly,
$T_{\rm{e},\text{[{O}{II}]}}$ and $T_{\rm{e},\text{[{S}{II}]}}$ have a higher
sensitivity to electron density inhomogeneities than
$T_{\rm{e},\text{[{N}{II}]}}$ (Rubin, 1989; Osterbrock & Ferland, 2006;
Méndez-Delgado et al., 2023b). [O II] and [S II] start to be dependent on
density around $n_{\rm{e}}\approx 10^{3}$ cm-3 while [N II] is independent of
density up to $n_{\rm{e}}\approx 10^{4}$ cm-3. We explore this dependence
further in Section 7.2.

####  7.1.2 Intermediate Ionization Zone

The $T_{\rm{e}}$–$T_{\rm{e}}$ relations between the low and intermediate-
ionization zone temperatures are shown in Figure 9. Given the lower ionization
potentials of O+, N+, and S+ with respect to S++, it is possible to have
differences in the temperature in the different ionization zones (Garnett,
1992). We compare our observed $T_{\rm{e}}$–$T_{\rm{e}}$ relations to
predictions from photoionization models for giant H ii regions produced by the
Bayesian Oxygen and Nitrogen abundance project BOND and found in the Mexican
Million Models database (Morisset, 2009). When available, we also compare the
observations to relations from CHAOS, MD23, and Z21.

The top panel of Figure 9 shows the temperature comparisons between
$T_{\rm{e},\text{[{S}{III}]}}$ and $T_{\rm{e},\text{[{O}{II}]}}$. The (BOND)
models predict that the temperature $T_{\rm{e},\text{[{S}{III}]}}$ should be
greater than $T_{\rm{e},\text{[{O}{II}]}}$ by a small constant offset across
the full temperature range. However, our data suggest
$T_{\rm{e},\text{[{O}{II}]}}$ rises faster than
$T_{\rm{e},\text{[{S}{III}]}}$. Our binned data favor the empirical trend line
from Z21. However, there is a grouping of points from NGC 5068 that show
contrasting behavior in both the individual and binned data comparisons.

The $T_{\rm{e},\text{[{S}{III}]}}$ and $T_{\rm{e},\text{[{S}{II}]}}$
comparison, comprised of 108 H ii regions, is shown in the middle panel of
Figure 9. We observe in the middle-left panel scatter $\approx$ $0.2\times
10^{4}$ K with no clear trend for higher or lower temperatures on either axis.
The binned data do not reveal a preference for any trend.

The bottom panel of Figure 9 shows the $T_{\rm{e},\text{[{S}{III}]}}$ and
$T_{\rm{e},\text{[{N}{II}]}}$ comparison. $T_{\rm{e},\text{[{S}{III}]}}$ and
$T_{\rm{e},\text{[{N}{II}]}}$ exhibit the smallest scatter, $\sim 0.1\times
10^{4}$ K, of the $T_{\rm{e}}$–$T_{\rm{e}}$ relations presented here. In the
left panel, we see that with increasing $T_{\rm{e},\text{[{N}{II}]}}$, the
trend of the individual points slightly favor hotter
$T_{\rm{e},\text{[{S}{III}]}}$. For the $T_{\rm{e},\text{[{S}{III}]}}$ and
$T_{\rm{e},\text{[{N}{II}]}}$ comparison there are available empirical
relations from CHAOS, Z21, and MD23 which we overlay in addition to the BOND
models. We observe a large fraction of the binned data that lie near the CHAOS
model, and Z21 relations at low $T_{\rm{e},\text{[{N}{II}]}}$. Regions with
hotter $T_{\rm{e},\text{[{N}{II}]}}$ and $T_{\rm{e},\text{[{S}{III}]}}$ would
be needed to to further differentiate between the models and empirical trends.

####  7.1.3 High Ionization Zone

We show $T_{\rm{e}}$–$T_{\rm{e}}$ relations between [O III] temperatures
(which trace the high ionization zone) and those from the low/intermediate
ionization zones in Figure 10. Additionally, we overlay
$T_{\rm{e}}$–$T_{\rm{e}}$ relations from CHAOS, Z21, and G92. Although we
observe some H ii regions with $T_{\rm{e}}$–$T_{\rm{e}}$ relations that agree
with literature relations, we also observe numerous H ii regions with much
higher $T_{\rm{e},\text{[{O}{III}]}}$ values of than what is predicted by
models with the given low and intermediate ionization zone temperatures.

The large excess in $T_{\rm{e},\text{[{O}{III}]}}$ for regions with cooler low
and intermediate ionization zone temperatures is unexpected. At the $\sim$
Solar metallicities of our sample, the relative flux [O III]$\lambda 4363$ is
expected to be $<10^{-2}\times$ H$\beta$ (Berg et al., 2015). The small number
of [O III]$\lambda 4363$ detections reflects this. Because we do not expect to
detect the line in most H ii regions, the ones we do detect may be unusual
cases or statistical outliers, especially given that the average S/N of the [O
III] detections is $\sim$4\. Despite this, there has been evidence (Peimbert
et al., 1991; Binette et al., 2012) that shock excitation can preferentially
enhance the high ionization zone temperature. This scenario is explored in
Section 7.2.5.

Figure 8: $T_{e}$–$T_{e}$ relationships for the low-ionization zone
temperatures measured from H ii regions in nearby galaxies. The left panel in
each row displays the individual temperatures and errors with each point
colored according to the host galaxy. The right panel in each row shows the
$T_{e}$–$T_{e}$ relations for temperatures binned in steps of 2000 K in the
x-axis, with minimum 2 H ii regions per bin, for each galaxy, compared to
$T_{e}$-$T_{e}$ trend lines from Z21 (blue dash), CHAOS (purple dot-dash), and
MD23 (green dot-dash). To aid the eye, we include in (black-dash) the 1–1
line.

Figure 9: $T_{\rm{e}}$–$T_{\rm{e}}$ relationships between the low and
intermediate ionization zone temperatures measured from H ii regions in nearby
galaxies. The temperatures measured from H ii regions in this work are color
coded by host galaxy. The left panel in each row displays the individual
temperatures and errors with each point colored according to the host galaxy.
The right panel in each row show the $T_{\rm{e}}$–$T_{\rm{e}}$ relations for
temperatures binned in steps of 2000 K in the x-axis, with minimum 2 H ii
regions per bin, for each galaxy, compared to $T_{\rm{e}}$–$T_{\rm{e}}$ trend
lines from Z21 (blue dash), CHAOS (purple dot-dash), BOND (green dot-dash) and
MD23 (red dot-dash). To aid the eye, we include in (black dash) the 1–1 line
Figure 10: $T_{\rm{e}}$–$T_{\rm{e}}$ relationships for the high ionization
zone. The temperatures measured from H ii regions in this work are color coded
by host galaxy. The (blue-dashed) line show the $T_{\rm{e}}$–$T_{\rm{e}}$
relationships from Z21. The (red-dashed) line show trends from CHAOS. The
(purple-dashed) and (green-dashed) line trends from the photoionization models
of G92 and BOND. To aid the eye, we include show in black line the 1–1 line in
each panel.

###  7.2 Temperature Differences Compared to H II Region, Stellar Population,
and Molecular Gas Properties.

Studies have shown that H ii region temperatures for different ionization
zones can be differently impacted by the ionization state of the gas (Berg et
al., 2020; Arellano-Córdova & Rodríguez, 2020; Yates et al., 2020) and by
temperature/density fluctuations (Méndez-Delgado et al., 2023a, b). However,
the physical processes that drive these differences in temperature remain
poorly constrained and uncertain (Garnett et al., 1991; Nicholls et al.,
2020). The KCWI+MUSE H ii region observations allow us to investigate
potential sources driving temperature differences,
$\Delta$($T_{\mathrm{ion,1}}$,$T_{\mathrm{ion,2}}$)=$T_{\mathrm{ion,1}}$$-$$T_{\mathrm{ion,2}}$,
between the low, intermediate, and high-ionization zone temperatures. We
compare temperature differences with H ii region properties derived from
emission line diagnostics; with the properties of the surrounding molecular
gas measured from ALMA (Leroy et al., 2021b); and with stellar population
masses/ages from SED fitting to HST photometry (Lee et al., 2022; Thilker et
al., 2022; Larson et al., 2023).

We summarize the comparisons to ionized gas properties in Table 5 and the
molecular gas/stellar population properties in Table 6. To gauge the
significance and degree of correlation between each comparison we calculate
the Spearman Rank Correlation Coefficient, $\rho$, and associated p-value,
$p$. A correlation is judged to be significant if it exhibits $p$ $\lesssim
10^{-3}$. The strength of the correlation is separated into the following
regimes: $1>|\rho|>0.8$ corresponds to a strong correlation, $0.8>|\rho|>0.4$
corresponds to a moderate correlation and $0.4>|\rho|>0$ identifies a weak or
no correlation.

We do not report any significant correlations between $\Delta T_{\rm{e}}$ with
any of the following properties: $n_{\rm{e}}$, integrated CO intensity, CO
peak temperature, cluster mass, cluster age, and association age. Figures
showing the $\Delta T_{\rm{e}}$ comparisons to these parameters are shown in
Appendix E. However, we do measure correlations with ionization parameter,
association mass, and molecular gas velocity dispersion. We discuss these
correlations in the following subsections.

Table 5: Summary of the p-values and Spearman Rank coefficients for comparisons between $\Delta T_{\rm{e}}$ and properties of the H ii region ionized Gas. $\Delta T_{\rm{e}}$ | [S III]/[S II] | [O III]/[O II] | $n_{\rm{e}}$
---|---|---|---
| ($\rho$, p) | ($\rho$, p) | ($\rho$, p)
$\Delta$($T_{\rm{e},\text{[{O}{II}]}}$, $T_{\rm{e},\text{[{N}{II}]}}$) | (0.00, 0 | (0.02, 0.79) | (-0.05, 0 |
$\Delta$($T_{\rm{e},\text{[{N}{II}]}}$, $T_{\rm{e},\text{[{S}{II}]}}$) | (-0.34, $<10^{-3}$) | (-0.03, 0.66) | (-0.28, 0 |
$\Delta$($T_{\rm{e},\text{[{S}{II}]}}$, $T_{\rm{e},\text{[{O}{II}]}}$) | (0.34, $<10^{-3}$) | (0.07, 0.40) | (0.27, 0 |
$\Delta$($T_{\rm{e},\text{[{O}{II}]}}$, $T_{\rm{e},\text{[{S}{III}]}}$) | (0.06, 0 | (0.02, 0.82) | (0.11, 0 |
$\Delta$($T_{\rm{e},\text{[{N}{II}]}}$, $T_{\rm{e},\text{[{S}{III}]}}$) | (0.26, 0 | (0.02, 0.81) | (0.19, 0 |
$\Delta$($T_{\rm{e},\text{[{S}{II}]}}$, $T_{\rm{e},\text{[{S}{III}]}}$) | (0.54, $<10^{-3}$) | (0.03, 0.70) | (0.24, 0 |
$\Delta$($T_{\rm{e},\text{[{O}{II}]}}$, $T_{\rm{e},\text{[{O}{III}]}}$) | (0.62, $<10^{-3}$) | (0.85, $<10^{-3}$) | (0.38, 0 |
$\Delta$($T_{\rm{e},\text{[{N}{II}]}}$, $T_{\rm{e},\text{[{O}{III}]}}$) | (0.72, $<10^{-3}$) | (0.82, $<10^{-3}$) | (0.13, 0 |
$\Delta$($T_{\rm{e},\text{[{S}{II}]}}$, $T_{\rm{e},\text{[{O}{III}]}}$) | (0.78, $<10^{-3}$) | (0.78, $<10^{-3}$) | (0.15, 0 |
$\Delta$($T_{\rm{e},\text{[{S}{III}]}}$, $T_{\rm{e},\text{[{O}{III}]}}$) | (0.59, 0 | (0.80, $<10^{-3}$) | (0.00, 1 |
Table 6: A summary of the p-values and Spearman rank coefficients for the whole set of comparisons between $\Delta T_{e}$ and H ii region stellar population and molecular gas properties. $\Delta T_{e}$ | AgeCluster | MCluster | AgeAssoc | MAssoc | $I_{\mathrm{CO}}$ | $T_{\mathrm{peak}}$ | $\sigma_{v,\rm{CO}}$
---|---|---|---|---|---|---|---
| ($\rho$, p) | ($\rho$, p) | ($\rho$, p) | ($\rho$, p) | ($\rho\ $, p) | ($\rho$, p) | ($\rho$, p)
$\Delta$($T_{\rm{e},\text{[{O}{II}]}}$, $T_{\rm{e},\text{[{N}{II}]}}$) | (-0.25, 0 | (-0.22, 0.13) | (0.0, 0 | (-0.11, 0.35) | (0.08, 0 | (0.20, 0.07) | (-0.16, 0 |  |  |
$\Delta$($T_{\rm{e},\text{[{N}{II}]}}$, $T_{\rm{e},\text{[{S}{II}]}}$) | (-0.01, 0 | (0.03, 0.79) | (0.04, 0 | (0.12, 0.20) | (-0.13, 0 | (-0.16, 0.67) | (0.01, 0 |  |  |
$\Delta$($T_{\rm{e},\text{[{S}{II}]}}$, $T_{\rm{e},\text{[{O}{II}]}}$) | (0.00, 0 | (0.21, 0.15) | (-0.12, 0 | (-0.06, 0.61) | (0.11, 0 | (0.03, 0.82) | (0.11, 0 |  |  |
$\Delta$($T_{\rm{e},\text{[{O}{II}]}}$, $T_{\rm{e},\text{[{S}{III}]}}$) | (-0.04, 0 | (-0.15, 0.39) | (-0.02, 0 | (0.11, 0.43) | (0.10, 0 | (0.15, 0.28) | (0.14, 0 |  |  |
$\Delta$($T_{\rm{e},\text{[{N}{II}]}}$, $T_{\rm{e},\text{[{S}{III}]}}$) | (-0.08, 0 | (0.17, 0.28) | (-0.23, 0 | (0.37, 0.001) | (0.05, 0 | (-0.21, 0.62) | (0.39,$<10^{-3}$) |  |  |
$\Delta$($T_{\rm{e},\text{[{S}{II}]}}$, $T_{\rm{e},\text{[{S}{III}]}}$) | (-0.32, 0 | (0.18, 0.24) | (-0.21, 0 | (0.13, 0.25) | (0.21, 0 | (0.08, 0.52) | (0.24, 0 |  |  |
$\Delta$($T_{\rm{e},\text{[{O}{II}]}}$, $T_{\rm{e},\text{[{O}{III}]}}$) | (0.19, 0 | (-0.42, 0.09) | (0.65, 0 | (-0.51, 0.02) | (0.10, 0 | (0.37, 0.16) | (-0.51, 0 |  |  |
$\Delta$($T_{\rm{e},\text{[{N}{II}]}}$, $T_{\rm{e},\text{[{O}{III}]}}$) | (0.30, 0 | (-0.40, 0.12) | (0.45, 0 | (-0.33, 0.12) | (0.23, 0 | (0.46, 0.05) | (-0.45, 0 |  |  |
$\Delta$($T_{\rm{e},\text{[{S}{II}]}}$, $T_{\rm{e},\text{[{O}{III}]}}$) | (0.30, 0 | (-0.37, 0.15) | (0.32, 0 | (-0.47, 0.03) | (0.17, 0 | (0.42, 0.90) | (-0.44, 0 |  |  |
$\Delta$($T_{\rm{e},\text{[{S}{III}]}}$, $T_{\rm{e},\text{[{O}{III}]}}$) | (0.48, 0 | (-0.39, 0.15) | (0.42, 0 | (-0.42, 0.63) | (0.05, 0 | (0.38, 0.15) | (-0.69,$<10^{-3}$) |  |  |

####  7.2.1 The Impact of Density Inhomogeneities on [S II] and [O II]
Temperatures

As expected, given that a majority of the $n_{\rm{e}}$ are in the low-density
limit, we do not observe any correlations between any $\Delta T_{\rm{e}}$
versus $n_{\rm{e}}$. However, recent studies have suggested that the
temperatures obtained from the auroral-to-nebular lines ratios of [O II] and
[S II] can be biased upwards due to the presence of density inhomogeneities,
even while the average density is underestimated by nebular doublet line [O
II]$\lambda\lambda 3726,3729$ and [S II]$\lambda\lambda 6716,6731$
diagnostics.

For example, in a sample of 190 high-signal-noise spectra of H ii regions and
other photoionized nebulae, Méndez-Delgado et al. (2023b) observed
systematically hotter $T_{\rm{e},\text{[{O}{II}]}}$ and
$T_{\rm{e},\text{[{S}{II}]}}$ relative to $T_{\rm{e},\text{[{N}{II}]}}$,
similar to what we observe, which they attribute to the presence of density
inhomogeneities. The atomic levels responsible for the nebular lines [S
II]$\lambda\lambda 6731,6716$ and [O II]$\lambda\lambda 3726,3729$, listed in
Table 3, have critical densities of order $10^{3}$ cm-3 which are at least two
orders of magnitude lower than the critical densities for the nebular levels
of [N II]. In addition, the auroral levels of the same ions have very high
critical densities. This makes the auroral-to-nebular temperature diagnostics
of $T_{\rm{e},\text{[{O}{II}]}}$ and $T_{\rm{e},\text{[{S}{II}]}}$ density
sensitive above $10^{3}$ cm-3, and therefore susceptible to biases if there
are important contributions to the line flux from gas above that density. The
auroral-to-nebular ratio of [N II], however, are not susceptible to such
sensitivity until much higher densities.

Both [S II]$\lambda\lambda 6731,6716$ and [O II]$\lambda\lambda 3726,3729$
nebular line doublet ratios serve as density diagnostics for densities
$10^{2}$ cm${}^{-3}<$ $n_{\rm{e}}$ $<$ $10^{3.5}$ cm-3. Furthermore, because
of the bias described above, Méndez-Delgado et al. (2023b) show that at fixed
temperature the auroral-to-nebular line ratios for [O II] and [S II] can serve
as a density diagnostic over a large range of electron density, $10^{2}$
cm${}^{-3}<$ $n_{\rm{e}}$ $<$ $10^{6}$ cm-3. For a uniform density H ii
region, the $n_{\rm{e}}$ returned from both of these diagnostics should be
identical, as long as $n_{\rm{e}}$ is within the sensitivity range of the
diagnostics. However, in the presence of density inhomogeneities, different
density diagnostics can return conflicting values. Even if high density gas
clumps make up a small fraction of the gas, such regions can continue to
contribute to the auroral line emission while no longer contributing
significantly to the nebular lines, since the effects of collisional de-
excitation on the nebular lines will reduce their emissivities relative to the
auroral lines (Rubin, 1989). Because of this, the nebular [S II] and [O II]
lines can reflect the dominant contribution of low-density gas, while the
auroral [S II] and [O II] lines will be sensitive to volume of high-density
gas (Peimbert, 1971; Rubin, 1989; Méndez-Delgado et al., 2023b).

To investigate if the presence of density inhomogeneities could bias our
measured $T_{\rm{e},\text{[{S}{II}]}}$ and $T_{\rm{e},\text{[{O}{II}]}}$ we
compare the observed auroral-to-nebular line ratios [O II] and [S II] to those
predicted using fixed $T_{\rm{e},\text{[{N}{II}]}}$. We show in Figure 11 the
measured [O II]$\lambda\lambda 7320,7330$/$\lambda\lambda 3726,3729$ and [S
II]$\lambda\lambda 4069,4076$/$\lambda\lambda 6716,6731$ line ratios versus
the region’s $T_{\rm{e},\text{[{N}{II}]}}$. We overlay the predicted trends of
auroral-to-nebular line ratios calculated using
$n_{\rm{e}}=10^{2}~{}\mathrm{cm}^{-3}$, $10^{2.5}~{}\mathrm{cm}^{-3}$, and
$10^{3}~{}\mathrm{cm}^{-3}$. We see in Figure 11 that under the assumption
that
$T_{\rm{e},\text{[{N}{II}]}}$$=$$T_{\rm{e},\text{[{S}{II}]}}$$=$$T_{\rm{e},\text{[{O}{II}]}}$,
the measured auroral-to-nebular line ratios could be consistent with
$T_{\rm{e}}$ traced by $T_{\rm{e},\text{[{N}{II}]}}$ but with a higher
electron density than that returned by [S II]$\lambda\lambda 6731,6716$ in the
low-density limit. This suggests that under inhomogenous conditions,
underestimated contributions from $>10^{3}$ cm-3 gas to the nebular [O II] and
[S II] lines could bias the density diagnostics and then the use of
underestimated densities in temperature calculations for ions with low
critical densities like [O II] and [S II] could lead to hotter estimated
temperatures inferred from auroral-to-nebular ratios.

Figure 11: The [S II] and [O II] auroral-to-nebular line ratios vs.
$T_{\rm{e},\text{[{N}{II}]}}$. We show auroral-to-nebular line ratios measured
from [S II] and [O II] against their measured $T_{\rm{e},\text{[{N}{II}]}}$.
We assume equality between the low- ionization zone temperatures and overlay
lines of predicted [S II] and [O II] auroral-to-nebular line ratios for fixed
electron densities $n_{e}=10^{2}\ \mathrm{cm}^{-3}$, $10^{2.5}\
\mathrm{cm}^{-3}$, and $10^{3}\ \mathrm{cm}^{-3}$, vs.
$T_{\rm{e},\text{[{N}{II}]}}$.

####  7.2.2 Correlations between Temperature Differences and Ionization
Parameter

We show in Figures 12 and 13 the comparisons between temperature differences
and the ionization parameter, $U$, traced with [S III]/[S II] or [O III]/[O
II].

Figure 12: Electron temperature differences compared to the H ii region
ionization parameter, $U$, traced by [S III]/[S II]. Top: The $\Delta T_{e}$’s
between the low ionization zone temperatures. Middle: The $\Delta T_{e}$’s
between the low and intermediate ionization zone temperatures. Bottom: The
$\Delta T_{e}$’s between the low, intermediate and high ionization zone
temperatures. We observe significant correlations with $U$, traced by [S
III]/[S II], between the low ionization zone temperatures differences
$\Delta$($T_{\rm{e},\text{[{N}{II}]}}$, $T_{\rm{e},\text{[{S}{II}]}}$) and
$\Delta$($T_{\rm{e},\text{[{S}{II}]}}$, $T_{\rm{e},\text{[{O}{II}]}}$);
between the low and intermediate ionization zones
$\Delta$($T_{\rm{e},\text{[{S}{II}]}}$, $T_{\rm{e},\text{[{S}{III}]}}$); and
the high ionization $\Delta$($T_{\rm{e},\text{[{N}{II}]}}$,
$T_{\rm{e},\text{[{S}{III}]}}$) and $\Delta$($T_{\rm{e},\text{[{S}{II}]}}$,
$T_{\rm{e},\text{[{O}{III}]}}$).

Figure 13: Electron temperature differences compared to the H ii region
ionization parameter, $U$, traced by [O III]/[O II]. The ionization zones
depicted in each row follow those in Figure 12. We observe strong correlations
between all the $\Delta T_{\rm{e}}$ involving $T_{\rm{e},\text{[{O}{III}]}}$.
The Correlations between the low and intermediate ionization zones and $U$,
traced by with [S III]/[S II] and shown in Figure 12, do not appear when
tracing $U$ with [O III]/[O II].

We observe a moderate correlation between
$\Delta$($T_{\rm{e},\text{[{S}{II}]}}$, $T_{\rm{e},\text{[{S}{III}]}}$) and [S
III]/[S II]. The $\Delta$($T_{\rm{e},\text{[{N}{II}]}}$,
$T_{\rm{e},\text{[{S}{II}]}}$) and $\Delta$($T_{\rm{e},\text{[{S}{II}]}}$,
$T_{\rm{e},\text{[{O}{II}]}}$) also show correlations with and [S III]/[S II],
although these are weaker than the correlation involving
$\Delta$($T_{\rm{e},\text{[{S}{II}]}}$, $T_{\rm{e},\text{[{S}{III}]}}$). In
addition, for all temperature differences with $T_{\rm{e},\text{[{O}{III}]}}$
we also observe strong correlations between $\Delta T_{\rm{e}}$ and [S III]/[S
II]. The largest temperature differences are associated with the smallest
values of [S III]/[S II]. As [S III]/[S II] increases, $\Delta T_{\rm{e}}$
converges to $\Delta T_{\rm{e}}=0$. Shown in Figure 13, we present the same
temperature differences but using [O III]/[O II] as a tracer for $U$. All the
correlations with ionization parameter that do not involve
$T_{\rm{e},\text{[{O}{III}]}}$ disappear when using [O III]/[O II] as a
tracer. The remaining $\Delta T_{\rm{e}}$ that include
$T_{\rm{e},\text{[{O}{III}]}}$ show similar correlations as before with
ionization parameter when using [O III]/[O II]. In both cases, the largest
temperature differences occur at the lowest values of ionization parameter,
traced by either [S III]/[S II] or [O III]/[O II].

Our correlations with ionization parameter tracers are similar to the results
presented in Yates et al. (2020). When comparing the temperatures of the low
and high ionization zone for oxygen, Yates et al. (2020) observed an increase
in the $T_{\rm{e},\text{[{O}{III}]}}$/$T_{\rm{e},\text{[{O}{II}]}}$ ratio that
is anti-correlated with the ratio of $\rm{O}^{2+}/\rm{O}^{+}$ (which closely
follows the ionization parameter traced by [O III]/[O II]). Our H ii regions
are all likely to be in the relatively high metallicity regime, where
$\rm{O}^{+}$ should be the dominant ionization state of oxygen. Here, the
average electron temperature will be best described by the auroral-to-nebular
line ratio of [O II] (with caution about density inhomogeneities, as
previously noted). Nevertheless, emission from [O III]$\lambda 4363$ can still
be produced, albeit more weakly, and given the exponential temperature
dependence, it will be biased towards hotter gas. Therefore, if there are
temperature inhomogenieities, the [O III] temperatures may reflect a small
amount of hot, high ionization gas and may not agree with the auroral-to-
nebular line ratio of [O II]. Yates et al. (2020) predict that regions with
$T_{\rm{e},\text{[{O}{III}]}}$ $>$ $T_{\rm{e},\text{[{O}{II}]}}$ will be
$\rm{O}^{+}$ dominant, i.e. $\rm{O}^{+}/\rm{H}^{+}>\rm{O}^{2+}/\rm{H}^{+}$ .
Because [O III]/[O II] $\propto$ $\rm{O}^{2+}/\rm{O}^{+}$, we would expect to
see the largest deviations in $T_{\rm{e},\text{[{O}{III}]}}$ at the lowest [O
III]/[O II]. To summarize, Yates et al. (2020) postulate that the differences
between hotter high ionization zone temperatures,
$T_{\rm{e},\text{[{O}{III}]}}$ and the low ionization zone temperature,
$T_{\rm{e},\text{[{O}{II}]}}$ will increase with decreasing ionization
parameter. These trends should in theory also be evident for sulfur, although
these correlations were not explored by Yates et al. (2020). Similar trends of
temperature differences associated with different ionization states of the gas
have been discussed by Berg et al. (2020).

Figure 14: The $T_{\rm{e}}$–$T_{\rm{e}}$ relation between (left)
$T_{\rm{e},\text{[{O}{III}]}}$–$T_{\rm{e},\text{[{O}{II}]}}$ and (right)
$T_{\rm{e},\text{[{S}{III}]}}$–$T_{\rm{e},\text{[{S}{II}]}}$ colored by the
value of [O III]/[O II] and [S III]/[S II]. In both panels, the largest ratios
of $T_{\rm{e},\text{[{X}{III}]}}$/$T_{\rm{e},\text{[{X}{II}]}}$, where X is
either S or O, occur for the lowest values of ionization parameter, as
indicated by the value of the line ratio, as expected according to the results
of Yates et al. (2020).

In Figure 14, we plot the $T_{\rm{e}}$–$T_{\rm{e}}$ between
$T_{\rm{e},\text{[{O}{III}]}}$–$T_{\rm{e},\text{[{O}{II}]}}$ and
$T_{\rm{e},\text{[{S}{III}]}}$–$T_{\rm{e},\text{[{S}{II}]}}$ colored by the
value of [O III]/[O II] and [S III]/[S II]. We also annotate the plot
according to the schematics from Figure 3 of Yates et al. (2020). For both
oxygen and sulfur, the largest ratios of
$T_{\rm{e},\text{[{X}{III}]}}$/$T_{\rm{e},\text{[{X}{II}]}}$, where X is
either S or O, occur for the lowest values of ionization parameter, as
indicated by the value of the line ratio. This suggests we are observing a
similar correlation with ionization parameter as postulated by Yates et al.
(2020).

If we consider the similarity between our results and Yates et al. (2020) as
evidence for correlations between $\Delta T_{\rm{e}}$ and $U$, then why is it
that the correlation between $\Delta$($T_{\rm{e},\text{[{S}{II}]}}$,
$T_{\rm{e},\text{[{S}{III}]}}$) and $U$ traced by [S III]/[S II]is not evident
when using [O III]/[O II] as a tracer for $U$? It might be the case that [S
III]/[S II] and [O III]/[O II] do not change with $U$ in similar ways, as [S
III] and [O III] arise from different ionization zones and conditions. It is
possible that density inhomogeneities may be playing a role, both [S II] and
[O II] have similar low critical densities, unlike either [S III] and [O III].
It remains unclear why some $\Delta T_{\rm{e}}$ vs. $U$ trends disappear.

####  7.2.3 Correlations between Temperatures Differences with Stellar Mass
and Age

Next we compare temperature differences to stellar mass and ages from compact
stellar clusters and associations matched to our H ii regions. As presented in
Appendix C Figures 32–34, we find no correlations between $\Delta T_{\rm{e}}$
with cluster mass, cluster age, or with association age.

We do observe a weak, positive correlation, between the most reliable $\Delta
T_{\rm{e}}$ indicator without density inhomogeneity issues,
$\Delta$($T_{\rm{e},\text{[{N}{II}]}}$, $T_{\rm{e},\text{[{S}{III}]}}$) with
the association mass in Figure 15. Here we see that
$T_{\rm{e},\text{[{N}{II}]}}$ is cooler than $T_{\rm{e},\text{[{S}{III}]}}$ at
lower association mass. As mass increases, we observe that the value of
$\Delta$($T_{\rm{e},\text{[{N}{II}]}}$, $T_{\rm{e},\text{[{S}{III}]}}$)
decreases until $T_{\rm{e},\text{[{N}{II}]}}$ becomes larger than
$T_{\rm{e},\text{[{S}{III}]}}$. This correlation with association mass spans a
small range in temperature differences $\sim$ [$-0.2$,$0.2$] $\times 10^{4}$ K
with association mass between $10^{3}$ M⊙ and $10^{6}$ M⊙.

Figure 15: Electron temperature differences compared to stellar association
mass. The ionization zones depicted in each row follow those in Figure 12.
Although the correlations are insignificant according to their p-value, we
observe potential, weak, correlations between
$\Delta$($T_{\rm{e},\text{[{N}{II}]}}$, $T_{\rm{e},\text{[{S}{III}]}}$), and
all $\Delta T_{\rm{e}}$ involving $T_{\rm{e},\text{[{O}{III}]}}$, with stellar
association mass.

As discussed in Section 6.5, we assigned the largest mass stellar population
to the H ii region if the region was matched to more than one association. To
see if this choice has any impact on the observed correlation, we plot in
Figure 16 the $\Delta$($T_{\rm{e},\text{[{N}{II}]}}$,
$T_{\rm{e},\text{[{S}{III}]}}$) against the sum of the matched association
masses. Using the total masses, we observe a negligible change in the
correlation. The p-value is increased from $p=0.001$ to $p=0.004$ and the
Spearman rank remains unchanged, $\rho=0.37$. To investigate the impact that
the error has on the correlation, we construct a distribution Spearman rank
coefficients by perturbing each data point by its associated error and
remeasuring the correlation. From this we measure a Spearman rank coefficient
$0.36\pm 0.05$, indicating that the error has negligible impact on the
correlation.

Figure 16: The $T_{\rm{e},\text{[{N}{II}]}}$$-$$T_{\rm{e},\text{[{S}{III}]}}$
temperature difference versus the total stellar association mass. The black-
points are individual H ii region association masses. In red we show the
average temperature differences and H ii region association masses calculated
in Log10($M_{\odot}$)=0.5 bins. We also show a fit and $1\sigma$ fit in grey.
Compared to using the largest mass measurement, see middle row Figure 15, the
p-value using total stellar mass p-value is higher, indicating a less
statistically significant correlation.

It is possible that the correlation between
$\Delta$($T_{\rm{e},\text{[{N}{II}]}}$, $T_{\rm{e},\text{[{S}{III}]}}$) and
association mass may result from biases introduced by an under-sampling of
initial mass function. Similar to our study, Scheuermann et al. (2023) matched
the stellar association catalog (Larson et al., 2023) to H ii regions in the
Nebular catalog (Kreckel et al., 2019; Groves et al., 2023). Instead of
including all masses measured for the associations, which we do in this work,
Scheuermann et al. (2023) implement a cut-off and assume masses
$<10^{4}M_{\odot}$, as masses below the threshold do not sample the IMF (da
Silva et al., 2012). We perform no such cutoff in this study. The correlation
between $\Delta$($T_{\rm{e},\text{[{N}{II}]}}$,
$T_{\rm{e},\text{[{S}{III}]}}$) and $M_{\mathrm{Association}}$ correlation
includes many regions with masses $<10^{4}M_{\odot}$, where the IMF may not be
fully sampled. Despite this, we do not fully dismiss this correlation;
however, the effects that undersampling the IMF could have on the correlation
with $\Delta$($T_{\rm{e},\text{[{N}{II}]}}$, $T_{\rm{e},\text{[{S}{III}]}}$)
warrants further investigation with a larger sample size on order to expand
the dynamic range of association mass.

####  7.2.4 Correlations between Temperatures Differences and Molecular Gas
Properties

We compare temperature differences to the properties of the molecular gas
derived from CO emission measured within the projected boundaries of our
sample of H ii regions. For the comparisons to $I_{\rm CO}$, shown in Figure
35 of Appendix E , and comparisons to $T_{\rm{peak}}$, shown in Figure 36 of
Appendix E, we observe no correlations with temperature differences between
the low, intermediate, and high ionization zones.

We show in Figure 17 the temperature differences compared to the CO velocity
dispersion ($\sigma_{v,\mathrm{CO}}$). We observe a moderate correlation
between $\sigma_{v,\mathrm{CO}}$ and $\Delta$($T_{\rm{e},\text{[{N}{II}]}}$,
$T_{\rm{e},\text{[{S}{III}]}}$). The values of
$\Delta$($T_{\rm{e},\text{[{N}{II}]}}$, $T_{\rm{e},\text{[{S}{III}]}}$) over
the range of $\sigma_{v,\mathrm{CO}}$ is small, only encompassing
$\Delta$($T_{\rm{e},\text{[{N}{II}]}}$, $T_{\rm{e},\text{[{S}{III}]}}$) from
$-$1000 K to 2000 K. Given the possible correlation with association mass, and
this correlation with $\sigma_{v,\mathrm{CO}}$, we find it intriguing that
these correlations between the low and intermediate ionization zone are only
seen in $\Delta$($T_{\rm{e},\text{[{N}{II}]}}$,
$T_{\rm{e},\text{[{S}{III}]}}$). It is possible that the low scatter in
$T_{\rm{e}}$–$T_{\rm{e}}$ relationship between $T_{\rm{e},\text{[{S}{III}]}}$
and $T_{\rm{e},\text{[{N}{II}]}}$ in the presence of density inhomogeneities,
contrary to the temperatures from [S II] and [O II], allow for better insight
into these underlying trends.

We observe that all the correlations of $\Delta T_{\rm{e}}$ involving
$T_{\rm{e},\text{[{O}{III}]}}$ and $\sigma_{v,\mathrm{CO}}$ have large,
negative, $\rho$ values, but are all insignificant according to their
p-values. The high Spearman rank values for the comparisons between $\Delta
T_{\rm{e}}$ and $\sigma_{v,\mathrm{CO}}$ appears to be largely driven by the
highest values of $T_{\rm{e},\text{[{O}{III}]}}$. Despite the fact that the
correlations are not strong, it is clear that the high
$T_{\rm{e},\text{[{O}{III}]}}$ regions go along with high CO velocity
dispersion.

Figure 17: Electron temperature differences compared to the CO velocity
dispersion, $\sigma_{v,\rm{CO}}$. The order of the panels follow those in
Figure 12. We observe a weak correlation between
$\Delta$($T_{\rm{e},\text{[{N}{II}]}}$, $T_{\rm{e},\text{[{S}{III}]}}$).
Although insignificant according to their p-values, we also observe that
highest $T_{\rm{e},\text{[{O}{III}]}}$ values are associated with the highest
$\sigma_{v,\rm{CO}}$.

The CO velocity dispersion can be enhanced by low-velocity shocks originating
from the interaction of molecular gas with late-time-expansion of supernovae
remnants (see Koo et al., 2001; Zhou et al., 2023), as well as from
interaction with low-velocity shocks from pressure-radiation driven H ii
region expansion (Hill & Hollenbach, 1978; Kothes & Kerton, 2002; Watkins et
al., 2023). These low-velocity shocks are also predicted to enhance
$T_{\rm{e},\text{[{O}{III}]}}$. Shock modeling have shown that outward
expanding low-velocity shocks can create conditions, such as pockets of high
post-shock temperature, where the strength of [O III]$\lambda 4363$ emission
will be enhanced compared to the little to no increase in emission from [O
III]$\lambda 5007$ (Peimbert et al., 1991; Binette et al., 2012; Méndez-
Delgado et al., 2021). The combination of large temperature differences with
$T_{\rm{e},\text{[{O}{III}]}}$ and high $\sigma_{v,\mathrm{CO}}$ suggests that
we may be observing the effects of low-velocity shocks. Motivated by this
scenario, we search the H ii regions for evidence of low-velocity shocks in
the following section.

####  7.2.5 Investigating H II Regions for Presence of Low-Velocity Shocks

The correlations of $\Delta T_{\rm{e}}$ involving the high ionization zone and
H ii region properties: [S III]/[S II], [O III]/[O II], and the CO velocity
dispersion appears to be driven by the presence of regions with high
$T_{\rm{e},\text{[{O}{III}]}}$, high $\sigma_{\mathrm{v,CO}}$ and low $U$. One
potential explanation is the presence of low velocity shocks enhancing
$T_{\rm{e},\text{[{O}{III}]}}$ (Peimbert et al., 1991; Binette et al., 2012;
Méndez-Delgado et al., 2021). To test this explanation, we search for evidence
of shocks in enhanced optical line ratios and line broadening.

When shocks collide with and compress gas, the ionization parameter of the gas
is reduced leading to partially ionized zones of enhanced nebular emission of
low-ionization species such as S+ and O relative to H$\beta$ and/or H$\alpha$
(Dopita & Sutherland, 1996; Allen et al., 2008). In Figure 18, we show the [S
II]/H$\alpha$ vs. [O I]/H$\alpha$ ratios, color-coded by
$T_{\rm{e},\text{[{O}{III}]}}$, for regions with measured [O III]$\lambda
4363$. We find that H ii regions with hotter $T_{\rm{e},\text{[{O}{III}]}}$
tend to populate a region with enhanced [S II]/H$\alpha$ and [O I]/H$\alpha$
ratios which suggests that these regions may host a partially ionized zone due
to shocks. We note here that it is possible that harder photons and X-rays
produced by X-ray binaries also enhance [S II] and [O I] relative to the
Balmer emission (Abolmasov et al., 2007; Grisé et al., 2008). Furthermore,
X-rays would provide high energy photons able to boost [O III]$\lambda 4363$.
However, lacking the high spatial resolution X-ray imaging of these H ii
regions, exploration of X-ray contributions to [O III]$\lambda 4363$ emission
is beyond the scope of this paper.

Another tracer that may indicate presence of shocks is the He II $\lambda
4686$/H$\beta$ ratio (Allen et al., 2008) though it is also sensitive to the
shape of the Lyman continuum below 228 Å (Garnett et al., 1991; Guseva et al.,
2000; Allen et al., 2008). Wolf-Rayet (WR) stars are capable of releasing
photons able to produce He II $\lambda 4686$. WR stars host stellar winds and
can be a source of high energy photons, E$>54$ eV, capable of doubly ionizing
Helium. We measured the He II$\lambda 4686$/H$\beta$ in H ii regions with
$T_{\rm{e},\text{[{O}{III}]}}$ detections and found that these regions exhibit
an average He II $\lambda 4686$/H$\beta=4.1\pm 1.6\%$. This value is within
the expected range of He II $\lambda 4686$/H$\beta$ values for H ii regions
with WR, 0.04%-7% (Guseva et al., 2000; Thuan & Izotov, 2005; Mayya et al.,
2023), and 100 km s-1 shocks, 4%-6% (Allen et al., 2008). We visually
inspected the spectra of regions with measured $T_{\rm{e},\text{[{O}{III}]}}$
for the characteristic red/blue bump associated with the presence of WR stars.
We found the blue bump in only 2 regions with measurable [O III]$\lambda
4363$, both of which have $T_{\rm{e},\text{[{O}{III}]}}$ $<10^{4}$ K.
Nevertheless, the He II $\lambda 4686$/H$\beta$ for the high
$T_{\rm{e},\text{[{O}{III}]}}$ regions lends evidence that these regions may
host shocks or undetected WR stars.

We also searched for kinematic signatures of shocks, but find no clear
kinematic evidence. Shocks can imprint asymmetries and/or broad emission near
the base of an emission line. We inspected the fit residuals of [O
III]$\lambda 5007$ in H ii regions with measured
$T_{\rm{e},\text{[{O}{III}]}}$ and found no evidence of line broadening due to
shocks. Next, we inspected the measured line widths of [O III]$\lambda 5007$
in these regions. We compared the line-spread-function corrected [O
III]$\lambda 5007$, as measured by MUSE, velocity dispersion,
$\sigma_{v,\lambda 5007}$ vs. the CO velocity dispersion,
$\sigma_{v,\mathrm{CO}}$. The line-widths of the optical and CO emission for
regions with high $T_{\rm{e},\text{[{O}{III}]}}$ are comparable to regions
with low, or not detected $T_{\rm{e},\text{[{O}{III}]}}$. Despite this, the
absence of these features may only exclude the presence of high-velocity
shocks.

If the regions do host low-velocity shocks, then their impact on the line
width of the optical emission may be too small to be resolved given the MUSE
resolution; $\sim 70$ km s-1 at $\lambda=5007$ Å. The high velocity
resolution, $2.5~{}\mathrm{km}\mathrm{s}^{-1}$ of the PHANGS-ALMA data, are
more sensitive than MUSE to low-velocity shocks, and are the strongest
evidence for their presence in the high $T_{\rm{e},\text{[{O}{III}]}}$ H ii
regions.

Figure 18: We compare the relative strength of the low-ionization species [O
I] and [S III] to Balmer emission. We color-code each data point by the value
of its measured [O III] temperature. The hottest
$T_{\rm{e},\text{[{O}{III}]}}$ regions also have the highest line ratios. The
higher line ratios are indicative of the existence of a partially-ionized zone
that may be due to the presence of shocks or harder photons such as X-rays
(Grisé et al., 2008).

##  8 Discussion

###  8.1 Electron Density Inhomogeneities

Similar to the results of Méndez-Delgado et al. (2023b), described in Section
7.2.1, we have found that the presence of electron density inhomogeneities may
cause the temperatures measured from [O II] and [S II] to be biased hotter
compared to those measured using the [N II] auroral lines.

The critical densities of the nebular lines of [S II]$\lambda\lambda
6716,6731$ and more so [O II]$\lambda\lambda 3726,3729$ are low enough that at
densities $n_{\rm{e}}>10^{3}\ \mathrm{cm}^{-3}$ the lines will undergo
increased collisional de-excitation. In the presence of density
inhomogoneities above this value, the low critical densities will reduce the
emissivity of the nebular emission lines from [S II] and [O II]. Because of
this, the nebular diagnostic lines of [S II] and [O II] will mainly describe
gas component with $n_{\rm{e}}<10^{3}\ \mathrm{cm}^{-3}$. This biases the [S
II] and [O II] density diagnostic to return electron densities that are lower
than the true average density of the H ii region. Since the auroral line
critical densities are far higher, this also makes the measured [S II] and [O
II] temperatures appear to reflect hotter values.

This effect has been observed in many studies. Densities measured in Milky Way
H ii regions (the Orion Nebula, NGC 3604 and NGC 3576) using [Cl III], which
is sensitive to greater densities than [S II], routinely show that the [S II]
diagnostic returns lower values than [Cl III]. (Pogge et al., 1992; García-
Rojas & Esteban, 2007; Núñez-Díaz et al., 2013; Weilbacher et al., 2015).
Although densities derived from the [S II] doublet are commonly used in the
literature due to their strengths and their insensitivity to dust extinction
the [S II]$\lambda 6731/\lambda 6716$ ratio is less sensitive to density than
[Cl III], [Ar IV] and [Fe III] diagnostics when $n_{\rm{e}}>10^{3}$ cm-3 (see
Figure 2 in Méndez-Delgado et al., 2023b). If the equality
$T_{\rm{e},\text{[{O}{II}]}}$$\approx$$T_{\rm{e},\text{[{S}{II}]}}$$\approx$$T_{\rm{e},\text{[{N}{II}]}}$
predicted by photoionization models is true, then the auroral-to-nebular line
ratios would suggest a factor of 10 higher density than the [S II] doublet,
but neither may well represent the true average density of the region.

Electron density variations may arise from shocks, turbulence, and pre-
existing non-uniform structure in the ISM (Hill & Hollenbach, 1978; Dopita &
Sutherland, 1996; Allen et al., 2008). Jin et al. (2022) have extended
photoionization modeling of ionized nebulae to more complex geometries.
Starting with an initial clumpy ISM, ionizing photons will pass through
diffuse regions more readily than denser clumps. The resulting photoionized
region will exhibit fluctuations in density and irregular geometry as opposed
to a uniform density and spherical morphology. The ionization parameter in the
dense clumps will be relatively lower than other regions of the nebulae due to
the increased density. At these locations, the emission of low-ionization
species, including [S II] and [O II], will be enhanced compared to higher
ionization species. Due to the higher critical densities of the auroral lines
of these ions, the emissivity of the auroral lines will be greater than those
of the nebular lines in the high-density portion the nebula. In this scenario,
the nebular density diagnostics can return an average density that traces the
low-density portion of the nebula. In doing so, the value of $n_{\rm{e}}$
returned by the nebular diagnostics may inaccurately describe the ionizing
conditions of the high-density clumps where the auroral line emissivities are
greater than the nebular lines, and overestimate [S II] temperatures.

Density inhomogeneities have been reported in studies of highly resolved local
HII regions like Orion where the inhomogeneities can be spatially resolved
(Baldwin et al., 1991; Pogge et al., 1992; Weilbacher et al., 2015; McLeod et
al., 2016; O’Dell et al., 2017). Weilbacher et al. (2015) mapped the spatial
variation of density in the Orion nebulae and found variations of density
between 500 cm-3 and in excess of 10,000 cm-3. A maximum of 25,000 cm-3 is
measured using [Cl III] at the location of the ionization front in the “Orion
S” area. Density inhomogeneities, associated with turbulence driven velocity
fluctuations, have been invoked as one mechanism to generate surface
brightness fluctuations within the Orion Nebula (Kainulainen et al., 2017).

The Orion Nebula is not particularly comparable to the H ii regions studied
here due to the difference in scales (i.e., Orion is more compact) and
resolution. A subset of our H ii regions, as shown in the comparison between
the MUSE H ii regions mask in Appendix C, are unresolved clusters of
individual regions. Measurements of density inhomogeneities using density
diagnostics besides the nebular [O II] and [S II] doublets for extragalactic H
ii regions that more closely match our sample are rare and require deep, high
$S/N$ spectra (Méndez-Delgado et al., 2023b). One consequence of the lack of
different diagnostics is that many studies will often assume a fixed density
of $100$ cm-3 when either [S II], or [O II], return a density in the low-
density limit (e.g. Kreckel et al., 2019). Studies using mid-infrared
observations have shown this latter assumption could be incorrect, as
densities up to $1000~{}\mathrm{cm}^{-3}$ have been measured using the [S
III]$\lambda 18.7/33.5\ \mu$m density diagnostic (see Rubin et al., 2016), and
indicate that density inhomogeneities are present in extra-galactic H ii
regions.

The consistency of auroral-to-nebular ratios with $n_{\rm{e}}\sim 1000$ cm-3
assuming $T_{\rm{e}}=$ $T_{\rm{e},\text{[{N}{II}]}}$ in Figure 11 supports
this picture that inhomogeneities must be considered when deriving
temperatures from the [O II] and [S II]. As a result, we consider
$T_{\rm{e},\text{[{N}{II}]}}$ a more reliable indicator of the low-ionization
zone temperature due to its relative insensitivity to density.

###  8.2 $T_{\rm{e},\text{[{N}{II}]}}$ and $T_{\rm{e},\text{[{S}{III}]}}$ as
Accurate Tracers of HII Region Temperatures

The low scatter and overlap with literature $T_{\rm{e}}$–$T_{\rm{e}}$ trends
between $T_{\rm{e},\text{[{N}{II}]}}$ and $T_{\rm{e},\text{[{S}{III}]}}$ in
Figure 9 indicate that these temperatures may be optimal tracers of H ii
region temperatures, more so than $T_{\rm{e},\text{[{S}{II}]}}$,
$T_{\rm{e},\text{[{O}{II}]}}$ and $T_{\rm{e},\text{[{O}{III}]}}$.

The tight relationship between $T_{\rm{e},\text{[{N}{II}]}}$ and
$T_{\rm{e},\text{[{S}{III}]}}$ measured in this study has been observed in the
some studies (Berg et al., 2015; Croxall et al., 2016; Berg et al., 2020;
Zurita et al., 2021), especially for the coolest, highest metallicity regions.
It is this metallicity regime where it is expected that the widely used [O
III] auroral line temperatures will be most difficult to observe, as evident
in the small number of $T_{\rm{e},\text{[{O}{III}]}}$ measurements relative to
$T_{\rm{e},\text{[{N}{II}]}}$ and $T_{\rm{e},\text{[{S}{III}]}}$ measurements
in our study. We showed that $\Delta$($T_{\rm{e},\text{[{N}{II}]}}$,
$T_{\rm{e},\text{[{S}{III}]}}$) does have a small correlation with the H ii
region’s matched stellar association mass and molecular gas velocity
dispersion, though, the magnitude of the temperature difference for these
correlations span a narrow range of $\Delta T_{\rm{e}}$, indicative that
$T_{\rm{e},\text{[{N}{II}]}}$ and $T_{\rm{e},\text{[{S}{III}]}}$ change in
similar ways with H ii region properties, unlike correlations with [O III]
temperatures, which show very large differences with other ionization zones.
In fact, it is tantalizing that the observed correlations between the low and
intermediate ionization zones mostly involve $T_{\rm{e},\text{[{N}{II}]}}$ and
$T_{\rm{e},\text{[{S}{III}]}}$. These underlying trends may be more easily
revealed due to the well behaved relationship between
$T_{\rm{e},\text{[{N}{II}]}}$ and $T_{\rm{e},\text{[{S}{III}]}}$ in
$T_{\rm{e}}$–$T_{\rm{e}}$ relations in the presence of inhomogeneous
conditions. Unlike [O II] and [S II], [N II] is less sensitive to density
inhomogeneities due to the the higher critical density of the atomic levels
involved. The [N II] auroral line is also less contaminated by blended
emission (Berg et al., 2020; Méndez-Delgado et al., 2023b).
Additionally,temperatures from [N II] were shown to strongly correlate with O
II RL temperatures measured for the low-ionization zone Méndez-Delgado et al.
(2023a).

Judged from Oxygen CEL and RL temperatures, the high ionization zone is
expected to be most affected by temperature fluctuations due to its proximity
to sources of feedback (Méndez-Delgado et al., 2023a). To what degree the
intermediate ionization zone temperatures are affected is unclear (Méndez-
Delgado et al., 2023b). In Figure 9 we showed that the relationship between
$T_{\rm{e},\text{[{N}{II}]}}$ and $T_{\rm{e},\text{[{S}{III}]}}$ generally
follow trends predicted from photoionization models (Morisset, 2009; Vale
Asari et al., 2016) with $t^{2}=0$. The overlap with our data indicates that
$T_{\rm{e},\text{[{S}{III}]}}$ may be minimally affected by temperature
inhomogeneities. Díaz & Zamora (2022) provide additional arguments for the use
of $T_{\rm{e},\text{[{S}{III}]}}$ over $T_{\rm{e},\text{[{O}{III}]}}$,
including: 1) the emission lines of [S III] have a lower exponential
dependence on electron temperature and 2) because [S III] overlaps gas volumes
containing both $\mathrm{O}^{2+}$ and $\mathrm{O}^{+}$ that
$T_{\rm{e},\text{[{S}{III}]}}$ can be representative of the entire H ii
region.

The behaviour between $T_{\rm{e},\text{[{N}{II}]}}$ and
$T_{\rm{e},\text{[{S}{III}]}}$ across several H ii region properties and
photoionization model $T_{\rm{e}}$–$T_{\rm{e}}$ suggests prioritization of [N
II] and [S III] temperatures over those from [O II], [S II], and [O III] for
regions in our sample. We will test this recommendation in an upcoming work,
where we will compare “direct” metallicities, derived using [N II] and [S III]
temperatures, to several calibrated methods (Rickards Vaught et al., in prep).

###  8.3 The High-Ionization Zone Temperature Excess

We found H ii regions with high $T_{\rm{e},\text{[{O}{III}]}}$, enhanced
velocity dispersion in the surrounding molecular gas and low-ionization
parameter. We investigated these H ii regions for enhanced low-ionization
species emission line ratios [S II]/H$\alpha$, [O I]/H$\alpha$ and the high-
ionization ratio He II/H$\beta$. We found that the high
$T_{\rm{e},\text{[{O}{III}]}}$ regions exhibited enhanced low-ionization line
ratios suggestive of shock ionization (or possibly X-ray ionization). We also
found that the He II/H$\beta$ ratios for these regions are within the range of
those expected from shock velocities $<100$ km s-1 (Allen et al., 2008) and/or
WR stars (Guseva et al., 2000; Thuan & Izotov, 2005; Mayya et al., 2023),
though we only find characteristic red/blue bumps associated with the presence
of WR stars, in 2 H ii regions neither of which had elevated
$T_{\rm{e},\text{[{O}{III}]}}$. Absent the WR signatures, we are motivated to
explore shocks as enhancers of the [O III] temperature.

We did not find any kinematic signatures of shocks in the optical emission
line profiles. We determine that if shocks are present and broadening the CO
emission, then the shock velocities are too low to be resolved by the MUSE
spectral resolution. Either way, high-velocity shocks are not expected to
effectively boost $T_{\rm{e},\text{[{O}{III}]}}$ (Méndez-Delgado et al.,
2021). Despite the uncertainty in whether shocks are present in the high [O
III] temperature regions, we can discuss the plausibility that low-velocity
shocks are the cause for excess $T_{\rm{e},\text{[{O}{III}]}}$.

####  8.3.1 Shock-enhanced [O III] Temperature

Low velocity, $<$ 100 km s-1, shocks can increase [O III]$\lambda 4363$ while
leaving [O III]$\lambda 5007$ unchanged (Peimbert et al., 1991; Binette et
al., 2012; Méndez-Delgado et al., 2021). For a sample of giant H ii regions,
Binette et al. (2012) measured $T_{\rm{e},\text{[{O}{III}]}}$ up to 6000 K
higher than $T_{\rm{e},\text{[{S}{III}]}}$ in regions with
$T_{\rm{e},\text{[{S}{III}]}}$ $<10^{4}$ K. To explore if shocks were boosting
their [O III] temperatures relative to [S III], they model outward expanding
shocks, mimicking those generated by stellar winds, by combining
shock$+$photoionization models, with increasing shock velocities (analogous to
increasing the post-shock temperature from 1.6$\times 10^{4}$ K to 7.2$\times
10^{4}$ K) between 20 and 60 km s-1. Additionally, the models span 5 different
metallicities between Z$=$0.01 Z⊙ to Z$=$1.6 Z⊙. Comparing the average
properties of between the lowest and highest shock velocity models, Binette et
al. (2012) found that mean doubly ionized fraction of oxygen,
$\mathrm{O}^{2+}$/$\mathrm{O}$ , decreases while at the same time leaving the
doubly ionized fraction of sulfur, $\mathrm{S}^{2+}$/$\mathrm{S}$, unchanged.
The imbalance the ionization fraction between oxygen and sulfur means that
hotter, post-shock gas contributes proportionally more to the observed [O III]
emission than [S III]. Because of the exponential sensitivity to temperature,
the [O III] auroral line will be enhanced, tracing the hotter post-shock
temperature rather than the local photoionized nebula temperature returned by
$T_{\rm{e},\text{[{S}{III}]}}$. Furthermore, the highest metallicity models
show the largest temperature differences between
$T_{\rm{e},\text{[{S}{III}]}}$ and $T_{\rm{e},\text{[{O}{III}]}}$, up to $\sim
7000$ K. Although this difference is 3000 K lower than our largest measured
$\Delta$($T_{\rm{e},\text{[{S}{III}]}}$, $T_{\rm{e},\text{[{O}{III}]}}$), a
complete understanding of the degree of enhancement of [O III]$\lambda 4363$
with velocity shocks require more complex 3D hydro-dynamical simulations
(Binette et al., 2012). Despite this, one extreme example of shock impact on
electron temperatures was observed in the outflow of the H ii region Sh 2-129
(Corradi et al., 2014). This outflow, with velocity $\approx 100$
$\mathrm{km}\mathrm{s}^{-1}$, exhibits $T_{\rm{e},\text{[{O}{III}]}}$
$=55,000$ K and $T_{\rm{e},\text{[{O}{II}]}}$ $\sim 20,000$ K. Much larger
than $\Delta T_{\rm{e}}$ observed in this study.

It has also been shown that shocks driven by pressure–radiation H ii region
expansion, and SNR, impact the surrounding cold molecular and ionized gas. If
the velocity of the expansion is greater than the sound speed ionized gas
$\sim 10$ km s-1, a layer of shocked H gas will form in-between the expanding
ionization front and a surrounding molecular gas (Hill & Hollenbach, 1978;
Kothes & Kerton, 2002; Tremblin et al., 2014; Watkins et al., 2023). The
impact of shock interaction with molecular gas has been studied in 18 galaxies
observed as part of PHANGS-ALMA (Leroy et al., 2021a, b), where Watkins et al.
(2023) identified hundreds of super bubbles (i.e. pockets of expanding gas
arising as byproducts of feedback). The superbubbles were identified using
spatial correspondence between CO shells and stellar populations contained in
PHANGS-HST catalogs (Lee et al., 2022; Thilker et al., 2022; Larson et al.,
2023). Due to the ALMA spatial resolution (50–150 pc), Watkins et al. (2023)
measured the expansion velocity only for the largest superbubbles. Assuming
the CO expansion velocity is equal to the shock velocity, Watkins et al.
(2023) measure the velocity of the approaching/receding CO shells and
determine an average, line-of-sight, expansion velocity of $v_{\rm exp}=$9.8
$\pm$ 4.3 km s-1. Although this value is similar to the sound speed of ionized
gas, superbubbles exhibit asymmetries in their morphology, and the velocity
can potentially reach up to a few tens of km s-1 depending on the conditions
of the gas, source or energy injection, and age (Watkins et al., 2023). 1D
models indicate that, at minimum, H ii regions exhibit expansion speeds of a
$\sim$few km s-1 (Tremblin et al., 2014). This suggests that H ii regions
expand with a large range of velocities.

As for ionized gas, Egorov et al. (2023) identify in the PHANGS-MUSE galaxies
more than 1400 regions of ionized gas with elevated intrinsic H$\alpha$
velocity dispersions $>45$ km s-1, and, under the assumption that these
regions are undergoing expansion, Egorov et al. (2023) infer expansion
velocities between $v_{\rm{exp}}=$ 10–40 km s-1 (see also Egorov et al., 2014,
2017; Cosens et al., 2022). The ubiquity of H ii region expansion, as well as
their effects on the surrounding molecular gas, make them good candidates for
drivers of low-velocity shocks capable of boosting [O III] temperatures.

####  8.3.2 Temperature Inhomogeneities

Temperature fluctuations within H ii regions are another potential explanation
for the high [O III] temperature. As discussed in Section 1, RL emissivities
have a linear sensitivity to temperature rather than the exponential
sensitivity of CELs. In the presence of temperature fluctuations CELs will
return higher estimates of temperature. Méndez-Delgado et al. (2023a) have
recently shown that differences between [O III] and [N II] temperatures are
strongly correlated with the temperature fluctuations parameter, $t^{2}$, of
the highly ionized gas. This suggests that $T_{\rm{e},\text{[{O}{III}]}}$ is
likely to overestimate the representative H ii region temperature. The
observed excess in $T_{\rm{e},\text{[{O}{III}]}}$ are likely to be produced by
phenomena other than those commonly observed in H ii regions. Simple models
from Binette et al. (2012) and the observations from Méndez-Delgado et al.
(2023a) show that the effect from $t^{2}>0$ is more pronounced in lower-
metallicity/high ionization parameter regions. Furthermore, the temperature
excess $T_{\rm{e},\text{[{O}{III}]}}$ relative to the other ionization zones
is too large in our observations to be caused solely by inhomogeneities. Using
our measured $\Delta$($T_{\rm{e},\text{[{N}{II}]}}$,
$T_{\rm{e},\text{[{O}{III}]}}$), we can infer from Méndez-Delgado et al.
(2023a) that $t^{2}>0.2$ for $\Delta$($T_{\rm{e},\text{[{N}{II}]}}$,
$T_{\rm{e},\text{[{O}{III}]}}$) $>5000$ K. These values are much higher than
what has been observed for nearby H ii regions (Peña-Guerrero et al., 2012;
Binette et al., 2012; Méndez-Delgado et al., 2023a). Even in the presence of
classical temperature inhomogeneities, a secondary effect would also need to
be included to explain our [O III] temperatures.

####  8.3.3 Potential Observation Bias

Another possibility is that the high [O III] temperatures are statistical
outliers. The galaxies in our sample exhibit strong line oxygen abundances
$\gtrsim$ 8.3 (Kreckel et al., 2019). Because electron temperature is anti-
correlated to the metallicity of the gas, the [O III]$\lambda 4363$
temperatures for these galaxies are expected to be low. For [O III]$\lambda
4363$ to be detectable, the temperature would need to be high, otherwise, we
would likely not detect the auroral line. [O III] temperatures from a lower-
metallicity sample may be compatible with photoionization models and
$T_{\rm{e}}$–$T_{\rm{e}}$ relations. Finally, Rola & Pelat (1994) have shown
that emission line measurements with S/N $<5$ can potentially overestimate the
true intensity by $80\%$. The average S/N [O III]$\lambda 4363$ measured from
our sample is $\sim 4.5$. Nevertheless, why we would measure high [O
III]$\lambda 4363$, even after removal of [Fe II]$\lambda 4360$ contamination,
that also exhibit enhanced molecular gas velocity dispersions is difficult to
explain purely with statistical outliers.

##  9 Conclusion

We presented combined KCWI and MUSE observations of the [N II]$\lambda 5756$,
[O II]$\lambda\lambda 7320,7330$, [S II]$\lambda\lambda 4069,4076$, [O
III]$\lambda 4363$, and [S III]$\lambda 6312$ auroral lines in a sample of 421
H ii regions in 7 nearby galaxies. We compared the derived electron
temperatures and temperature differences between multiple H ii region
ionization zones to several H ii region properties such as electron density,
ionization parameter, molecular gas velocity dispersion, stellar mass and age
obtained from PHANGS observations. We found that:

  * •

Temperatures obtained from [S II] and [O II] are consistent with being
overestimated due to the presence of density inhomogeneities in the H ii
regions. Because of these potential biases, we recommend the use [N II]
temperatures to trace the low ionization zone.

  * •

We found that the [N II] and [S III] temperatures exhibited the lowest scatter
of the $T_{\rm{e}}$– $T_{\rm{e}}$ relations and follow trends predicted from
photoionization models. The well-behaved relationship between [N II] and [S