The impact of pre-supernova feedback and its dependence on environment

Integral field units enable resolved studies of a large number of star-forming regions across entire nearby galaxies, providing insight on the conversion of gas into stars and the feedback from the emerging stellar populations over unprecedented dynamic ranges in terms of spatial scale, star-forming region properties, and environments. We use the VLT/MUSE legacy data set covering the central $35$ arcmin$^{2}$ (${\sim}12$ kpc$^{2}$) of the nearby galaxy NGC 300 to quantify the effect of stellar feedback as a function of the local galactic environment. We extract spectra from emission line regions identified within dendrograms, combine emission line ratios and line widths to distinguish between HII regions, planetary nebulae, and supernova remnants, and compute their ionised gas properties, gas-phase oxygen abundances, and feedback-related pressure terms. For the HII regions, we find that the direct radiation pressure ($P_\mathrm{dir}$) and the pressure of the ionised gas ($P_{HII}$) weakly increase towards larger galactocentric radii, i.e. along the galaxy's (negative) abundance and (positive) extinction gradients. While the increase of $P_{HII}$ with galactocentric radius is likely due to higher photon fluxes from lower-metallicity stellar populations, we find that the increase of $P_\mathrm{dir}$ is likely driven by the combination of higher photon fluxes and enhanced dust content at larger galactocentric radii. In light of the above, we investigate the effect of increased pre-supernova feedback at larger galactocentric distances (lower metallicities and increased dust mass surface density) on the ISM, finding that supernovae at lower metallicities expand into lower-density environments, thereby enhancing the impact of supernova feedback.

Many numerical studies have shown that including stellar feedback in galaxy simulations is required to recover global observational properties such as star formation rates and star formation efficiencies (e.g. Agertz et al. 2013;Hopkins et al. 2014;Crain et al. 2015;Fujimoto et al. 2019;Keller & Kruĳssen 2020), to replicate well known scaling relations like the Kennicutt-Schmidt relation (Schmidt 1959;Kennicutt 1998), and to reconcile the long global depletion time (∼1 Gyr; e.g. Bigiel et al. 2008;Leroy et al. 2008) with the short gravitational collapse time-scale (∼10 Myr; e.g. Heyer & Dame 2015;Utomo et al. 2018;Schruba et al. 2019). Different numerical studies implement different (combinations of) feedback mechanisms, and where multiple ones are included, efforts go towards understanding the relative importance of individual feedback mechanisms. For example, supernovae (SNe) have long thought to be the dominating source of feedback, but both simulations (e.g. Haid et al. 2018;Lucas et al. 2020;Semenov et al. 2021;Keller et al. 2021) and observations (e.g. Kruĳssen et al. 2019;Chevance et al. 2020a) now clearly show that early stellar feedback (i.e. radiation pressure, ionisation, and stellar winds) play a major role in regulating the impact of SN feedback by altering the conditions of the interstellar medium (ISM) prior to the first SN events. These results indicate that not only is the environment being affected by stellar feedback, but environmental properties in turn set the effectiveness of feedback. While the impact of feedback on the environment is an established area of active research, the impact of the environment on feedback is now starting to be explored from the observational perspective (e.g. Lopez et al. 2011Lopez et al. , 2014Chevance et al. 2016; Barnes et al. 2020;Olivier et al. 2021). Outstanding questions include: What are the dominant feedback mechanisms from massive stars as a function of their stellar properties (e.g. mass, chemical composition, rotation rate, binarity, evolutionary phase)? How does feedback change with environment (e.g. metallicity, ambient gas density, location within a galaxy)? How does our knowledge of feedback change with physical scale, from small (clouds) to large (galaxies) scales?
Over the past decade, the increasing availability of large fieldof-view, large wavelength coverage, and medium spectral resolution integral field unit (IFU) instruments has enabled the simultaneous study of the feedback-driving stellar populations and the feedbackaffected matter. For example, this has led to optical studies of the stellar and ionised gas properties and kinematics of entire spatiallyresolved star-forming regions in the Milky Way (e.g. McLeod et al. 2015;Weilbacher et al. 2015;McLeod et al. 2016;Flagey et al. 2020), the Magellanic Clouds (e.g. Castro et al. 2018;McLeod et al. 2019b), and nearby galaxies (e.g. Monreal-Ibero et al. 2011, 2012Westmoquette et al. 2013;McLeod et al. 2020). While these regions are observationally convenient for detailed multi-wavelength studies of feedback on small scales or in select galactic hosts, they are not broadly representative of star formation and feedback over all parameter space and do not consider the effects of feedback in the larger context of their galactic hosts. The need for large region samples spanning a vast parameter space and being spatially-matched with available multi-wavelength ISM observations has produced large nearby galaxy IFU surveys such as SIGNALS 1 (Rousseau-Nepton et al. 2019) and PHANGS 2 -MUSE (Schinnerer et al. 2019;1 Star formation, ionised Gas and Nebular Abundances Legacy Survey with SITELLE. 2 Physics at High Angular Resolution in Nearby Galaxies. Emsellem et al. 2021, for early results see Kreckel et al. 2019;Pessa et al. 2021). For tens of thousands of regions, these surveys deliver spatially-resolved (for the nearest systems up to a few Mpc) and integrated (beyond a few Mpc) stellar and ionised gas properties, enabling environmental studies of stellar feedback in orders of magnitude more star-forming regions than previously possible. These are also well-matched with recent advancements made in computational galaxy evolution models.
Here, we study the environmental dependence of ionised gas properties and feedback-related pressure terms in the nearby galaxy NGC 300, and study their implications in the context of early pre-SN and subsequent SN feedback. This galaxy has been the focus of two recent feedback-related studies upon which this present paper builds. (1) Kruĳssen et al. (2019) use a novel statistical method based on the spatial decorrelation between young stars and molecular gas to infer feedback-related quantities across the galactic disc of NGC 300 (i.e. molecular cloud lifetimes, feedback time-scales, outflow velocities, mass-loading factors, and star formation efficiencies). They show that star formation in NGC 300 is rapid and inefficient, with giant molecular clouds (GMCs) having integrated star formation efficiencies of only 2−3 per cent, but being dispersed by feedback from massive stars within 1.5 ± 0.2 Myr. (2) McLeod et al. (2020, henceforth referred to as M20) study 5 H II regions in NGC 300 in a spatially resolved manner and find that their expansion is governed by the pressure of the ionised gas and by the winds from the massive stars within them.
In this paper, we use a legacy value data set covering the inner star-forming disc of NGC 300 taken with the VLT/MUSE instrument (Bacon et al. 2010), consisting of a contiguous 7 × 5 mosaic (∼4 kpc × 3 kpc) and covering the galaxy out to galactocentric radii of about 0.45 25 (with 25 ∼ 5.33 kpc being the optical radius; Paturel et al. 2003). NGC 300 is an ideal target to study stellar feedback: it is the closest ( ∼ 2 Mpc; Dalcanton et al. 2009), non-interacting, star-forming disc galaxy that can be mapped at the necessary spatial resolution (i.e. 1 , which corresponds to ∼10 pc, resolving individual star-forming regions and supernova remnants). The large spatial coverage (to cover most of the star-forming disc) available not only in the optical with MUSE but throughout the electromagnetic spectrum (e.g. Helou et al. 2004;Westmeier et al. 2011;Riener et al. 2018;Kruĳssen et al. 2019;Schruba et al. 2019) makes NGC 300 the ideal target for simultaneous resolved feedback, stellar population, and ISM studies. Closer galaxies like the Magellanic Clouds, M31 or M33 do not allow similar large-scale multi-wavelength mapping due to their large angular sizes, while more distant galaxies (beyond a few Mpc) do not allow multi-wavelength studies with similar spatial resolution across the optical, infrared, mm/sub-mm, and radio. NGC 300 perfectly bridges between ∼100 pc resolution IFU surveys of nearby galaxies like PHANGS, and upcoming (sub-) pc scale resolution IFU surveys of the Milky Way and the Magellanic Clouds like SDSS-V/LVM (Kollmeier et al. 2017). Further, NGC 300 offers a favorable inclination of ∼40 • (Puche et al. 1990), it is actively forming stars at a rate between ∼0.08 and ∼0.30 M yr −1 (Kang et al. 2016, and references therein), and has a well-studied population of H II regions (e.g. Deharveng et al. 1988;Bresolin et al. 2009;Faesi et al. 2014), planetary nebulae (PNe; e.g. Soffner et al. 1996;Peña et al. 2012;Stasińska et al. 2013), and supernova remnants (SNRs; e.g. Blair & Long 1997;Millar et al. 2012;Vučetić et al. 2015).
This paper is organised as follows. After a brief overview of the VLT/MUSE observations in Section 2, in Section 3, we describe the methods used to identify and classify emission line regions. In Section 4, we compute ionised gas properties and feedback-related pressure terms for the detected H II regions and discuss environmental dependencies. In Section 5, we analyse the environment in which the covered SN events occurred. Finally, we summarise our findings and conclusions in Section 6.

