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

H ii region electron temperatures are a critical ingredient in metallicity determinations, and recent observations have revealed systematic variations in the temperatures measured using different ions. We present electron temperatures (T e ) measured using the optical auroral lines ([N ii]λ5756, [O ii]λ λ7320, 7330, [S ii]λ λ4069, 4076, [O iii]λ4363, and [S iii]λ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) and the Very Large Telescope (Multi-Unit Spectroscopic Explorer). We compare the different T e measurements with H ii region and ISM 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 overestimated due to the presence of electron density inhomogeneities in H ii regions. 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. In agreement with previous studies, the T e–T e between [N ii] and [S iii] temperatures have the lowest observed scatter and 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].

ization 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].

INTRODUCTION
The characterization of abundance variations within galaxies provides insight into the physical processes that drive galaxy and chemical evolution. A galaxy's gasphase 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. 2019Kreckel et al. , 2020van 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 e ) of the gas. Due to their exponential dependence on T e , one of the ways to infer electron temperature is through the auroral-tonebular 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 ∼ 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 e include the ratio between recombination line (RL) emission from H and other species, and which exhibit linear sensitivity to T 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. * ARC DECRA Fellow For ions with optical auroral lines studied in this work, we can measure the temperatures for each ion using the following line ratios: 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. 2015Berg et al. , 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, ∼ 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. 2015Berg et al. , 2020Arellano-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 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. 2015Berg et al. , 2020Wei 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 A 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 3.5. 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.

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

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 sample 2 (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, δ > −30 • . Table 1 presents general properties of these galaxies, including distances, masses, sizes, and the angular resolution of the MUSE data.

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 A with a spectral resolution R ∼ 900, corresponding to a full width at half maximum (FWHM) ∼ 5.1Å (or ∼ 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.

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 ′ × 1 ′ pointings. The median spatial resolution across all pointings is ∼ 50 pc (or ∼ 0.80 ′′ ) with a typical spectral FWHM of ∼ 2.5Å (but varying with wavelength). The data were reduced using the pymusepipe 3 and spectral fitting and analysis was performed using the Data Analysis Pipeline 4 packages described in Emsellem et al. (2022). The individual MUSE pointings were homogenized to a uniform Gaussian point-spreadfunction (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.

ALMA
3 https://pypi.org/project/pymusepipe/ 4 https://gitlab.com/francbelf/ifu-pipeline 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 12 CO, 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 ∼ 1.3 ′′ , similar to the resolution of both KCWI and MUSE. The velocity resolution is 2.5 km s −1 .

HST
The PHANGS-HST survey ) observed 5 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 ) 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.

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) 6 . It was built using the Keck Data Reduction Pipeline Framework package 7 and is a port of the initial IDL reduction pipeline 8 (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Å. 5 Lee et al. (2022)  a From the compilation of Anand et al. (2021) b Derived by Leroy et al. (2021a), using GALEX UV and WISE IR photometry.
d The FWHM of the Gaussian PSF for the homogenized COPT mosaic from PHANGS-MUSE (Emsellem et al. 2022). e The average FWHM of the KCWI PSF for the set of a galaxy's observed pointings. * Denotes galaxies observed with MUSE using ground based adaptive optics.
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 ∼ 9% at λ = 4600Å. The maximum standard deviation within the wavelength range containing the lines used in this analysis is ∼ 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 documentation 9 . 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σ uncertainty cubes, as well as a bad-pixel mask cube. 9 https://kcwi-drp.readthedocs.io/en/latest/sky subtraction.html

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 ′′ × 0.29 ′′ . We reprojected the cubes onto a square pixel grid using the astronomical mosaic software Montage 10 . 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 ′′ × 0.29 ′′ . We validated the flux-conservation in our data before and after reprojection.

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 cat-alogs. 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 widefield 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β 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 ± 17 pixels (or 5 ′′ ) 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 (±1 ′′ 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.

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 KCWI vs. SB 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 ± 0.02 and b = −7.6 ± 1.7 × 10 −20 erg s −1 cm −2Å−1 . This ∼ 3% deviation from a 1-1 slope and low level of additive offset show that the calibrations were already in good agreement.

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 ′′ × 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, PSF k , where the PSF k 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 ′′ ± 0.2 ′′ 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.

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  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 ′′ 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, ∆SB/SB SDSS , is between −1.1% and 0.7% with a median value across all galaxies of −0.1% ± 4%.

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 gband 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 λ (x, y), with the g-band transmission curve, T g (λ), according to the following equation: The SDSS imaging is presented in units of nanomaggies or f ν = 3.631×10 −6 Jy. To compare with the KCWI data, with native units of flux density, f λ , we converted to flux density with the following expression (Tokunaga & Vacca 2005), where λ p is the pivot wavelength of the band-pass. The pivot wavelength for the g-band filter is λ p = 4702Å. The surface brightness for both KCWI and SDSS are measured inside 3 ′′ radius apertures. The mean and 1σ 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, ∆SB/SB 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 ± 7% which suggests that the KCWI calibration is in good agreement with SDSS.

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β 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α 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. 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 ∼ −0.1%.

Construction of Hβ Maps
Hβ 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β 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β map, pixels with weak or no Hβ 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β maps by integrating the continuum-subtracted spectra between 4856-4876 A. The final maps are suitable for H ii region identification using HIIphot.

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β (or Hα), then iteratively grows each "seed" and terminates only when the gradient of the Hβ (or Hα) 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, ∆ = 5 EM pc −1 or 2.43 × erg s −1 arcsec −2 pc −1 , as the recent PHANGS-MUSE work by Santoro et al. (2022).
HIIphot uses the PSF to convolve the input Hβ 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β 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 spu-rious small and pixelated regions or runaway growth. In total, HIIphot identifies ∼ 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 threefold decrease in strength of Hβ emission relative to Hα. In Figure 20 of Appendix C we show histograms of the Hβ 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.

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.
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, σ propagated , to the median standard deviation of the H ii region spectrum, σ measured , in the emission line free wavelength range 5400 − 5450Å. We measured an average ratio between the propagated and measured error of σ measured /σ propagated = 1.1 ± 0.2, implying that we are appropriately propagating the error. We also perform a similar comparison for KCWI, and find agreement, σ measured /σ propagated = 1.9 ± 0.5, between the . 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 (solidblack ). Across the sample, the offset between the KCWI and SDSS surface brightness is ∼ 1%. 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.

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β 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 mask-ing 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, [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β flux for KCWI and MUSE is −2.3 ± 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 χ 2 ≈ 1. Together the resulting fits for both the KCWI and MUSE have an average reduced χ 2 reduced ≈ 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 χ 2 reduced (Cappellari & Emsellem 2004;Cappellari 2017;Emsellem et al. 2022).

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α and Hβ emission, we construct a distribution of the MUSE Hα/Hβ ratios by sampling the error for each line. Next, using PyNeb, we calculate the extinction by comparing the average of the Hα/Hβ distribution to the theoretical value assuming the Case B recombination conditions n e =10 3 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×10 4 K and find that the standard deviation of A V for a fixed Balmer decrement is ∼ 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 ∼ 0.9 mag. We find a negligible difference when using KCWI Hβ in place of the MUSE Hβ 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/overestimated extinction may have on the measured electron temperature in Section 7.1.1.  For each emission line, we place a 10 ′′ × 10 ′′ aperture around each H ii region and measure the median DIG flux in pixels with S/N > 3 and Hα surface brightness Log 10 (SB Hα /[erg s −1 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.

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.  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  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, σ 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 σ 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 continuumsubtracted spectrum. After completing the N trials, we calculate the average integrated line flux, F avg , and the standard deviation of the measured fluxes, σ avg . If S/N > 3 we consider the auroral line detected. In Figure 28 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]λ4363 can be blended with [Fe II]λ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]λ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 We assess the ionized gas physical conditions-n 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.

Electron Density
Using PyNeb we calculate the electron density n e for each H ii region using the [S II]λλ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]λλ3726, 3729 doublet is also commonly used to estimate n e . The atomic levels responsible for [S II]λλ6716, 6731 and [O II]λλ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]λλ6716, 6731 and [O II]λλ3726, 3729 doublets are sensitive to densities 10 2 cm −3 < n e < 10 3.5 cm −3 . However, in the KCWI measurements, the [O II]λλ3726, 3729 doublet is unresolved, and therefore we do not use it to measure n e . With the exception of a handful of regions, the measured electron densities are in the low-density limit, n e < 100 cm −3 .

H II Region Electron Temperatures
The auroral lines allow the determination of the electron temperatures for the O + , O 2+ , N + , S + and S 2+ ions, or T e, [O II] , T e, [O III] , T e, [N II] , T e,[S II] , and T e,[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 e , to convert an auroral-tonebular ratio to temperature. The critical densities of the requisite lines, see Table 3, used to estimate T e, [O III] , T e,[N II] , and T e,[S III] are high enough such that the auroral-to-nebular lines ratios are insensitive to choice of n e < 10 4 cm −3 . The density sensitivity in the lines used for T e, [O II] and T e,[S II] begins at smaller densities n e ≈ 10 3 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.

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πR 2 n 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 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 as resolved studies show that H ii regions exhibit a range of non-spherical morphologies and non-uniform densities (Wood & Churchwell 1989 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 cospatial 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 CO , peak temperature, T peak and the  (2004) Tayal (2011)   velocity dispersion, σ v,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σ 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 σ v,CO span a similar range, up to Log 10 (σ v,CO /[km s −1 ])=1.5, with an average and standard deviations of Log 10 (σ v,CO /[km s −1 ])=0.88 ± 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.

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.

RESULTS
We present electron temperatures derived using the auroral-to-nebular line ratios from [O II] We construct T e -T e diagrams to compare any multi-ionization zone T 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 e , and ionization parameter, U . We also relate the temperatures to properties of the molecular gas and stellar populations. We observe in the T e -T e relations that T e, [O II] and T e,[S II] are systematically hotter than T e, [N II] . Based on photoionization models (e.g., Campbell et al. 1986;Garnett 1992) which is not compatible with the effects of overestimated extinction correction. Despite the above, our temperatures hierarchy could potentially be produced by underestimated extinction as this would lead to over-estimates of T e, [O II] and T e,[S II] while again leaving T e,[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]λλ7320, 7330 could be an additional systematic error, but is unlikely to be significant in our data because the [O II]λ7320/λ7330 ratio measured for our sample is 1.27 ± 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 start to be dependent on density around n e ≈ 10 3 cm −3 while [N II] is independent of density up to n e ≈ 10 4 cm −3 . We explore this dependence further in Section 7.2.

Intermediate Ionization Zone
The T e -T e relations between the low and intermediateionization 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 (Gar-nett 1992). We compare our observed T e -T 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 e, [S III] and T e, [O II] . The (BOND) models predict that the temperature T e,[S III] should be greater than T e,[O II] by a small constant offset across the full temperature range. However, our data suggest T e,[O II] rises faster than T e,[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 e, [S III] and T e,[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 ≈ 0.2×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 and T e,[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 e, [N II] . Regions with hotter T e, [N II] and T e,[S III] would be needed to to further differentiate between the models and empirical trends.

High Ionization Zone
We show T e -T 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 e -T e relations from CHAOS, Z21, and G92. Although we observe some H ii regions with T e -T e relations that agree with literature relations, we also observe numerous H ii regions with much higher T e, [O III] values of than what is predicted by models with the given low and intermediate ionization zone temperatures.
The large excess in T e, [O III] for regions with cooler low and intermediate ionization zone temperatures is unexpected. At the ∼ Solar metallicities of our sample, the relative flux [O III]λ4363 is expected to be < 10 −2 × Hβ (Berg et al. 2015). The small number of [O III]λ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 ∼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. 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, ∆(T ion,1 ,T ion,2 )=T ion,1 −T 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 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, ρ, and associated p-value, p. A correlation is judged to be significant if it exhibits p ≲ 10 −3 . The strength of the correlation is separated into the following regimes: 1 > |ρ| > 0.8 corresponds to a strong correlation, 0.8 > |ρ| > 0.4 corresponds to a moderate correlation and 0.4 > |ρ| > 0 identifies a weak or no correlation.
We do not report any significant correlations between ∆T e with any of the following properties: n e , integrated CO intensity, CO peak temperature, cluster mass, cluster age, and association age. Figures showing the ∆T 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. As expected, given that a majority of the n e are in the low-density limit, we do not observe any correlations between any ∆T e versus n e . However, recent studies have suggested that the temperatures obtained from the auroral-to-nebular lines ratios of [ 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 e, [O II] and T e,[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]λλ6731, 6716 and [O II]λλ3726, 3729 nebular line doublet ratios serve as density diagnostics for densities 10 2 cm −3 < n 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-tonebular 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 e < 10 6 cm −3 . For a uniform density H ii region, the n e returned from both of these diagnostics should be identical, as long as n 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  sus the region's T e,[N II] . We overlay the predicted trends of auroral-to-nebular line ratios calculated using n e = 10 2 cm −3 , 10 2.5 cm −3 , and 10 3 cm −3 . We see in Figure 11 that under the assumption that 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)   T e,[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).
In Figure 14, If we consider the similarity between our results and Yates et al. (2020)

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 31-33, we find no correlations between ∆T e with cluster mass, cluster age, or with association age.
We do observe a weak, positive correlation, between the most reliable ∆T e indicator without density inhomogeneity issues, ∆(T e, [N II] , T e,[S III] ) with the association mass in Figure 15. Here we see that T e,[N II] is cooler than T e,[S III] at lower association mass. As mass increases, we observe that the value of ∆(T e,[N II] , T e,[S III] ) decreases until T e,[N II] becomes larger than T e, [S III] . This correlation with association mass spans a small range in temperature differences ∼ [−0.2,0.2] ×10 4 K with association mass between 10 3 M ⊙ and 10 6 M ⊙ .
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 ∆(T e, [N II] , T e,[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, ρ = 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 ± 0.05, indicating that the error has negligible impact on the correlation.
It is possible that the correlation between ∆(T e,[N II] , T e,[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)  We perform no such cutoff in this study. The correlation between ∆(T e, [N II] , T e,[S III] ) and M Association correlation includes many regions with masses < 10 4 M ⊙ , 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 ∆(T e, [N II] , T e,[S III] ) warrants further investigation with a larger sample size on order to expand the dynamic range of association mass.

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 CO , shown in Figure 34 of Appendix E , and comparisons to T peak , shown in Figure 35 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 (σ v,CO ). We observe a moderate correlation between σ v,CO and ∆(T  The CO velocity dispersion can be enhanced by lowvelocity 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 pressureradiation driven H ii region expansion (Hill & Hollenbach 1978 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 lowionization species such as S + and O relative to Hβ and/or Hα (Dopita & Sutherland 1996;Allen et al. 2008). In Figure 18,  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, Xrays would provide high energy photons able to boost [O III]λ4363. However, lacking the high spatial resolution X-ray imaging of these H ii regions, exploration of X-ray contributions to [O III]λ4363 emission is beyond the scope of this paper.
Another tracer that may indicate presence of shocks is the He II λ4686/Hβ 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 (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]λ6731/λ6716 ratio is less sensitive to density than [Cl III], [Ar IV] and [Fe III] diagnostics when n e > 10 3 cm −3 (see Figure 2 in Méndez-Delgado et al. 2023b). If the equality T e,[O II] ≈T e,[S II] ≈T e,[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 re- gions 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 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 cm −3 have been measured using the [S III]λ18.7/33.5 µ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 e ∼ 1000 cm −3 assuming T e = T e, [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 e,[N II] a more reliable indicator of the low-ionization zone temperature due to its relative insensitivity to density.

T e,[N II] and T e,[S III] as Accurate Tracers of HII Region Temperatures
The low scatter and overlap with literature T e -T e trends between T e, [N II] and T e,[S III] in Figure 9 indicate that these temperatures may be optimal tracers of H ii region temperatures, more so than T e, [S II] , T e, [O II] and T e, [O III] .
The tight relationship between T e, [N II] and T e,[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 [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 e, [N II] and T e,[S III] generally follow trends predicted from photoionization models (Morisset 2009  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).

The High-Ionization Zone Temperature Excess
We found H ii regions with high T e,[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α, [O I]/Hα and the high-ionization ratio He II/Hβ. We found that the high T e,[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β 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 e, [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 e, [ ], 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×10 4 K to 7.2×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 low- In red we show the average temperature differences and H ii region association masses calculated in Log10(M⊙)=0.5 bins. We also show a fit and 1σ 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.
est and highest shock velocity models, Binette et al. (2012) found that mean doubly ionized fraction of oxygen, O 2+ /O , decreases while at the same time leaving the doubly ionized fraction of sulfur, S 2+ /S, unchanged. The imbalance the ionization fraction between oxygen and sulfur means that hotter, post-shock gas contributes proportionally more to the observed It has also been shown that shocks driven by pressureradiation 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 ∼ 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 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 exp =9.8 ± 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 ∼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α 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 exp = 10-40 km s −1 (see also Egorov et al. 2014Egorov et al. , 2017Cosens 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.

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.

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 ≳ 8.3 (Kreckel et al. 2019). Because electron temperature is anti-correlated to the metallicity of the gas, the [O III]λ4363 temperatures for these galaxies are expected to be low. For [O III]λ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 lowermetallicity sample may be compatible with photoionization models and T e -T e relations. Finally, Rola & Pelat (1994)  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: • In a follow-up paper, Rickards Vaught et al. (in prep), We will test temperature recommendations for measuring "direct" metallicities using our full set of measured auroral line temperatures. This work, along with the PHANGS-MUSE survey, demonstrate the power of integral field spectrographs on 10m class telescopes for measuring faint auroral emission lines from large samples of H ii regions in nearby galaxies. Future efforts with deeper observations or expanded samples will be critical for further elucidating the temperature and ionization behavior of these regions, particularly as [O III]λ4363 and other faint lines are now being routinely detected in galaxies at high redshift with JWST and used in metallicity determinations.

ACKNOWLEDGEMENTS
The data presented herein were obtained at the W. M. Keck Observatory, which is operated as a scientific partnership among the California Institute of Technology, the University of California and the National Aeronautics and Space Administration. The Observatory was made possible by the generous financial support of the W. M. Keck Foundation. The authors wish to recognize and acknowledge the very significant cultural role and reverence that the summit of Maunakea has always had within the indigenous Hawaiian community. We are most fortunate to have the opportunity to conduct observations from this mountain. We also wish to thank all the Keck Observatory staff, including Gregg Doppman, Luca Rizzi, Sheery Yeh, and Rosalie McGurk for observational support.
This research is also based on observations collected at the European Southern Observatory under ESO programmes 094.C-0623 (PI: Kreckel) This research made use of Montage. It is funded by the National Science Foundation under Grant Number ACI-1440620, and was previously funded by the National Aeronautics and Space Administration's Earth Science Technology Office, Computation Technologies Project, under Cooperative Agreement Number NCC5-626 between NASA and the California Institute of Technology. This research has made use of NASA's Astrophysics Data System. This research made use of Astropy, a community-developed core Python package for Astronomy (Astropy Collaboration et al. 2013), as well as substantial use of the nebular diagnostics toolkit Pyneb (Luridiana et al. 2015).
The distances in Table 1 were compiled by Anand et al. (2021) and are based on Freedman et al. (2001); Nugent et al. (2006) Table 7 and 8. Additionally, we describe below the measurement of seeing from the set of standard stars.
To measure seeing, we fit a 2D Gaussian to each standard star observation. We find an average FWHM of ∼ 1.2 ′′ , however, at some points during the nights of 10-17-2018, 03-27-2019, and 03-28-2019 the seeing was poorer with values between 1.6 ′′ -2 ′′ . Aside from these portions of the nights, the seeing was stable near the average FWHM. The FWHM measurements for each standard star observation are summarized in Table 8. Table 7. Summary of KCWI Observations. This table reports the the galaxy and the field observed (Field). Note the repeated fields are offset by 1/2 slice width, with the exception of all the fields in NGC3627. The table also reports: the original coordinates of the field center in right ascension and declination (RA, DEC), the applied astrometric correction in right ascension and declination (RAcorr, DECcorr) determined from the image registration, the UTC start-time and total integration the exposure in seconds (Date, Exposure) as well as the airmass (Airmass), position angle (Angle) and the angular FWHM of the pointing.     a The field numbering is from internal lists of potential pointings which were not performed in numerical order, therefore the field number does not reflect the sequence of observations. b The astrometric correction is described in Section 3.2.
c The measured angular FWHM for the pointing, see Section 3.4.  a The name of each standard star as displayed in the KDRP list of standards.
b The angular FWHM was measured from a 2D Gaussian fit to a white-light image of the standard star.

B. KCWI LINE SPREAD FUNCTION
The KCWI line spread function in the large slicer is non-Gaussian. Deviations from a Gaussian profile have been observed in KCWI data before. For example, using a combination of the medium slicer and grating van Dokkum et al. (2019). Following van der Marel & Franx (1993), in order to parameterize the degree of deviation from a Gaussian profile, we fit the spectra from a pipeline reduced arc lamp exposure a Gauss-Hermite function of the form: with X = (λ − λ 0 )/σ, amplitude (γ) and spectral width (σ). The anti-symmetric and symmetric deviations from pure Gaussian profiles are captured by the constants, h 3 and h 4 . A Gaussian function is recovered by setting h 3 = h 4 = 0.
In Figure 19 we show histograms of the fitted values from h 3 and h 4 after fitting Eq. B1 to approximately 20 isolated arc lamp emission lines. We find that the emission lines in the arc lamp images are consistently flat topped, with an average value of h 4 = −0.14 ± 0.01 across all spatial pixels. The degree that the line profile is non-Gaussian due to asymmetric deviations is small compared to symmetric deviations with an average value of h 3 = −0.007 ± 0.01. The authors state that the root cause of the deviation from a Gaussian profile is due to the slit width-limited resolution of the medium slicer (see, Casini & de Wijn 2014). We posit that a similar limitation for the large slicer is responsible for the deviations measured in this work.

C. COMPARISON OF H ii REGION IDENTIFICATION BETWEEN KCWI AND MUSE
We identify potential H ii regions with Hβ emission maps constructed from the KCWI galaxy mosaics using HIIPhot. H ii regions, for the same galaxies, in the PHANGS Nebular Catalog were identified with HIIPhot and MUSE Hα maps (Kreckel et al. 2019;Groves et al. 2023). Given the higher resolution of the MUSE imaging, as well as the three-fold brightness increase of Hα relative to Hβ, we expect our H ii region catalog to be less sensitive to the faintest and smallest H ii regions.
We also observe in Figure 20 that regions with radii less than the KCWI angular resolution are missed KCWI-Hβ region sample. This can be seen clearly in Figures 21-27, where we compare the boundaries of the MUSE-Hα and KCWI-Hβ regions. In these figures, we see that many of the missed regions are small, and unresolved in the KCWI Hβ map. Also owing to the larger number of detections is that many of the larger KCWI-Hβ regions are resolved into smaller structures in the MUSE-Hα regions.       In this section we show example auroral line fits. In Figure 28 we show fits to high S/N auroral lines from an H ii region in NGC 5068. We also include in Figure 28 annotations describing the the standard deviation of the fit residuals, σ res , the S/N for a single emission line (or in the instance of simultaneous double line fits the S/N of the red and purple Gaussian fit: S/N r and S/N p ), the continuum noise, σ cont and reduced χ 2 . In this particular example, the auroral line in these fits are isolated from contaminating sky lines or nearby emission, especially in the case of [O III] Figure 29. . Auroral Line fits for an H ii region in NGC 5068. We show in each panel a summary of auroral line fits for a single H ii region in NGC 5068. The top frame in each panel shows the data, (solid-black ), for the fitting and continuum wavelength ranges. The (shaded-grey) show the 1σ ranges of the trial spectra whilst the (shaded red and purple for the double Gaussian fits), show the 1σ ranges of the fitted models. The bottom frame in each panels shows a histogram of the residuals. We also print text summarizing the S/N and average reduced χ 2 of the fits as well as the 1σ of the residuals and value of σcont.   Low-Intermediate-High Ionization Zones Figure 33. Electron temperature differences compared to the stellar association age. The order of the panels follow those in Figure 30.    Low-Intermediate-High Ionization Zones Figure 35. Electron temperature differences compared to the CO peak temperature, T peak . The order of the panels follow those in Figure 30.