OBSERVATIONS
This work is based on the VLT/MUSE data set of NGC 300 first presented in M20. The data were taken in the nominal wavelength range of the MUSE instrument (∼4750−9350 Å) and using its widefield mode (∼1 × 1 per pointing), as part of the observing program 098.B-0193(A) (PI McLeod). As opposed to M20, where only 2 of the NGC 300 MUSE data cubes are analysed, here we exploit the full coverage of the in total 35 individual mosaic pointings which cover a 7 × 5 contiguous mosaic of the central star-forming disc of NGC 300. The data were taken prior to the availability of the Adaptive Optics system for MUSE's large field-of-view, and seeing-limited angular resolutions in a range of about 0. 45−1. 3 were achieved (see Table 1). Each individual mosaic pointing was observed three times in a 90 • -rotation dither pattern with an exposure time of 900 seconds per rotation. The full mosaic is shown in Fig. 1, which consists of a three-colour composite of the 35 pointings with three emission lines tracing the ionised gas (red = [S II] 6717, green = H , blue = [O III] 5007), overlayed on an optical ESO-DSS image for reference. Observational details of the individual pointings are given in Table 1.
As described in M20, we proceed in reducing the data with the MUSE pipeline (Weilbacher et al. 2012) in the environment with the standard static calibration files. For each specific observing block we use the available calibration files from the ESO archive, and subtract the sky lines according to Zeidler et al. (2019) (for brevity we do not describe sky subtraction details here and refer the interested reader to M20). The flux calibration (also performed with the MUSE pipeline) is done for each observing block using the matched nightly standard star observations. The three exposures of each individual pointing are combined into single cubes with the built-in exposure combination recipes of the MUSE pipeline.
The extended nature of the types of objects analysed here (i.e. mainly H II regions and SNRs) necessarily means that some of these significantly overlap between different cubes, making it imperative to mosaic these cubes such that an integrated spectrum for the regions that overlap between two or more individual pointings can be extracted. Two issues arise prior to producing a mosaic of the cubes. First, combining 35 individual MUSE cubes, each several GBs in size, would result in a single, giant cube of unreasonably large (and thus unwieldy) file size. Second, the individual pointings were observed during different nights and at varying observing conditions, leading to relative flux offsets across the field-of-view. To overcome these two problems we proceed in the following way: • The 35 cubes are first resampled to a common wavelength grid (to overcome slight wavelength offsets between cubes) and divided into sub-cubes spanning 500 wavelength elements each (i.e. 625 Å; the full wavelength range of each MUSE cube is 4750−9350 Å and the sampling is 1.25 Å).
• Individual 2D slices from each sub-cube are then mosaicked using the A (Price-Whelan et al. 2018) Montage wrapper with a background match.
• The 2D mosaics are then recombined into data cubes which now cover the entire field-of-view and span across a manageable wavelength range. Spectra of regions of interest can then be extracted from each one of the large cubes spanning about 625 Å each, and combined to cover the entire wavelength range. To assess the performance of combining the cubes as described, we compare fluxes obtained from the mosaicked MUSE data to those reported in previous spectroscopic studies of regions in NGC 300While the comparison with fluxes from photometric studies (e.g., Faesi et al. 2014) is feasible, it requires carefully reproducing the filter parameters of the used instrument. The H II regions used for this comparison are spread across the MUSE FOV to compare fluxes across different pointings. We compare MUSE fluxes with those reported in Toribio San Cipriano et al. (2016), who use VLT/UVES data to derive abundances of 7 H II regions in NGC 300, 3 of which overlap with the multi-night MUSE data (specifically, their R20, R23, and R76a are in MUSE fields H8, H12, and L4, respectively). Crucially, for each of these regions, Toribio San Cipriano et al. report the central coordinates and area covered by their observations (their Table 1), and we extract integrated spectra from the MUSE data accordingly. The comparison, summarised in Table 2, shows excellent agreement (within errors) between the MUSE and UVES fluxes, with the exception of the H flux of R23, which is likely due to a repetition typo in Table 1 Bresolin et al. (2009) or Stasińska et al. (2013), as both of these studies use VLT/FORS2 data but do not specify the slit lengths adopted for individual regions, thus hindering a direct comparison. Further, for the purpose of the subsequent analysis, we note that when making the mosaic no PSF matching was performed, as the region size constraints described below ensure that all analysed objects are resolved regardless of the seeing achieved in a particular field, and because all analyses are performed on integrated spectra (and thus do not retain PSF information).
Emission line maps (such as those shown in Fig. 1) are obtained by collapsing the cubes over ±3 Å (about ±140 km s −1 at H ) around the lines of interest. The analysis described in the next section is based on a continuum-subtracted H map. For this, we first produce a continuum map by summing (i.e., collapsing) over a wavelength range equal in width to that used for the H line but covering a portion of nearby continuum (i.e. free of emission/ absorption features, centered on 6540 Å), and then subtract this from the H map. The analyses presented here are not sensitive to small (∼1 −2 ) WCS shifts relative to, e.g. HST coordinates often observed in MUSE observations, and we therefore do not perform a WCS correction here. However, this is necessary when wanting to combine these with other data, and we note that data releases for this legacy MUSE data set of NGC 300 (which is part of several ongoing student projects) are scheduled to commence late 2022, in the form of fully reduced and WCS-corrected cubes, catalogs, and spectra.

REGION IDENTIFICATION AND CLASSIFICATION
While the objects of interest in this study are H II regions and SNRs, the emission line regions as traced by the ionised gas also include sources of different nature (e.g. PNe, emission line stars, microquasars, ultra-luminous X-ray sources). It is therefore necessary to isolate the regions of interest from the general population of emission line sources. As already mentioned in Section 1, NGC 300 is a well-studied galaxy and substantial catalogs of emission line regions have been compiled by previous studies (e.g. Deharveng et al. 1988;Stasińska et al. 2013;Vučetić et al. 2015). This is the first time, however, that a large-scale IFU data set exists for this galaxy, meaning that we can now perform spectro-photometric studies of these regions rather than relying on either only photometry, or targeted (slit) spectroscopy of selected regions. Given the rapid rise and wealth of already existing IFU observations of nearby galaxies, we use this MUSE data set of NGC 300 to explore new and efficient empirical identification and classification methods which can be readily adapted and applied to other optical IFU data of similar spatial resolution (e.g. the closest galaxies of the SIGNALS survey, Rousseau-Nepton et al. 2019, or NGC 7793, Della Bruna et al. 2020. The MUSE NGC 300 data set gives access to simultaneous, spatially-resolved, photometric and spectroscopic information of >100 regions that are bright in the main nebular emission lines. We can therefore analyse population trends, such as radial abundance gradients, with improved number statistics. For example, the largest spectroscopic study of H II region abundances in NGC 300 is that of Bresolin et al. (2009), who obtained deep spectra for 28 H II regions. While our data are not as deep (e.g. we do not obtain sufficient signal-to-noise on faint auroral lines needed to determine electron temperatures and temperature-based ionic and elemental abundances) and do not extend to as large galactocentric radii, we are able to spectroscopically analyse a factor ∼7 more H II regions. In what follows, we describe how emission line regions are identified and how they are separated into the three main categories, these being H II regions, SNRs, and PNe.

Emission line region identification
In a first step, we proceed in identifying and isolating emission line regions across the MUSE mosaic footprint, regardless of their nature. The input for this step is a continuum-subtracted integrated H map, which reliably traces ionised gas in H II regions and SNRs in the optical regime. PNe on the other hand are known to sometimes have very weak Balmer emission (Zhang et al. 2004), and to identify all of them an additional tracer of highly ionised gas (e.g. [O III] 5007) should be included. We will perform a detailed analysis of the PNe present in the MUSE data in a forthcoming study, and note that the goal of this paper is not that of obtaining a complete PNe census, but rather correctly classifying those PNe that are picked up by the emission line region identification algorithm and excluding them from the H II region and SNR catalogs.
A variety of different methods to identify emission line structures are found in the literature, examples include HIIPHOT (Thilker et al. 2000) and Clumpfind (Williams et al. 1994), the latter having been applied to MUSE data to identify H II regions in NGC 628 (Kreckel et al. 2016). Here, we proceed in a similar fashion to Della Bruna et al. (2020), who used the Python package -3 to identify bright regions in a MUSE H map of NGC 7793 ( ∼ 3.4 Mpc), but we add an extra spectral clustering step to the classical dendrogram hierarchy 4 of 'trunks', 'branches', and 'leaves' which, as described below, is not ideal at the high spatial resolution achieved with MUSE in NGC 300. Hence, to identify emission line regions in the H map, we use the following simple two-step approach: (i) Regions are first broadly isolated by computing a dendrogram of the H map, exploiting the fact that dendrograms divide the emission in a 2D map into hierarchical structures based on user-defined minimum-flux (with respect to a background) and size thresholds.
(ii) The dendrogram structures are then fed into the algorithm (Colombo et al. 2015) which groups these into coherent and relevant regions based on a spectral clustering paradigm.
Dendrograms are heavily dependent on so-called user-defined pruning parameters. On opposite extremes, different pruning parameter choices can lead to either very few large regions encompassing what clearly are individual structures or unrealistically many small fragments. Here, initial pruning parameters are set such that the catalog resulting from the two-step approach approximately contains the expected number of regions in the MUSE footprint based on both visual inspection and on the literature, i.e. ∼83 H II regions (Deharveng et al. 1988), ∼20 PNe (Stasińska et al. 2013), and ∼12 SNR (Millar et al. 2012). The detection flux threshold is set to three times the standard deviation of the H flux map, i.e. about 5 × 10 −18 erg s −1 cm −2 or about twice the mean H flux of the map. This comfortably includes the faintest structures identifiable by eye and corresponds to an emission measure of about 60 pc cm −6 , which is somewhat lower than the cutoff value of 82 pc cm −6 used by Della Bruna et al. (2020) in NGC 7793 and closer to 50 pc cm −6 set by Hoopes et al. (1996) in NGC 300. We tested different thresholds, finding that higher values do not return the contours of by-eye identified individual regions, whereas lower thresholds lead to large structures encompassing multiple individual regions. The minimum number of pixels is set to 25, which, for circular apertures, corresponds to diameters of ∼1 , i.e. approximately the highest seeinglimited resolution achieved by our observations. With these settings, the resulting dendrogram contains 794 structures, which are then passed to the clustering algorithm.
The package was originally developed to identify giant molecular cloud structures. At its core is an unsupervised pattern recognition algorithm which, in very general terms, groups together pixels in an image which are considered to be similar to each other. As noted in Colombo et al. 2015, this approach to structure identification overcomes the problem introduced by high spatial resolution observations which causes tools such as dendrograms to overestimate the number of regions (i.e. over-divide the input image into too many small structures). We run on the dendrogram obtained with the pruning parameters described above, and set the clustering to be performed based on H flux and to return isolated leaves as well as grouped ones as independent structures. The latter ensures that unresolved emission line regions, e.g. PNe, are included in our structure identification.
With the above described pruning and clustering settings, we obtain an initial catalog of 204 emission line regions and we extract integrated spectra for all of these based on the identified region contours (see Fig. 2). While some of these regions encompass what are likely multiple regions, these are the minority and only marginally affect the H II region analysis described below, given that our main interest lies in radial trends. As no method is 100 per cent accurate in recovering individual structures, in large field-of-view data sets like the one considered here, the vastly greater efficiency of our semi-automated region identification method is certainly preferred over a subjective by-eye identification. Upon testing different pruning parameters within a sensible range (i.e., as not to return unrealistically many small or just a few large regions), the results from the analyses are unchanged. From an initial visual inspection of the resulting integrated spectra we eliminate 16 objects which are either of insufficient quality (this is particularly true for regions in cube L15 which is heavily contaminated by a bright foreground star), or correspond to known stellar sources such as emission line stars (e.g. Roth et al. 2018), WR stars (e.g. Schild et al. 2003), as well as PNe (e.g. Stasińska et al. 2013) that have very low signalto-noise MUSE spectra. The list of eliminated spectra also includes the ultra-luminous X-ray pulsar ULX-1 (Vasilopoulos et al. 2019;Binder et al. 2020) as well as a known microquasar (McLeod et al. 2019a;Urquhart et al. 2019), both of which were picked up in the emission line region identification step due to their interaction with the surrounding ISM (i.e. shock-excited gas). The remaining 188 spectra are passed to the Gaussian fitting routine and subsequent classification scheme described in the next section.

Emission line fitting and region classification
Before fitting the emission lines, the integrated spectra are corrected for Balmer absorption, caused by the unresolved stellar background, using the Python implementation of PXF (Cappellari 2017) together with spectral templates from the MILES library (Sánchez-Blázquez et al. 2006) (as well as a standard list of emission lines to mask). We perform the PXF fit in the range 4750−7200 Å, thus, excluding redder wavelengths from the fit that contain contaminating residual sky emission. We adopt a MUSE LSF parametrization as in Guérou et al. (2017) and assume initial guesses of 146 km s −1 for the systemic velocity (from eq. 8 in Cappellari 2017 and a redshift of ∼ 0.000487 for NGC 300) and 20 km s −1 for the velocity dispersion (the latter based on a preliminary inspection of the spectra). While the main output of PXF fitting consists of the kinematics of the unresolved stellar population, the purpose here is purely to subtract the best-fitting spectral template from the observed integrated emission line spectra and, thus, to obtain a continuum-subtracted and absorption-corrected nebular spectrum for each region. A spatially-resolved study of the stellar kinematics across the entire MUSE mosaic will be presented in a forthcoming publication. We fit the emission lines in each spectrum with the Python package P S K (Ginsburg & Mirocha 2011) assuming single component Gaussians, and correct the obtained line fluxes for extinction using the PyNeb package (Luridiana et al. 2015) based on the Balmer decrement (with an intrinsic H /H ratio of 2.86), assuming = 3.1 and a Galactic extinction curve (Cardelli et al. 1989). Uncertainties on the ionised gas properties obtained from the Gaussian fits described in the following sections are derived by propagating the errors on the best-fit parameters from P S K . As most of the analyses rely on emission line ratios, intrinsic uncertainties are minimised. A further note on the contribution of diffuse ionised gas (DIG) to the integrated spectra. While this paper is not aimed at studying the DIG, it is important to assess whether (and to what extent) the integrated region spectra are contaminated by DIG emission, i.e how much DIG is likely included in the contours we extract spectra from. To this end, we use the region contours (see Fig. 2) to produce a rough DIG map by simply masking all the pixels within the region contours in the continuum-subtracted H map. We then crudely estimate the amount of DIG in the covered portion of the galaxy by summing the pixel values of the DIG map dividing by the corresponding value of the summed H map. We obtain a DIG fraction of ∼ 47±2%, which is in good agreement with Hoopes et al. (1996) who find a DIG fraction of 53±5% for their emission measure threshold of 50 pc cm −6 . Together with the completeness tests described in Section 3.2.3, we therefore conclude that the DIG contribution to the extracted region spectra is negligible.
The catalog of 188 regions obtained as described in Section 3.1 mainly consists of H II regions, SNRs, and PNe. For the H II region and SNR analyses described in the following sections of this paper, the catalog objects therefore need to be classified by type. While several catalogs of H II regions, SNRs, and PNe already exist in the literature, we intentionally do not cross-match our initial emission line region catalog with known sources. This is for two reasons, the first one being that with ongoing large surveys delivering IFU coverage of entire nearby galaxies (e.g. Rousseau-Nepton  Emsellem et al. 2021), it is in the community's interest to use well-studied galaxies like NGC 300 to find empirical classification methods specifically tailored to the capabilities of the instrument, which can then be used for less well-studied targets. The second reason is that determining the center, size, and boundaries of H II regions at spatial resolutions of roughly 7−10 pc (i.e. what is achieved with MUSE in NGC 300) is very subjective and very much dependent on the method that is being used, particularly where regions are in close vicinity or even overlap, meaning that H II region catalogs of the same galaxy but from different papers can vary significantly. This is particularly noticeable for H II region catalogs that have been defined by eye versus those that result from more sophisticated approaches relying on structure identification algorithms. We do perform a literature cross-match for SNRs and PNe in a second step, which serves the purpose of validating our classification method.
In what follows, we describe the classification scheme used to disentangle between three main groups of objects, i.e. H II regions, SNRs, and PNe. For illustration purposes, Fig. 3 shows normalised spectra representative of an H II region, a SNR, and a PN classified as per the below. The spectra in this figure are cropped to wavelength ranges covering relevant emission lines, i.e. the H and

Supernova remnants
The most commonly used line ratio diagnostic to identify SNRs in external galaxies is [S II]/H , as it gives a relative measure of the ionisation stages of sulphur within a given region. In H II regions, where the gas is mostly photoionised, a larger fraction of sulphur is expected to be in the higher excitation state (S ++ ), and [S II]/H ratios are typically around 0.1 (Long et al. 2018). In SNRs, where shocks contribute to populating S + , [S II]/H values are expected to be 0.4 or higher (Allen et al. 2008). Blair & Long (1997) Long et al. 2018), are expected to be resolved in our observations. Older SNRs with line widths <100 km s −1 might therefore be misclassified, however, these also are typically very faint and would likely fall below our detection threshold.
As is shown in Fig Millar et al. (2012) sources that lie within the MUSE field, and find that 7 out of the 8 sources identified here indeed correspond to previously known SNRs. The unmatched source (id #538, see Table B1) is the only one of the 8 to be unresolved in our data, and we therefore exclude it from further analyses. We do however assign this a SNR candidate flag in our catalog. Thus, 7 of the 11 Millar et al. sources  Of the remaining 4 Millar et al. sources that we do not recover (S08, S09, S19, and S22), S08 and S09 do not show line broadening and are consistent with being H II regions in our sample (based on their line ratios, line widths, and upon visual inspection), S22 is a faint structure that lies below our detection threshold (and upon further inspection does not show line ratios and line widths consistent with a SNR), and S19 is grouped into a contour with an adjacent H II region (we will discuss the implication of this in Section 3.2.3).
The above method to distinguish SNR from other emission line regions is thus robust, but we recommend that, if the purpose is that of identifying unknown SNRs rather than removing them from an H II region sample, the detection threshold in the dendrogram emission line region identification step is lowered to also include potentially low-brightness SNRs. We also recommend using a [S II] map in addition to the H map when specifically identifying SNRs.

Planetary nebulae
With the SNRs being removed from the emission line catalog according to the empirical separation described in Section 3.2.1, we now seek to disentangle H II regions from PNe. Countless pairs of line ratios and other empirical relations have been used in the literature, e.g. the traditional BPT diagnostic plot (Baldwin et al. 1981), the H /[S II] versus H /[N II] diagnostic (Riesgo & López 2006), or more recently, with MUSE data of NGC 628, empirical narrow-band criteria (Kreckel et al. 2017).
Here, we exploit the fact that PNe typically exhibit large [O III]/H ratios (e.g. Baldwin et al. 1981) and, due to their hot central stars, large degrees of ionisation, together with the fact that at the distance of NGC 300 PNe are spatially unresolved and, therefore, appear as point sources (with PNe typically having radii 1 pc, Jacob et al. 2013). This is shown in the lower panel of Fig. 5 solved). This selects the sources residing in the teal box in Fig. 5.
Several more pairs of quantities to disentangle H II regions from PNe were explored in NGC 300 by Stasińska et al. (2013), who found a clear sequence in the H line luminosity, (H ), versus the ionised gas mass, ion , ranging from PNe to compact H II regions and to evolved H II regions. While we clearly recover this sequence, shown in the upper panel of Fig. 5 (data points are colour-coded by region radius), the large uncertainties in estimating the electron density, e , and the rather coarse assumptions on stellar initial mass function (IMF) sampling when converting the H luminosity to an ionising flux, (H ) (see Section 4.2), do not allow reliable estimates of ion , which we compute as in M20, with the proton mass p and the Case B coefficient B (corresponding to optically thick nebulae), and assuming negligible dust absorption. Electron densities, e , are computed via from the [S II] 6717/[S II] 6731 line ratio assuming a temperature of 10 4 K 5 , and while H II regions in our sample have measured electron density uncertainties, e , on the order of ∼20 per cent, the PNe have e ranging from about 25 per cent to over 200 per cent. Given the consequently large uncertainties on ion and given that there is no clear separation between PNe and H II regions, we do not use this parameter pair to disentangle the two types of regions.
To confirm the validity of our PNe selection criteria, we crossmatch the resulting list of 13 PNe candidates with the 18 PNe from the Stasińska et al. (2013) catalog falling within the MUSE field-ofview, confirming all but one of our PNe candidates, i.e. 12 out of the expected 18 PNe from Stasińska et al. are correctly identified. The unconfirmed source of our 13 candidates (id #156, see Table B2) which is not in the Stasińska et al. catalog is consistent with being a PN based on its line ratios which place it well above the extreme starburst lines in the BPT diagram (Fig. 6), and we therefore assign it a PN flag. This leaves 6 of the Stasińska et al. PNe that we do not recover with our selection. Of these, 4 are below the detection threshold and are therefore not recovered in the dendrogram, one falls within an H II region contour, and the last is consistent with being an H II region based on its emission line ratios and its radius of ∼11 pc (which is therefore well resolved in the observations unlike the other PNe).
This shows that our PNe selection method is very robust, considering that with the source detection threshold and other dendrogram pruning parameters we recover all of the known PNe in the field that were picked up in the emission line region identification step. Again, we note that the goal of this study is not to compile a complete census of PNe, but rather to remove the ones that fall within our detection algorithm from the H II region sample. By lowering the detection threshold in the dendrogram analysis and by including the [O III] map to identify regions, the four undetected PNe would likely have been correctly matched. Conversely, this implies that we also do not detect other fainter H II regions, which however does not further impact our analysis given their expected low signal-to-noise ratios and thus high uncertainties.

H II regions
With the SNR and PNe selection criteria described above, the initial emission line region catalog of 188 objects is reduced to 103 spatially resolved sources (i.e. after removing SNRs and PNe we place an additional constraint on the size of the regions by requiring a radius > 7 pc to only include spatially resolved sources, bringing the number down to 103 regions), which are therefore classed as bona fide H II regions and used for the subsequent analyses. Their line ratios are consistent with and place them in the expected BPT diagram space occupied by H II regions. This is illustrated in Fig. 6, where, in addition to the extreme starburst lines from Kewley et al. (2001) and Kauffmann et al. (2003) that are widely used to separate star formation-from AGN-dominated galaxies, we also show the separation line proposed by Stasińska et al. (2006) for local star-forming galaxies.
As mentioned above, cross-matching the resulting H II region catalog with catalogs from the literature is not as straightforward as for PNe and SNRs. While a spatial comparison with the location of the ∼83 H II regions from Deharveng et al. (1988) that fall within the MUSE field-of-view shows a qualitative good agreement, it is clear that the dendrogram+clustering structure identification recovers smaller and more compact regions that do not appear in the Deharveng et al. catalog on the one hand, but on the other hand tends to group together a handful of regions into larger complexes. A more quantitative measure for the robustness of the identification approach in recovering NGC 300 H II regions is to compare population statistics with previous studies, e.g. the H II region H luminosity function, as is shown in Fig. 7. This shows good agreement with Deharveng et al. (1988), in particular in the 10 37 erg s −1 < (H ) < 10 38 erg s −1 regime. At the high luminosity end, the dendrogram approach used here leads to a slight horizontal shift towards higher luminosities due to the grouping of some regions, while at the low luminosity end, the much higher spatial resolution of the MUSE data (compared to Deharveng et al.) leads to a flatter tail. Also shown in Fig. 7  Based on the above, the H II region catalog is likely complete in the regime of small region radii down to the spatial resolution limit of about 7 pc. However, we are likely incomplete in the regime of large region radii, due to the fact that larger and brighter regions are less easily detected than smaller, fainter regions, i.e., for given fluxes (respectively, for given radii), regions with increasing radii (respectively, decreasing fluxes) will eventually fall below the surface brightness cut. To further investigate the completeness, we proceed in a similar fashion to Rosolowsky et al. (2010) and produce a simulated H map, inject Gaussian sources, and run our detection algorithm to assess the source recovery. While Gaussian sources are not ideally matched to the flux distributions observed in the real H map (in which H II regions have a variety of different morphologies, from roughly circular with a central peak to open and filamentary), this remains a worthwhile test to do to assess completeness. The simulated map has the same dimension (i.e., number of pixels, WCS, pixel size, etc.) and noise properties as the observed image used for the structure identification. We then run the region identification algorithm on the synthetic map by injecting Gaussian sources (with radii evenly distributed between the observed values) with increasing surface brightnesses into a non-crowded FOV, to assess the impact of low surface brightness regions. We do not vary pruning parameters (minimum number of pixels for independent structures and the leaf significance) in this test, as these are not driving the (in-)completeness at large region sizes. To estimate the degree of completeness we compare the number of injected to the number of detected sources as a function of input parameters from 10 bootstrap iterations to estimate uncertainties on the detection fraction. The completeness fraction as a function of input Gaussian source flux is shown in Fig. 8, indicating that the catalog is >99% complete for source fluxes >2.5 , with = 1.7 × 10 −18 erg cm −2 s −1 being the mean noise level of the H map. This means that with the pruning parameters set as described in Section 3.1, we recover most sources brighter than >2.5 . In terms of purity, we note that contamination from PNe to the resulting H II region catalog is unlikely, as PN contaminants would have been removed with the size constraint. SNR contamination is more likely, given that very old  Figures 4 and 5) are NGC 300 specific and are likely different in different environments as the strength of emission lines is affected by factors such as metallicity and the resolved stellar population within the regions. To make the presented classification scheme widely applicable, it would need to be augmented with photoionisation models, and tested on spectroscopic data sets of nearby galaxies that differ not only in terms of galaxy properties (e.g. mass, metallicity, type, etc.) but also in terms of distance which affects resolution.

FEEDBACK-DRIVEN GAS IN H II REGIONS
Having compiled a catalog of H II regions, we now exploit the full range of nebular emission lines covered by the MUSE observations  to characterise the feedback-driven gas in the regions. With the aim of exploring the existence of a relation between H II region properties and the impact of massive feedback-driving stars that have formed within them, of particular interest are gas-phase abundances and ionisation properties. For this, in Section 4.1 we first derive key gas properties, and then link these to feedback-related pressure terms within the regions in Section 4.2.

Gas-phase abundances and ionisation properties
To derive gas-phase abundances in H II regions, the direct temperature-based method is generally preferred (e.g. Peimbert et al. 2017). It has been used to derive abundances in NGC 300 by Bresolin et al. (2009, 28 H II regions), Stasińska et al. (2013, 9 compact H II regions), and more recently by Toribio San Cipriano et al. (2016, 7 H II regions), confirming that NGC 300 has a negative metallicity gradient, as discussed in Deharveng et al. (1988). However, this  Bresolin et al. 2002). The MUSE observations used in this work are not sufficiently deep to obtain reliable auroral line detections, and we therefore use the so-called strong line method to determine abundances. This method relies on calibrations obtained either from theoretical calculations or by fitting observed relations between different strong line ratios and abundances derived from the direct method.
There is a rich history of both auroral and theoretical calibrations in the literature for a variety of different metallicity-sensitive strong line ratios, and a review of these is given in Kewley et al. (2019). Here, we use the N2 ratio (≡ [N II] 6584/H ) which, among the possible strong line ratios covered by MUSE, is the least sensitive to reddening and flux calibration issues due to the vicinity of the two lines. More importantly, it makes use of only two emission lines, therefore minimising possible dependencies on other emission lines of correlations discussed below. Another possible strong line ratio available from the MUSE coverage is O3N2 (≡ ([O III] 5007/H )/([N II] 6584/H )), which is valid across a larger metallicity regime than N2, but the line ratio itself suffers from a strong dependence on the ionisation parameter which would propagate to the derived metallicities. The ionisation parameter dependence of the N2 ratio is lower, with the N2-derived metallicity varying by about 1 order of magnitude with the ionisation parameter (Kewley et al. 2019). Composite diagnostic strong line ratios have been proposed to overcome the ionisation parameter dependence (e.g. Dopita et al. 2016, which combines the N2 and [N II]/[S II] ratios), as well as theoretical calibrations which include an ionisation parameter correction (e.g. Kewley et al. 2019). The former introduces a non-negligible amount of scatter when compared to abundances derived from N2 only, and we therefore do not use it. Because most H II region properties computed in this paper are based on strong lines, we refrain from using metallicity calibrations that use many lines in order to avoid introducing systematic dependencies when analysing relations (e.g. the pressure of the ionised gas, H II , is proportional to the [S II] lines, see Section 4.2). Fig. 9 shows the abundance gradient (derived from the N2 ratio) as traced by the H II regions, which, in the central part of NGC 300 covered by the MUSE mosaic, have abundances ranging from ∼2/3 solar to ∼1/2 solar. We perform a linear fit, obtaining 12 + log ( Stasińska et al. (2013), and about 0.15 dex lower than that of Bresolin et al. (2009). We note that differences in slope and intercept with the literature are likely due to two key factors. Firstly, we use the strong line method to derive abundances, while the values we compare our results to have been derived using the direct method. As discussed in Bresolin et al. (2009), different strong line ratios can sometimes lead to drastically different intercepts and slopes, with N2 the one showing the best agreement with the abundance gradient from the direct method. Indeed, the 5007, is temperature-sensitive. Secondly, the MUSE data only cover the inner part of NGC 300 (for reference, the Bresolin et al. study extends to about 1 25 ), and additional MUSE data covering the outer parts of the disc would be needed to better constrain the abundance gradient.
Another key property that can be determined using optical emission line ratios is the degree of ionisation of the gas within the H II regions. With the observed metallicity gradient, this is indeed an interesting quantity to explore, as lower metallicities imply that stars are hotter and have higher photon fluxes, which therefore translate into higher degrees of ionisation. Here, we use the ratio of the two sulphur ionisation states within the MUSE wavelength coverage, S32 ≡ [S III]/[S II], which is considered a good tracer of the degree of ionisation since the actual line ratio itself is significantly less dependent on metallicity than, e.g. We further explore the ionisation properties of H II regions within the covered portion of NGC 300 in terms of the radiation  (2009) but is not covered with sufficient signal-to-noise by our observations, and we therefore distinguish the modified softness parameter as defined by Eq. 2 from the one used by Pérez-Montero & Vílchez by adding the subscript N II. The softness parameter, does not dependent on local effects of the ionisation parameter and density, it is sensitive to the effective temperature of the stars ( ∼ 1/ eff ) within a given region, considering the difference of the ionising potentials involved in the ratio (i.e. ion (S + ) ∼ 10.4 eV, ion (S 2+ ) ∼ 23.3 eV, ion (N + ) ∼ 14.5 eV, ion (O 2+ ) ∼ 35.1 eV), and its value decreases with increasing hardness of the radiation field. Pérez-Montero & Vílchez find that individual nearby galaxies have negative gradients with varying slopes. Bresolin et al. (2009) find a metallicity trend in the hardness of the ionising radiation in NGC 300. While a quantitative comparison with the linear least-square fit found by these authors is not possible due to the modification of the softness parameter used in this work, a qualitative comparison confirms the trend of N II with metallicity ( Fig. 10, lower panel). Together with the trend found for S32, this indicates that stars have higher effective temperatures in the outer H II regions of NGC 300 probed by the MUSE data, producing harder ionising photons and resulting in higher degrees of ionisation. Fig. 10 shows an overall picture of the radial and metallicity trends of the ionisation state and radiation hardness in the H II regions, and we additionally show the [O III]/H ratio as another tracer of the degree of ionisation (middle panel). While linear fits are shown in the various panels, the purpose of this figure is mainly qualitative, and it highlights that a negative trend with galactocentric radius corresponds to a positive trend with metallicity (derived from the N2 ratio) and vice versa. A more quantitative study of the H II regions based on photoionisation modeling will be discussed in a forthcoming publication, but we detail the fit parameters in Table 3. The forthcoming quantitative H II region study will also contain a detailed comparison between abundances derived via the strong line method from the MUSE data and corresponding values obtained via the direct temperature method from Bresolin et al.. These types of comparisons are crucial for refining strong line calibrations (e.g., Curti et al. 2017), in particular given the rise of spatially resolved observations of orders of magnitude more regions in nearby galaxies.

Feedback-related pressure terms
With the picture emerging from the previous section in which lowermetallicity H II regions tend to have harder radiation fields, we now ask the question of whether this has measurable consequences on the impact of stellar feedback in these regions. To assess stellar feedback from the optical IFU data, we proceed as in M20 in computing feedback-related pressure terms (see also Lopez et al. 2014 andOlivier et al. 2021). Specifically, we focus here on the pressure of the ionised gas, H II , and the direct radiation pressure, dir , but do not quantify the effect of stellar winds. Compared to the pressure of the ionised gas, stellar winds have been shown to be increasingly less effective towards evolved H II regions (e.g. Lopez et al. 2014;Barnes et al. 2020;McLeod et al. 2020), and they are expected to be weaker at lower metallicities (Kudritzki 2002), such that the ionised gas pressure will therefore likely dominate the total feedback pressure, in particular at H II region radii 10 pc like the ones observed here . In M20, we demonstrate that stellar wind pressures of NGC 300 H II regions can be derived from the resolved population of O-type stars in the regions, and we are currently working towards identifying individual massive stars from the MUSE data to then be able to quantify the relative importance of stellar winds and ionisation towards total H II region pressures as well. Further, even though dir is only marginally contributing to the current expansion of the H II regions observed here (as is shown later in this section), its impact during the embedded and compact stages of the regions' evolution contributed towards the present day region properties and is therefore certainly worth investigating.
The ionised gas pressure is computed from the electron density, e , as where we assume an H II region temperature of e = 10 4 K as well as singly-ionised helium. With H II region temperatures in NGC 300 ranging between ∼8000 K and ∼12000 K (Bresolin et al. 2009), this is a reasonable assumption. The direct radiation pressure is evaluated from the combined (bolometric) luminosity, bol , of the massive stellar population in each region, where is the H II region radius (which is derived from as the geometric mean of the major and minor axes of the projection onto the position-position plane, i.e. of an ellipse that describes the identified region, and for which we assume a measurement uncertainty of 20 per cent). To compute the bolometric luminosity, we assume that the IMF and stellar age distribution are fully sampled in the H II regions and convert the extinctioncorrected H luminosity to the bolometric luminosity according to Kennicutt & Evans (2012), bol = 138 (H ). As discussed in Lopez et al. (2014), there are caveats to the conversion from H to bolometric luminosity for young star-forming regions and stellar populations. First, because we are likely tracing stellar populations that are younger than ∼5 Myr, for which the ratio of H to bolometric luminosity is higher than the Kennicutt & Evans value, we are likely overestimating the bolometric luminosity, and this overestimation will be different for regions of different ages. Second, the assumption of fully sampling the IMF is no longer valid for star clusters with masses below ∼10 4 M , introducing stochastic variations of the H to bolometric luminosity ratio, which for a randomly selected cluster would imply an underestimation of bol (e.g. Krumholz et al. 2015;Haydon et al. 2020). Similar to Lopez et al. (2014), the selection of the H II region sample investigated here is not random but based on H emission, mitigating the stochastic effect to a factor ∼2 in the level of uncertainty. Fig. 11 explores the impact of environment on the feedbackrelated pressure terms computed above. As already found by Lopez et al. (2011), McLeod et al. (2019a, and in M20, the predominantly more evolved H II regions studied here are dominated by the pressure of the ionised gas, and we find H II / dir > 1 (the median across all H II regions is ∼56). While a clear correlation cannot be drawn and correlation coefficients reveal weak relations at best given the large scatter (with correlation coefficients of ∼0.20 and 0.18 for H II and dir , respectively), both H II and dir tend to increase with increasing galactocentric radius, along the negative metallicity gradient. To further illustrate the trend we divide the regions into five radial bins, the mean value of which are plotted against galactocentric radius in Fig. 11, and these having correlation coefficients of ∼0.72 and ∼0.64 for H II and dir , respectively. Here, we use the distance from the centre of NGC 300 as an indicator of the varying ISM conditions, given that we have shown in Section 4.1 that along with the gas-phase abundance, the degree of ionisation and the hardness of the radiation field also show radial trends. We do not, however, explicitly plot the pressure terms against these, because similar emission lines and line ratios go into computing the different quantities, thus contaminating any resulting relations. For example, dir and metallicity are both dependent on the H flux, while H II depends on the [S II] lines, as do N II and S32. Independent measurements are required to attempt to quantify the metallicity dependence of dir and H II , for example by obtaining deeper H II region spectra and derive abundances from the direct method, by using different line ratios, or by combining emission line measurements with photoionisation modeling.
An increase in H II with decreasing metallicity can be explained with the harder ionising radiation and increased photon fluxes. Further, while we have assumed a fixed temperature for all H II regions, the negative metallicity gradient in NGC 300 leads to a positive temperature gradient due to line cooling being less efficient at lower metallicities (Bresolin et al. 2009). If we were to take this into account, the increase of H II towards larger galactocentric distances would be further enhanced. The increase of dir is not as easily understood. The well established relation between the gas-to-dust ratio (G/D) and metallicity shows that the dust content tends to decrease with decreasing metallicity (e.g. Rémy-Ruyer et al. 2014, and references therein), and we would therefore expect a decrease in radiation pressure at lower metallicities due to less dust to impart momentum to. This is indeed the case when comparing dir values of H II regions in the LMC and SMC studied by Lopez et al. (2014): the SMC H II regions (where the metallicity and dust content are lower than in the LMC, e.g. Roman-Duval et al. 2014) have systematically lower radiation pressure values than the LMC regions. We assess the amount of dust towards the H II regions in NGC 300 from the reddening correction, as returns both the colour excess, ( − ), and the extinction, , based on the measured H /H ratios. Fig. 12 reveals a positive H II region extinction gradient in the part of NGC 300 covered by the MUSE mosaic, which is in agreement with Casasola et al. (2017) who found that the dust mass surface density in NGC 300 is slightly increasing up to about ∼0.5 25 (their fig. A.2). At the larger galactocentric radii not covered by MUSE the dust mass surface density decreases (Casasola et al. 2017), and additional MUSE data probing regions beyond 0.45 25 are required to provide additional insight on the dependence of dir on the dust content. We note that this does not go against the G/D-metallicity relation, given that the metallicity range probed by the H II regions in the inner ∼2.5 kpc of NGC 300 is rather small, i.e. ∼0.2 dex, range in which the scatter around mean G/D values can be large (Galametz et al. 2011;Rémy-Ruyer et al. 2014) such that a correlation within the small metallicity range cannot be established. We suggest that the dir − relation shown in Fig. 11 is driven by the dust content, together with the fact that at lower metallicity result in higher photon fluxes.
We further explore the role of metallicity and dust content in regulating pressure terms with the radiation-hydrodynamical models of Ali (2021). These models used the Monte Carlo radiative transfer code (Harries et al. 2019) to explore stellar feedback in clusters at four different metallicities (2, 1, 0.5, 0.1 Z ). This method calculated photoionisation equilibrium and radiation pressure (direct plus dust-processed pressure) step-by-step with hydrodynamics in 3D. Fig. 13 shows the mean thermal pressure and radiation pressure in H II regions extracted at different evolutionary times. The left panel shows results from the original models of Ali (2021) in which the gas-to-dust ratio is proportional to metallicity (G/D ∝ 100/( /Z )), therefore serving as the control run where the dust content and the metallicity both decrease. The right panel shows post-processed snapshots where the radiative transfer has been recalculated with a fixed G/D = 100 for all metallicities (i.e. boosting the dust content at sub-solar metallicity and lowering it at super-solar metallicity), therefore emulating environments similar to those observed in the covered portion of NGC 300 in terms of dust mass surface density. Together with the observations, these preliminary results indicate that the radial variation of the dust content is indeed likely playing the dominant role in regulating the dir − relation: an increased amount of dust results in greater extinction of UV photons, making the H II regions smaller and the radiation pressure larger. With simulations in general adopting the canonical G/D-metallicity relation, we therefore note that across spatially resolved scales within observed systems the dust content can in fact remain constant or even increase despite the presence of a negative metallicity gradient, which therefore can impact the relative importance of radiation pressure. This is likely even more important in young, compact H II regions, and should be taken into account when simulations are tailored to match specific observed systems.
We note that the radial pressure trends observed here are however a combination of several intrinsic and environmental factors, e.g. gas/dust content, metallicity, star formation rate, evolutionary stage of the H II regions, age of the stellar population, the relative importance of which needs to be assessed. For example, the star formation rate (a relative increase of which over a given period of time would result in overall increased stellar feedback) is observed to be approximately constant in the portion of the galaxy covered by our observations (out to about 0.45 25 ; Gogarten et al. 2010;Williams et al. 2013;Casasola et al. 2017), and we therefore do not consider it to be a dominant source driving the increase of stellar feedback here. In terms of evolutionary stage, one would expect the radiation pressure to be enhanced in younger, more compact regions when compared to more evolved regions (Olivier et al. 2021). With a radius cutoff of about 7 pc as described we are not able to probe compact and ultra-compact H II regions, and an evolutionary dependence as traced by the ages of the stellar populations within the regions will be explored in a forthcoming paper. Further, the wealth of existing and upcoming IFU data from nearby galaxy surveys like SIGNALS or PHANGS which contain targets with substantial ancillary, co-spatial multi-wavelength coverage is ideal to study environmental dependencies of stellar feedback. Current substantial effort towards this will certainly lead to further insight in this area in the near future.

SUPERNOVA REMNANTS IN THE CONTEXT OF EARLY STELLAR FEEDBACK
Supernova feedback has long been considered to be the main mechanism responsible for driving turbulence and regulating star formation rates and efficiencies on galaxy-wide scales. However, in recent years simulations have started to show that early, pre-SN feedback (radiation pressure, ionisation, stellar winds) need to be accounted for because the energy deposited by SNe is not sufficient to disrupt and disperse dense molecular clouds (see Section 1). For NGC 300, we applied the statistical methodology of Kruĳssen & Longmore (2014) and Kruĳssen et al. (2018) to characterise the GMC lifecycle and found that the GMCs in this galaxy are dispersed on short time-scales of 1.5 ± 0.2 Myr, requiring pre-SN feedback . This result has been generalised in recent observational studies by Chevance et al. (2020a,c), who have quantified feedback time-scales across nine nearby disc galaxies, finding that GMCs are dispersed within 1−5 Myr after the emergence of massive stars from their dust-enshrouded birth places. This analysis has been extended further by Kim et al. (2021) to also encompass the dust embedded stages of star formation than those probed by Chevance et al., and both studies conclude that early (pre-SN) stellar feedback (in the form of stellar winds and photoionisation in particular) is a major component driving the GMC disruption. These results are consistent with optical/near-ultraviolet studies of the young cluster population, which find that clusters become unassociated with their natal clouds after just a few Myr (e.g. Hollyhead et al. 2015;Grasha et al. 2018;Hannon et al. 2019). Both Chevance et al. and Kim et al. do not find significant correlations between environmental properties and the time-scales over which feedback acts.
As described in Section 4, the H II regions in the inner portion of NCG 300 reveal a trend of increasing pre-SN feedback (traced by the direct radiation pressure and the pressure of the ionised gas) with increasing galactocentric radius, i.e. along NGC 300's negative metallicity and positive extinction gradient. The massive stellar populations residing in the H II regions are not only the drivers of early stellar feedback, but are also the progenitors of SN events. In other words, (core collapse) SNe typically occur in H II regions within a few Myr of the onset of pre-SN feedback, which has already started to affect the surrounding ISM and alter the environment into which SNe expand into (e.g. Haid et al. 2018;Lucas et al. 2020;Keller et al. 2021). Hence, we now ask the question of whether there are any systematic differences in terms of the environment of the seven detected SNRs, and if yes, what the likely driver of the environmental differences is. To probe the pre-SN environment we compute the pre-shock ISM density, ISM (the density of the ISM in which SNe went off into, which differs from the density as one would measure from, e.g. the ratio of the [S II] lines, as this would deliver the electron density of the shocked matter), proceeding as in McLeod et al. (2019a): We use the relation between the flux of the H line from a surface element of a radiative shock with velocity, s , and ISM density, ISM , as given in Dopita & Sutherland (1996) With the shock velocities measured from the [S II] FWHM of the SNRs ( s ≈ FWHM; Heng 2010), we derive pre-shock ISM densities on the order of a few particles per cm −3 (see Table 4), values that are consistent with shocks in a pre-ionised medium (Dopita & Sutherland 2017).
With an estimate of the pre-shock ISM density for each SNR, we now explore whether there is any trend with galactocentric radius and/or metallicity, thus testing the scenario in which the enhanced early stellar feedback at these galactocentric radii (i.e. metallicities) could have created less dense environments in which subsequently the supernova events occurred. Fig. 14(a) shows the inferred preshock ISM density, ISM , as a function of galactocentric radius, together with two linear fits (with and without taking into account the uncertainties on ISM in magenta and orange, respectively). While there appears to be a relation in the expected direction, the limited number of SNRs in the MUSE mosaic is likely hindering a more robust interpretation. To further explore the significance (or lack thereof) of the relation between ISM and , we compute the Pearson correlation coefficient which, for the ISM − values shown in Fig. 14(a) is about −0.83, indicating a relatively strong negative relation. We assess the uncertainty of by bootstrapping while: (i) randomising and (ii) varying the ISM values within their respective uncertainties.
With (i) we test whether reshuffling the -axis yields the original relation between ISM and , hence testing the null-hypothesis of there being no relation between ISM and . By bootstrapping this with 10 4 iterations, we thus obtain a statistical significance of nearly 3 for the relation. Or, in other words, by randomising we do not recover the relation 99.64 per cent of the times. With (ii) we quantify the relation given the uncertainties on ISM , thus computing the uncertainty of . The correlation coefficient analysis for the ISM − relation is shown in Fig. 14(b), where the black histogram corresponds to randomising , while the teal histogram corresponds to varying ISM : within 1 of , the correlation is 2 away from the null-hypothesis representing no relationship. We further analyse the environmental dependency of the pre-shock ISM density in Fig. 14(c) and Fig. 14(d), which show that a corresponding positive correlation is found with metallicity, as would be expected if increased pre-SN feedback at lower metallicities created lower density environments. While the correlation between the pre-shock ISM density and the gas-phase metallicity, 0.71, is not as strong as the one with galactocentric radius above ( −0.83), we still recover a 2 significance within 1 of .
These findings are in excellent agreement with the results of Lucas et al. (2020), who simulate supernovae in star-forming molecular clouds with and without early stellar feedback (in the form of ionisation and stellar winds). In their study, Lucas et al. find that early stellar feedback creates pre-SN cavities with systematically lower densities. Further, compared to the control run, they find that including pre-SN feedback leads to simulated clouds with enhanced low column density channels via which the SN ejecta (and shockheated gas) can escape more easily, thus affecting the lower density ISM at greater distances from the natal cloud. This means that because of early stellar feedback, SNe can deposit more energy on galactic scales. Conversely, Smith et al. (2021) find that including early stellar feedback in their simulations disrupts molecular clouds and leads to less SN clustering, reducing outflow rates in a substantial manner. Keller et al. (2021) apply an empirically-motivated, early feedback model based on the observations of Kruĳssen et al. (2019) and Chevance et al. (2020a,c) in isolated disc galaxy simulations, and also find that the inclusion of early feedback reduces SN clustering, but next to a moderate decrease of the time-averaged outflow, its burstiness drops precipitously, leading to a considerably 'smoother' galaxy-scale baryon cycle. While our findings do not allow us to quantify the impact of the SNe themselves, these numerical studies clearly demonstrate the importance of the measurement in Fig. 14.
The environmental dependencies found here are consistent with the results of Chevance et al. (2020c) and Kim et al. (2021), who observed no statistically significant radial trends for the derived feedback time-scales within the uncertainties. Indeed, the increase of feedback pressure reported here does not directly imply accelerated feedback time-scales, as the time-scales depend on other environmental properties that set how rapidly the imparted feedback (which is what we measure) disperses the natal GMCs. In fact, our findings support the weak relation between metallicity and feedback time-scales found by Chevance et al. where lower metallicity galaxies have shorter feedback time-scales. Furthermore, this underpins the results of Kruĳssen et al. (2019), who showed that the GMC lifetimes and feedback time-scales both decrease with increasing galactocentric radius in NGC 300 and that these are regulated by early stellar feedback. An accelerated feedback timescale towards larger galactocentric radii would imply that one would observe less evolved regions at larger radii, which could in turn contribute to the radial increase of dir as it is enhanced in younger regions.
We caution that while we show direct evidence for an environmental dependency of the pre-shock ISM density in the form of galactocentric radius and metallicity, a causal connection between enhanced dir and H II and lower pre-shock ISM densities, while tempting, cannot be directly made from the observations, so that this connection is purely inferred. We also stress that while the correlation coefficient analysis indicates that the ISM -environment relation is relatively robust, additional MUSE coverage of more SNRs in NGC 300 at larger radii and lower metallicities would be needed to further strengthen the result. Two additional points need to be made. First, while the SNR identification is robust, the sample is potentially missing old, evolved SNRs with low surface brightness and unresolved shock velocities. Second, with the focus here being on SNRs, this analysis does not allow the investigation of SNe that potentially already occurred within the H II regions but are not identified here due to the dominating H II region spectra. More clarity towards this will be obtained by analysing the age of the stellar populations within the regions.

SUMMARY AND CONCLUSIONS
In this paper, we have used optical IFU data from the VLT/MUSE instrument to identify and classify emission line regions, and study the environmental properties of H II regions and SNRs. The emission line region identification procedure exploits dendrograms, which are complemented with a spectral clustering algorithm which ensures that the spatially coherent structures are not over-fragmented. The main drawback is, like it is the case for most region identification methods, that the algorithm can sometimes group together into a single complex what are likely separate regions, particularly where several emission line regions are crowded together and spatially overlap. Given that we are mostly analysing radial dependencies, grouping together neighbouring regions only has a marginal impact on the analysis performed here.
We separate regions into SNRs, PNe, and H II regions based on emission line ratios and emission line widths, and show that this method is robust and can be readily used for nearby galaxy data sets of similar spectral and spatial coverage. We do however recommend that if the aim is to specifically identify SNRs and PNe, additional emission line maps to the H one should be used (e.g. [S II] for SNRs and [O III] for PNe) to ensure that objects with low H surface brightness are also recovered. We stress that this method is purely empirical, and a comparison with other methods is planned for the future.
For the identified H II regions we then derive several quantities associated with the ionised gas within them. These are the (oxygen) abundance, the degree of ionisation, and the hardness of the radiation field, and we recover the known negative metallicity gradient of the galaxy. We show that, within the portion of the galaxy covered by the MUSE mosaic (i.e. out to ∼0.45 25 ), the degree of ionisation increases and the radiation fields become harder with increasing galactocentric radius and decreasing metallicity. This is consistent with the stellar populations residing in the H II regions having higher photon fluxes and higher effective temperatures at lower metallicities, as inferred from higher degrees of ionisation and harder radiation fields. We then compute feedback-related pressure terms for ionisation and radiative feedback and discuss their dependence but also impact on their environment. With the strength of early stellar feedback seemingly increasing as a function of distance from the galactic centre and decreasing metallicity, we then analyse the environment into which SNe in NGC 300 occurred.
The main results from this paper are summarised as follows: • and the hardness of the H II region ionisation fields, such that at lower metallicities (i.e. larger galactocentric radii) we find harder radiation fields and increased degrees of ionisation indicating stars with higher effective temperatures and higher ionising photon fluxes. Both photoionisation modeling and auroral line measurements are needed to further disentangle intrinsic dependencies of the used strong line ratios to characterise the ionised gas.
• Linked to the previous point, we observe weakly increased pre-SN feedback (traced by the direct radiation pressure and the pressure of ionised gas) in H II regions at larger galactocentric radii and lower metallicities.
• We suggest that the increased radiation pressure in H II regions at larger galactocentric radii and lower metallicities is likely due to an increase in the dust content towards the outer regions probed by the data. Additional MUSE observations of the outer parts of NGC 300 (beyond the presently available coverage which extends to ∼0.45 25 ) are needed to further quantify the dependence of dir on the dust content.
• Preliminary results from dedicated simulations of star-forming molecular clouds support our conclusion that dust regulates the relative importance of radiation pressure. We suggest that, where available, the known dust mass surface density should be taken into account if simulations are tailored to reproduce specific galactic systems.
• Given the radial and metallicity trends of pre-SN feedback we further study the density of the environment the detected SNR expanded into, and find that SNe at larger galactocentric radii (and at lower metallicities) expanded into lower-density environments. We tentatively suggest that this could be a consequence of the increased pre-SN feedback at low metallicities, which has contributed to creating lower-density environments, in excellent agreement with simulations.
In conclusion, we reiterate that ongoing and upcoming IFU surveys of nearby galaxies enable H II region and SNR studies analogous to what is presented in this paper, while providing a much wider range of different environmental properties to explore. These will then provide an observational quantification of the environmental impact and dependence of stellar feedback, which can be used to improve feedback prescriptions in simulations. Table 4. Properties of SNRs (see text Section 5). Columns correspond to the dendrogram id (1), the galactocentric radius (2), the H intensity (3), the inferred shock velocity (4), the pre-shock ISM density (5) and the SNR gas-phase abundance (6).

DATA AVAILABILITY
The data underlying this article will be shared on reasonable request to the corresponding author.  (8.43, 8.46, 8.37, and 8.34, respectively) are indicated with vertical lines.

APPENDIX A: ABUNDANCES FROM STRONG LINE RATIOS
In Section 4.1, we derive strong line oxygen abundances for the detected H II regions. For this we use the N2 ratio and the Marino et al. (2013) calibration. However, the N2 ratio itself has a slight dependence on the ionisation parameter, and here we explore the difference between using the empirical Marino et al. calibration and using the theoretical calibration by Kewley et al. (2019), which corrects for the ionisation parameter dependence. In addition, we also compare these two to the Marino et al. (2013) O3N2 and the Pilyugin & Grebel (2016) calibrations. Fig. A1 shows the cumulative distributions of the oxygen abundance derived from the two N2 calibration, for which a KS test yields a -value of 1.4 × 10 13 . Compared to Marino et al.,

APPENDIX B: EMISSION LINE REGION PROPERTIES
In the following tables, we report the properties of the identified SNRs (Table B1), PNe (Table B2), and H II regions (Table B3).  This paper has been typeset from a T E X/L A T E X file prepared by the author.