MUSE crowded field 3D spectroscopy in NGC300 I. First results from central fields

Aims. As a new approach to the study of resolved stellar populations in nearby galaxies, our goal is to demonstrate in NGC300 that integral field spectroscopy with high spatial resolution and excellent seeing conditions reaches an unprecedented depth in severely crowded fields. Methods. MUSE observations with seven pointings in NGC300 have resulted in datacubes that are analyzed in four ways: (1) PSF-fitting 3D spectroscopy with PampelMUSE yields deblended spectra of individually distinguishable stars. The technique also provides samples of planetary nebulae that are complete down to m5007=28. (2) pseudo-monochromatic images, created at the wavelengths of the most important emission lines and corrected for continuum light by using the P3D visualization tool, provide maps of HII regions, SNR, and the diffuse ISM at a high level of sensitivity, allowing for the discovery of planetary nebulae, WR stars etc. (3) The use of the P3D line-fitting tool yields emission line fluxes, surface brightness, and kinematic information for gaseous objects, corrected for absorption line profiles of the underlying stellar population. (4) Visual inspection of the datacubes is demonstrated to be effcient for data mining and the discovery of background galaxies and unusual objects. Results. We present a catalogue of luminous stars, rare stars such as WR and other emission line stars, carbon stars, symbiotic star candidates, planetary nebulae, HII regions, supernova remnants, giant shells, peculiar diffuse and filamentary emission line objects, and background galaxies, along with their spectra. Conclusions. The technique of crowded-field 3D spectroscopy is capable of deblending individual bright stars, the unresolved background of faint stars, gaseous nebulae, and the diffuse component of the interstellar medium, resulting in unprecedented legacy value for observations of nearby galaxies with MUSE.


Introduction
The quest for understanding the formation and evolution of galaxies has provided us with a wealth of data from imaging and spectroscopic surveys, most notably the Sloan Digital Sky Survey (SDSS), to quote just one prominent example (York et al. 2000;Alam et al. 2015). As a more recent development, integral field spectroscopy (IFS) enables us to obtain spatially resolved information on stellar populations, gas, dust, and the star formation history, including spatially resolved kinematics across the entire face of a galaxy, rather than the limited apertures of a single fibre, or a slit (see Sanchez et al. 2012). The CAL-IFA survey is the first effort to fully exploiting this capability with a reasonably large sample of galaxies of all Hubble types (Walcher et al. 2014;Husemann et al. 2013;Garcia-Benito et al. 2015). Surveys with deployable integral field units (IFU) like SAMI (Cortese et al. 2014) and MaNGA (Bundy et al. 2015) are carrying this approach further to much larger sample sizes. However, given the integration of light from typically hundreds of thousands of individual stars over a spatial resolution element (spaxel) of the respective integral field units, it is not entirely clear whether the technique of stellar population synthesis for the resulting spectra is capable to uniquely determine stellar populations and their star formation history. The study of resolved stellar populations in nearby galaxies, where stars can be measured Article number, page 1 of 37 arXiv:1806.04280v1 [astro-ph.GA] 12 Jun 2018 A&A proofs: manuscript no. NGC300_MUSE_2column_arXiv accurately on an individual basis, would be able to lift any degeneracies and allow to calibrate the integral light measurements of the more distant galaxies. Such work has been attempted with Hubble Space Telescope (HST) photometry, e.g. in the ANGST survey (Dalcanton et al. 2009), for the galaxy NGC 300 studied in this paper see specifically Gogarten et al. (2010). However limitations of photometry, e.g. the age-metallicity degeneracy, sensitivity to extinction, or the lacking ability to distinguish objects such as Wolf-Rayet (WR) stars from normal O stars, carbon stars from M stars, etc. would rather call for spatially resolved spectroscopy. Except for conventional slit spectroscopy of individual bright giants or supergiants in local volume galaxies, e.g. Bresolin et al. (2001), Smartt et al. (2001) in the optical, or Urbaneja et al. (2005), Gazak et al. (2015) in the NIR, there exists to date no comprehensive spectroscopic dataset for stars in heavily crowded fields of nearby galaxies beyond the Magellanic Clouds. Pioneering attempts to utilize the technique of IFS for crowded field spectroscopy, i.e. fitting the point spread function (PSF) to point sources through the layers of a datacube, analogous to CCD photometry in a single image, e.g. DAOPHOT (Stetson 1987), have demonstrated that it is indeed possible to obtain deblended spectra for individual point sources in heavily crowded fields in local group galaxies , e.g. Becker et al. (2004), Roth et al. (2004), Fabrika et al. (2005). These first attempts were limited by the small field-of-view (FoV) and the coarse spatial sampling provided by first generation IFUs, e.g. 8 × 8 arcsec 2 at 0.5 arcsec sampling in the case of PMAS (Roth et al. 2005), the light collecting power of 4m class telescopes, and less than optimal image quality in terms of seeing.
The advent of MUSE as a second generation instrument for the ESO Very Large Telescope (VLT) has completely changed the situation: the IFU features a 1 × 1 arcmin 2 FoV with spatial sampling of 0.2 arcsec, high throughput, excellent image quality and a one octave free spectral range with a spectral resolution of R=1800 . . . 3600 (Bacon et al. 2014). In preparation for new opportunities with MUSE, Kamann (2013) developed the PampelMUSE software and validated the PSF fitting technique on the basis of PMAS data obtained at the Calar Alto 3.5m telescope (Kamann et al. 2013), measuring for the first time accurate radial velocities of individual stars in the very central crowded regions of the globular clusters M3, M13, and M92, thus putting stringent limits on the masses of any putative intermediate-mass black holes in the cores of those clusters (Kamann et al. 2014). These first experiments were taken to another level with commissioning data from MUSE, where the globular cluster NGC6397 was observed as a 5 × 5 mosaic of ≈ 60s snapshot exposures, resulting in a spectacular dataset of spectra for more than 12000 stars. By fitting to a library of PHOENIX spectra (Husser et al. 2013) and with photometry, good results for the atmospheric parameters T eff and log g, as well as metallicities were obtained, such that a spectroscopic Hertzsprung-Russell diagram (HRD) could be presented for a globular cluster for the first time (Husser et al. 2016). In a kinematic analysis of the same dataset, Kamann et al. (2016) study again the radial velocities and velocity dispersion of individual stars, however with a sample size of more than 12000 objects. Despite the low spectral resolution of MUSE, it is shown that for spectra with sufficient signal-to-noise ratio (S/N), the radial velocities are measured to an accuracy of 1 km s −1 .
Given these encouraging first results, we are now attempting to extend the technique of crowded field 3D spectroscopy to nearby galaxies with the goal of providing a new element to the current state of the art of stellar population synthesis, namely individual spectroscopy of the brightest stars: AGB and RGB stars, the most massive O stars, blue and red supergiants, lu- minous blue variable (LBV) and WR stars. As opposed to the Milky Way, where severe selection effects are at work, e.g. dust extinction, small numbers of known objects, lack of accurate parallaxes (see Crowther 2007), one can hope to obtain complete samples down to a given limiting magnitude (although some of these stars are rare), and to then study their distribution across the face of the galaxy. Unlike the classical technique of broad-/narrowband filter imaging, and then follow-up spectroscopy (see Massey et al. 2015), IFS allows one to obtain images and spectra homogeneously in a single exposure. Filters with arbitrary transmission curves can be devised in the process of data analysis after the observation to reconstruct images from the data cube, which not only enables one to define extremely narrow filter bandwidths and then obtain a very high sensitivity for continuum background limited emission line objects, but also to create efficient notch filters to mask out disturbing emission lines from reconstructed broadband images. Next to the stars, we are then also targeting planetary nebulae (PN) as potentially useful indicators of the underlying stellar population (Ciardullo 2010; (Fesen et al. 1985), the distribution and physical properties of H ii regions as well as the discovery of compact H ii regions, in particular, and the diffuse and filamentary interstellar medium (ISM) are further pieces of information that can be retrieved from the datacube.
As part of the coordinated guaranteed observing time (GTO) program of the MUSE consortium, we have launched a pilot study on the galaxy NGC 300, that is located in the foreground of the Sculptor group, in order to demonstrate the feasibility of the proposed approach. Basic parameters of this galaxy are listed in Table 1. Here we report the first results of our GTO observations in ESO periods P93, P94, P95, also demonstrating new legacy value opportunities provided by MUSE. The reader shall be convinced that any datacube hosts copious amounts of information in addition to the objects of the science case proper that has lead to the observation in the first place. For the purpose of a proof-of-principle, we present first selected results in order to underpin this claim. In future follow-up papers we shall explore the completed data set with a focus on different science cases (see section 6).
The paper is organized as follows: section 2 presents details of the observations and data reduction, section 3 describes the data analysis, section 4 summarizes the results with catalogues of stars, planetary nebulae, H ii regions, SNRs, and serendipitous discoveries, along with a discussion within each class of objects. The paper closes with conclusions and an outlook to the next steps of this project.

Observations and data reduction
Observations were made with the multi unit spectroscopic explorer instrument (MUSE; Bacon et al. 2014), which is placed at the Nasmyth focus of the UT4 8.2m telescope at the Very Large Telescope observatory (VLT) of the European Southern Observatory (ESO) in Chile. NGC 300 was observed as part of guaranteed time observations of the MUSE instrument-building consortium during the three periods P93, P94, and P95. All observations were performed in the extended mode, which includes the wavelength range 4650-9300 Å at a resolution of about 1800-3600 and a dispersion of 1.25 Å pixel −1 . Effects of second-order stray light are clearly seen with this mode in continua with blue emission for wavelengths λ > ∼ 8000 Å.
The fields were chosen to cover as much as possible of the center regions of NGC 300, while avoiding the brightest H ii regions seen in the Hα image of Bresolin et al. (2009). Our aim was to observe each field for a total of 5400 s, which was achieved with three observation blocks (OBs) of 1800 s each, and each OB was in turn split in two exposures of 900 s. To allow for an accurate normalization of the data, the IFU was shifted off the center and rotated between the exposures. Only the first four (two) values were used when only two (one) OB was observed. The sky was observed for two minutes in each OB using the same field on the sky, centered 7 West and 8 41 South of field a. Observations of the standard stars EG 21, Feige 110, EG 274, and GJ 754.1A allowed a spectrophotometric calibration of the data. Daytime calibrations of the morning after the observing night were used. The footprint of the selected fields is outlined in Fig. 1, and the observations journal is shown in Table 2.
Data were reduced with the MUSE pipeline version 1.0 (Weilbacher et al. 2012) through the MUSE-WISE framework (Vriend 2015). The standard procedure was followed for the ba-sic calibration. The first three steps are: ten bias images are combined to make a master bias image, five continuum exposures are combined to make a master flat-field image, and one exposure of each of the three arc lamps (HgCd, Xe, and Ne) are used to derive a wavelength solution. Eleven flat-field exposures of the sky, taken during the evening twilight preceding the science exposures, were combined and used to create a three-dimensional correction of the illumination where 4650 < ∼ λ < ∼ 8000 Å. Detector defects were removed by means of a bad-pixel table.
We applied the calibration products to the science exposures as well as the standard star exposures. The average atmospheric extinction curve of Patat et al. (2011; included in the pipeline) was used when creating the sensitivity functions. The sensitivity functions and the astrometric calibration were applied to each exposure individually. The pipeline automatically corrects the data for atmospheric refraction (using the equation of Filippenko 1982) and the barycentric velocity offset. The sky exposures in each OB were used to create a sky spectrum as an average of the total of 90000 spectra of the sky cube, that was subtracted from the extracted data set of the respective OB. The two to six extracted data sets of each field were combined into a single cube. Owing to the IFU rotation of the observing strategy, the data were affected by the derotator wobble, which is why each exposure needed to be repositioned. Each field contains plenty of stars that could be used to this purpose. The resulting cubes were created using the standard sampling of 0 . 2 × 0 . 2 × 1.25Å.

Analysis
The data analysis was done in two complementary approaches: automatic source detection of point sources (stars, PNe) on the one hand, and visual inspection for the discovery and characterization of spatially extended objects and emission line sources on the other hand. The former was done with the PampelMUSE PSF-fitting code for IFS (Kamann 2013), whereas the latter was based on inspection and processing using the P3D tool, and maps individually extracted from datacubes using DS9, respectively.
P3D is an open source IFS software package developed at the AIP and publicly available under GPLv3 from SourceForge 1 (Sandin et al. 2010). Previous applications of P3D for the VI-MOS and FLAMES IFUs at the VLT were presented in Sandin et al. (2011). In the following paragraphs we describe in more detail how these approaches were accomplished.

Mapping stars and gaseous nebulae
The inventory of stars, gas, and dust is most intuitively mapped by coadding a number of suitable wavelength bins of the data cube to then form broadband or narrowband images. As the MUSE free spectral range spans (in the extended mode) one octave from 465 to 930 nm, there is an almost unlimited choice of synthetic filter curves that can be realized. We have chosen to define three broadband transmission curves V, R, I, that are similar to the Bessell filter curves (Bessell 1990) in that they cover approximately the same wavelength intervals. We have chosen top-hat functions for simplicity and did not attempt to mimic the exact slopes of the transmission curves. Fig. 2   wide enough to cover the instrumental profile of MUSE, respectively, in order to not lose any flux. The filter width, however, was chosen small enough to suppress continuum background light as much as possible. With a typical narrowband filter width of 5 spectral bins (6.25 Å) MUSE is therefore superior in sensitivity to any conventional narrowband imaging camera based on interference filters with a full width at half-maximum (FWHM) of typically ≈ 30 Å. In order to account for the residual continuum flux collected by these filters, a continuum correction was applied to the emission line fluxes on the basis of continuum estimates shortward and longward of the line center. The synthetic filter parameters for the emission lines and corresponding continuum estimates used for our analysis, as well as for the broad band filters, are listed in Table 3.

Extracting point source spectra
To extract spectra of bright individual stars and emission-line objects in our MUSE data, we used the tool PampelMUSE (Kamann et al. 2013). The tool relies on a reference catalog with relative positions and (initial-value) brightnesses of the objects. In a first approach, we used the available stellar-source catalog of the ACS nearby galaxy survey (ANGST 2 ). While some of our fields are completely (fields b and c) -or nearly completely (fields a and j) -covered by the ANGST catalog, three fields are only partially covered (fields d, e, and i), which is why an alternative approach was needed to extract spectra in all parts of all fields. We used the FIND tool of DAOPHOT (Stetson 1987) to determine the centroids of stellar images and create a source catalog. We again used FIND in both a blue and a red part of the continuum to locate as many point sources as possible in each data cube. We set the intensity of each object to a faint value and man- Fig. 2. VRI images reconstructed from MUSE datacubes for pointings (a), (b), . . ., (j). The seeing FWHM for each pointing is indicated, the best image quality with 0.6" (green frame) in (i). Fields (b) and (c) are affected by poor seeing (1.2" and 1.0", respectively). Owing to a reddish colour, an old stellar population is apparent in the nuclear region sampled in (a) and (b), whereas the spiral arm fields (c), (d), (e), and (i) are dominated by blue stars. The interarm region in (j) features extinction by dust lanes. The background image is from the ESO/MPI 2.2m WFI (credit: ESO).
ually added a set of additional objects with a brighter value that were used as PSF stars. All object coordinates were tied to those of the ANGST catalog, which covers at least a few sources in each field. Notably, the alternative approach only allows a spectrum extraction of the brightest blue and red giant stars, while all fainter appearing stars are part of the background. Concerning emission line point sources, we again used find to specify the positions of potential planetary nebulae or compact H ii regions after summing up five pixels on the dispersion axis about [O iii] λ5007 in the respective data cube; additional faint blobs that find missed -which could also be sources -were added to the resulting catalog manually after a visual inspection of the respective data cube.
The extraction of the stellar spectra is performed in a procedure that includes several steps. The analysis begins with an initial guess of the MUSE data PSF, which is modeled as an analytic Moffat profile with up to four free parameters: the FWHM, the kurtosis β, the ellipticity e, and the position angle θ, all of which may depend on the wavelength. Initially, the FWHM is set to the seeing of the observations, β = 2.5, and e = 0. By combining the PSF model with the stellar positions and bright-nesses in the reference catalog, a mock MUSE image is created for a pre-defined filter. Another image is created by integrating the MUSE cube in wavelength direction over the same filter curve. By cross-correlating the two images, initial guesses for all catalogued sources are obtained. The subset of those sources for which meaningful spectra can be extracted is then identified. This is done by estimating the S/N of each source based on its magnitude in the reference catalog, the PSF, and the variances of the MUSE data. In addition, the density of brighter sources around the source in consideration is determined. Only those sources are used in the further analysis where S/N> 5, and where the density determination yields less than 0.4 sources of similar or greater brightness per resolution element. The brightest and most isolated of the selected sources are flagged as PSF sources. They are used in the actual extraction process to model the PSF parameters and the coordinate transformation from the reference catalogue to the cube as a function of wavelength.
The spectrum extraction is carried out in a layer-after-layer approach, starting at the central wavelength and then progressing alternatingly to the red and blue ends of the cube. In each layer, a sparse matrix is created, containing the model fluxes of one star (according to the current estimate of the PSF and its position) per column. Next to the stars, local background estimates are included in the matrix in order to account for the non-negligible surface brightness of unresolved stars or diffuse nebular emission. PampelMUSE allows for the definition of such background elements as square tiles with a user-selectable size. Via matrix inversion, the fluxes of all (stellar and background) sources are fitted simultaneously. Afterwards, all sources except those identified as isolated enough to model the PSF are subtracted and the parameters of the PSF and the coordinate transformation are refined. The new estimates are then used in another simultaneous flux fit. The procedure is iterated until convergence is reached on the source fluxes and the analysis of the next wavelength bin is started. Each new wavelength bin uses the resulting values of the previous bin as an initial guess.
After all wavelength bins are processed this way, a final PSF model is derived for all of the data cubes. To this aim, the values of the PSF parameters obtained in the individual wavelength bins are fitted with low-order polynomials. The object coordinates are also fitted with polynomials along the dispersion axis to reduce the effect of small random jumps between wavelength bins and thereby increase the S/N. The use of polynomials for this task is justified because ambient characteristics such as atmospheric refraction or the seeing should result in a smooth change of the PSF and the source coordinates with wavelength.
Notably, while the FWHM should always show a monotonic decrease with wavelength, which is the theoretically expected behaviour, we found that it instead increases where λ > ∼ 8000 Å. However, such a behavior is expected owing to contamination from second-order scattered light of bluer wavelengths, when using the extended mode of MUSE. Meanwhile, β did not vary strongly with wavelength. In the last step, the final spectra are extracted by traversing all bins of the cube once more, using the fitted estimates of the PSF and the object coordinates. This was done again by a simultaneous flux fit to all stars and background components. After convergence of the fitting process is reached, the stellar fluxes and the local background estimates are available as individual spectra for further analysis. The background spectra turned out to be very useful for the measurement of extremely faint diffuse gas emission, in particular at the wavelengths of Hα and Hβ where the nebular emission coincides with stellar absorption lines. Stars with successfully extracted spectra are referenced in what follows through the PampelMUSE input catalogue number per field as listed in column (1) of Table 6.

Fitting stellar spectra
Two major goals of the project were to demonstrate, as a proof of principle for crowded fields in NGC 300, that PSF-fitting IFS is capable of extracting spectra of individual stars with sufficient quality to derive trustworthy spectral type classifications, and to measure radial velocities, even with moderate to low signalto-noise ratios (S/N). To this end, we have fitted the extracted spectra both to an empirical library of stellar spectra, in what follows MIUSCAT, and to a grid of models computed with the Phoenix code (Husser et al. 2013), henceforth GLIB (Göttingen Library). MIUSCAT is an (unpublished) extension of the MILES library (Sánchez-Blázquez et al. 2006;Cenarro et al. 2007;Falcón-Barroso et al. 2011) that was kindly provided to us by Alexandre Vazdekis. It was created to reach up to the calcium triplet region and covers the gap between the blue and the red spectral ranges of MILES and CaT with data from the indo-U.S. library (Valdes et al. 2004) as a basis for the stellar population synthesis models from Vazdekis et al. (2012). Unlike other empirical libraries, MIUSCAT covers the entire free spectral range of MUSE in the extended mode.
The stellar spectra that could be extracted successfully by PampelMUSE were fitted to MIUSCAT by means of the ULySS code 3 , that was originally developed to study stellar populations of galaxies and star clusters, as well as atmospheric parameters of stars (Koleva et al. 2009). An advantage for our application is the fact that ULySS, as adapted from pPXF (Cappellari & Emsellem 2004), fits a spectrum as a linear combination of nonlinear components, multiplied by a polynomial continuum, thus helping to identify unresolved blended stars, as opposed to a single best guess for the spectral type in question. Especially for our application, where we are confronted with spectra of low S/N, the output of several (up to 10) library spectra with their respective weights supports a proper judgement of the quality of the fit and the identification of potential problems. As a practical disadvantage, the ULySS output requires a significant level of human interaction, i.e. a decision for each individual input spectrum in how far the linear combination of library spectra with different weights allows to make a good guess for the spectral type of the observed star. However, any less than plausible spectra can be ruled out immediately, e.g. foreground stars, or main sequence stars, the latter of which, at a distance of 1.9 Mpc, would be far below the detection limit of MUSE. In order to assist with this criterion, we searched the Simbad database for photometry of the MIUSCAT library stars and shifted their magnitudes to the most recent, cepheid-based NGC 300 distance modulus of (m-M) 0 =26.37 (Gieren et al. 2004 for comparison with the apparent magnitudes from the ANGST catalogue, neglecting extinction. An important aid in assessing the validity of the classification is the inspection of images to find any evidence of blending. To this end, we plotted for each object a reconstructed VRImap from the MUSE datacube as a post-stamp-like image of a size of 4×4 arcsec 2 , accompanied with an HST ACS image of the same region, colour coded from the filter combination F475W+F606W+F814W (Dalcanton et al. 2009), see Fig. 3. Together with a blend flag issued by PampelMUSE the probability and severity of blending effects was assessed and recorded as a quality flag.
As an alternative to ULySS fits to the MIUSCAT library, we used the technique of Husser et al. (2013), that was initially developed for globular cluster stars, to fit GLIB spectra to our objects. In this case, the result is a set of stellar parameters for a single locus in the HRD, with T eff , log g, and metallicity. Both the MIUSCAT and GLIB approaches also yielded measurements of the radial velocity.
For the final assessment of spectral type and radial velocity, we performed a visual comparison of the fits relative to the measured spectra, inspecting whether or not important absorption lines would be in accord with the fit and stand out from the noise, checked the resulting stellar parameters for plausibility, and listed the outcome with a set of quality flags to indicate (a) the quality and plausibility of the fits, based on the visual inspection of critical absorption lines and the photometry, (b) the agreement between the MIUSCAT and GLIB fits, any apparent effects of blending from nearby stars, and the plausibility of the measured radial velocities. In a final step, the most probable spectral type (or a range of spectral types) was determined, as well as the most probable radial velocity, and a global quality flag for the final result. The adopted outcome of the fits and a b c d e f visual inspection were recorded in individual log files and summarized in a catalogue file. In order to reduce the elements of subjective judgement on a steep learning curve, the whole procedure was actually exercised twice, and only the results of the final assessment were retained. The quality of spectral type classification is summarized in Table 5 and discussed in § 4.1. An excerpt of the catalogue is presented in Table 6.

Extracting non-stellar sources
As expected from the outset, visual inspection has shown indeed that our MUSE datacubes contain valuable information about objects that are not discovered as continuum point sources with the methods described above. This is particularly true for emission line objects like H ii regions, supernova remnants (SNR), supershells, diffuse ionised interstellar gas (DIG), planetary nebulae (PN), etc. Also faint background galaxies reveal their presence though redshifted emission lines that stand out from the foreground galaxy continuum. The search and classification of such objects was done by visually inspecting the emission line maps described in §3.1 with DS9. H ii regions were readily iden-tified on the basis of spatial extension and brightness. PNe were detected by blinking [O iii] vs. Hα, as well as He ii in order to find high excitation objects. SNR, superbubbles and supershells were distinguished from H ii regions by the technique of ionization parameter mapping of Pellegrini et al. (2012). Emission line point source objects were identified by measuring the FWHM of their images through a Gaussian fit and comparison with the PSF obtained for stars in the same datacube with PampelMUSE. By blinking Hα against continuum V, R, and I images the emission line point sources that were found to coincide with stars, yielded candidates for emission line stars as a valuable complement and cross-check for the detection of such objects made with PampelMUSE. All visually detected emission line objects are referenced with an ID given by a leading character for fields (a)...(i), followed by a serial detection number (regardless of type of object) as listed in Column (1) of Tables 7, 8, 9, 10, and 11. Subsequent analysis using the line fitting capability of the P3D visualization tool allowed to measure emission line fluxes and radial velocities for all kinds of emission line objects. As most of the objects exhibit several sufficiently bright emission lines suitable for the line fitting tool, we typically measured all sufficiently bright lines and determined a resulting Doppler shift for the radial velocity estimate from a flux-squared weighted average of all lines. The uncertainty estimate was obtained from the scatter of the contributing lines. Broad line profiles that are indicative of strong stellar winds were measured for point sources, which were thus identified as emission line stars. In extended sources, the [S ii]/Hα line ratio was used to discriminate SNR against H ii regions, which is particularly important for SNR candidates that are too faint for the ionization parameter mapping technique to be applied. A prerequisite for good flux measurements is an accurate subtraction of the background surface brightness that is composed of contributions from continuum light of faint unresolved stars and of diffuse or filamentary emission line flux arising from the DIG, ancient SNR shells, etc. Background estimates were obtained with P3D by defining an aperture for the object in question and a surrounding annulus for the background, where either strictly circular, or otherwise arbitrary user-defined geometries can be defined in order to accomodate complex surface brightness distributions. For the special case of recording DIG intensities, we used the unresolved background estimates as output from PampelMUSE to correct P3D flux measurements. Background galaxies were discovered by browsing the row stacked spectra available in the visualization tool of P3D and searching for emission features at unusual wavelengths. As a map is automatically displayed when the cursor moves through the suspicious wavelengths, a summed spectrum for the affected spaxels is readily created. This technique proved to be extremely efficient for data mining redshifted background galaxies.
In what follows, we describe how we have exploited the high level of sensitivity obtained with MUSE for the discovery of emission line objects like PNe, emission line stars, compact, normal, and giant H ii regions, SNR, superbubbles, giant shells, and DIG.
To this end, we first of all searched for point sources using DAOPHOT FIND as described in § 3.2 and extracted their spectra with PampelMUSE, similar to the procedure with stars. We also visually examined all of the fields in Hα and recorded extended and point sources down to very low contrast levels to create a provisional initial catalogue, without consulting the PampelMUSE catalogue to avoid any subjective bias in the detection process. By blinking against the [O iii] images we identified high excitation objects such as PN candidates, finding also objects that would not be bright enough to appear in the Hα image. The He ii maps allowed us to readily identify high excitation PNe and to discover a WR star (see below).
We then inspected each object from the initial catalogue one by one using the P3D visualization tool, and measured emission line fluxes with aperture spectrophotometry, assisted by an interactive background subtraction feature of P3D which allows one to define object and background apertures as standard regions of interest (ROI) of different geometries (circular, elliptical, rectangular), and a mouse-controlled editor to quickly add or remove spaxels. The P3D line fitting tool was also used to measure the central wavelength of emission lines and the corresponding lineof-sight radial velocity from the Doppler shift with respect to laboratory wavelengths.
The uncertainty of emission line fluxes was estimated using the equation by Gonzalez-Delgado et al. (1994), however adding an extra term to account for flat fielding and flux calibration uncertainties: where σ cont is the standard deviation of the continuum near the emission line, N is number of spectral bins used to measure the line, ∆ is the reciprocal dispersion in Å/bin, EW is the equivalent width of the line, and F l is the flux measured over the N spectral bins.
We also employed the PhAst tool (Mighell et al. 2012) to perform aperture photometry in the narrowband images to doublecheck the P3D flux measurements, and to estimate the FWHM for any point-like object that was found from the visual inspection. We finally merged the results from the two different approaches and classified the detected objects on the basis of emission line fluxes, line ratios, and the FWHM of point-like sources or otherwise the size of extended objects as further discussed in the following section.

Stars
Results: For the demonstration purpose of this paper, we selected field (i) as the best pointing in terms of seeing (FWHM=0.6"), as highlighted in Fig. 2, covering a fraction of the north-western spiral arm at a galactocentric distance of ≈1.5 kpc. We adopted both of the two procedures introduced in §3.2, i.e. the input of stellar centroid priors from high spatial resolution HST images, and alternatively from a search in the datacube using DAOPHOT FIND. PampelMUSE was run with both catalogues, resulting in a total of 3540 and 552 extracted spectra, respectively. It should be noted that the HST coverage is only 2/3 of the MUSE field (i). In order to reduce the number of poor quality spectra that would not allow conversion of the fitting procedure, we set a threshold of S/N=3 for the estimate computed by PampelMUSE and forwarded only spectra above the threshold to the ULySS code. From successful ULySS fits to the MIUSCAT library we obtained a total of 345 and 392 results for the HST and FIND input catalogues, respectively. For the same number of spectra, the fitting procedure was repeated using GLIB and the code from Husser et al. (2013). All of the spectra and corresponding fits were inspected visually, along with VRI images of the stars extracted from the datacube, as well as from HST, where available, for comparison and assessment of blending. Following the determination of spectral type, radial velocity v rad , and quality parameters as described in §3.2, we created the final catalogues for the two sets of spectra. An excerpt of the catalogue for field (i) with HST input is presented in Table 6. The complete catalogue for field (i) will be made publicly available through CDS. To illustrate the quality of spectra that can be obtained with MUSE, Fig. 4 presents three examples extracted from the stars in field (i) that are centered on the poststamp images (left panel: MUSE, right: HST). The spectrum extracted from the upper panel star numbered ID119 has a S/N=19 and is classified A1Ib with v rad =149±18 km s −1 , the one in the middle from ID260 has S/N=14 and is classified K3Iab with v rad =171±17 km s −1 , the one in the bottom panel from ID1379 has S/N=6 and is classified M6III with v rad =174±7 km s −1 . The plot range covers the full MUSE free spectral range from 460 nm to 930 nm. Black lines correspond to observed spectra, and the orange ones to the best MIUSCAT fits, that are for the most part almost indistinguishable from the measured data.   5 demonstrates how the PSF fitting technique works for an example of two stars in field (i) that are separated by 0.6" as determined from the corresponding HST image. The separation amounts to the FWHM of the MUSE PSF. The blend consists of a red star A to the north (ID1539, F606W=23.26), and a blue star B to the south (ID633, F606W=21.84). Deblending with PampelMUSE delivers the spectra as shown in the right panels: A is classified as an M4.5 subgiant, clearly identifiable through its prominent TiO absorption bands, wheras B is classified as A2Ia supergiant with relatively strong Balmer lines and a noticeable Paschen series. No obvious cross-talk is seen in any of the two deblended spectra, remarkably similar to the showcase for deblending of globular cluster stars shown in Fig. 2 of Husser et al. (2016). This capability is in stark contrast to the limitations of fiber-based spectroscopy, where an isolation criteria must be imposed, e.g. less then 20% contamination from neighboring stars (Massey et al. 2016). The HST input catalogue stars that have delivered useful spectra cover a magnitude range in the F606W filter of 20. . . 25 mag. Fig. 6 presents a case that is even more extreme: a faint blue star (A) (F606W=24.43) separated by 1.4" from K5III star (B) (ID682, F606W=22.49). Due to its faint magnitude and a blend with group (C) of five faint K stars, the case of star A is too difficult to allow for the extraction of a meaningful spectrum with PampelMUSE. The object, however, was discovered as an emission line point source (i101) close to the detection limit of the field (i) Hα map. Subsequent inspection using the P3D tool revealed that i101 is slightly offset to the north from group (C), thus most probably associated with the faint blue star B visible in the HST image. This is supported by the fact that the emission line is quite broad with FWHM(Hα)=6.5Å, i.e. not of nebular origin. The Hα flux was measured to 4.0 × 10 −18 erg/cm 2 /s, and the radial velocity as 144 km s −1 . The collapsed broad-band image from the MUSE cube shows merely a vague blue hue and illustrates the fact that for ground based observations only integral field spectroscopy opens a chance to discover such an object. In Fig. 7 we plot the S/N of extracted spectra of field (i) as a function of magnitude, broken down into three groups of different classification quality. The top panel illustrates the effect of the final estimate of spectral type classification with plausiblity flag PC2 . . . PC0, where apparently good fits are plotted with blue circles (PC2), marginal fits as green crosses (PC1), and uncertain cases as red circles (PC0). First of all, the latter ones are found mostly at S/N levels below 4. Secondly, marginally plausible fits show a similar behaviour, with only a handful of cases at brighter magnitudes and higher S/N. Thirdly, the distribution of the majority of good fits suggests completeness down to a F606W magnitude of 22.5. The lower panel shows the distribution broken down by effective temperature, with blue symbols corresponding to stars hotter than 5500 K, and red symbols to stars cooler than this temperature. This plot shows that in the range of F606W=22. . . 23 the cool stars tend to exhibit a higher S/N than the hot stars, reflecting the fact that the F606W magnitude is not a good measure for the red flux of cool giants, which show spectra with high equivalent width absorption lines, e.g. the calcium triplet, or in the case of M stars, very pronounced molecular bands, in a spectral region where MUSE is most sensitive, whereas the cut-off of MUSE in the extended mode at 460 nm limits the sensitivity for hot stars with regard to diagnostic lines in the blue. Moreover, the wavelength-dependence of seeing leads to a stronger susceptibility to blending in the blue than it does in the red, e.g. measured as FWHM=0.60" at 460 nm vs. 0.48" at 850 nm for field (i). One could also argue that dust extinction leads to a selection effect that is more important for blue stars, however it is beyond the scope of this paper to address this issue quantitatively.
We have tested the validity of our approach to spectral type classification in two ways. First of all, we checked whether the MUSE spectrophotometry correctly reproduces the ACS photometry. Artificial star tests as commonly used for CCD photometry were not found to be useful because PampelMuse works already with a catalogue of stars obtained at high angular resolution (HST) with the benefit of accurately pin-pointing stellar centroids, regardless of magnitude. We therefore used the information from that catalogue as a reference. To this end, we convolved our flux-calibrated MUSE spectra with the ACS F606W filter curve as provided at the Spanish Virtual Observatory website 4 . In Fig. 8 the resulting MUSE magnitudes (with a zeropoint of 36.5 mag) are plotted versus ACS magnitudes, whose errors amount to 0.01 mag for F606W= 22. . . 23, and a scattered distribution between 0.01 and 0.1 mag for fainter stars with F606W= 23. . . 25.5. The red plot symbols represent cool stars with col- Both datasets show a very similar distribution, with a tight correlation at bright magnitudes, and a rapid degradation towards large errors at faint magnitudes. Residuals against the 1:1 relation in magnitude bins of 0.2 mag are shown in the lower right part of the graph (with an offset of 20.5 mag). The error bars indicate the standard deviation in each bin, with values of order 0.1 mag for bright stars up to F606W=22.7, and an abrupt increase to 0.3 mag and larger for stars fainter than F606W=22.8. At this magnitude, the distribution begins to become asymmetrical, and there is the onset of a bias to larger magnitudes. With reference to Fig. 7, we interpret the branch towards bright magnitudes as the object photon shot noise dominated regime with a slope of -2, as expected. On the other hand, the branch towards fainter magnitudes must be source confusion limited, where the subtraction of blends and of the background of unresolved stars introduces errors at a level comparable with the poissonian noise of the object spectrum.
The second test was performed on the basis of seven template spectra from the MIUSCAT library, chosen such as to cover the relevant range of effective temperature and luminosity of the stars we would expect to discover in NGC300 (M1a-ab, K4III, K0III, G3III, F6Iab, A2Ia, B3Ib). These spectra were modulated with random noise resulting in S/N = 30, 20, 10, 7.5, 5, 4, 3, 2, 1, with 10 random realizations per S/N value, and a sample of in total 630 simulated spectra. The ULySS code was applied to each spectrum to recover the input spectral type from the noisy simulation. Fig. 9 shows the outcome of the exercise as the recovered T eff for each simulation versus the corresponding S/N value as provided as an output parameter by the ULySS code, color-coded for the different spectral types.
One can immediately see that at S/N ≥7.5 all of the spectra are perfectly well recovered (with the exception of two G3III outliers at S/N=7.5). Below that noise level, the reliability of the fit is strongly depending on the spectral type. For example, K0 and K4 subgiants show still agreement down to S/N=3, however with some outliers involving offsets of a few 100 K in T eff . The M1 supergiant spectra, owing to the pronounced molecular band features, are even recognized unequivocally down to S/N=2. This is in stark contrast to G3III spectra, that scatter across a range of almost 2000 K at S/N=4 and below. The F supergiant simulations are uniquely recovered again down to S/N=4, whereas the hot A and B supergiant spectra begin to scatter at S/N=4. The different robustness of spectral type identification against noise depends on the dominant strength of absorption features characteristic for different temperatures, e.g. the calcium triplet for cool stars, molecular bands for M stars, the Balmer and Paschen series for hotter stars, etc., and the spectral energy distribution with regard to those features. It is noteworthy that for the global input noise levels imposed by the simulations, ULySS returns S/N estimates that tend to be slightly lower than the input value, depending on spectral type. This points to the fact that a global S/N value is not a uniquely determining parameter for spectra over a broad free spectral range, such as with MUSE. For the sake of simplicity, we have not attempted to correct for these deviations. By combining our findings from photometry and the simulations, we conclude that the ULySS fits provide robust spectral type estimates for spectra with S/N ≥5, regardless of T eff . For a S/N below that level, hot and cool star show a different behaviour, chiefly in the sense that M and K stars still provide meaningful T eff estimates around S/N=3. . . 4, where A and B stars are already suffering significant scatter. From Table 4, we estimate a limiting magnitude for completeness at F606W≈22.5 for blue stars, and at F606W≈23.5 for red stars, in the sense that more than 50% of stars within a magnitude bin yield a reliable fit (nomenclature: N b 4 is the number of blue stars with SNR>4, N tot is the total number of stars within a magnitude bin, f indicates the useful fraction of the total per magnitude bin in %). These limits are dictated by the onset of crowding and the associated additional noise contributions for any given star -hence they are no sharp thresholds, but rather depend on the environment (i.e. amount of blending) in each individual case. Table 5 summarizes the distribution of spectral types as classified by the procedure explained above, comprising roughly 2/3 of the spectra extracted from the HST input catalogue, and complemented for the remaining 1/3 with DAOPHOT FIND detections. The total number of classifications is thus 517, of which 265 were assigned highly plausible (PC2, 51%), 160 marginal (PC1, 31%), and 92 uncertain (PC0, 18%). The distribution of stars (subgiants, giants, and supergiants) in terms of temperature reflects the stellar population of a spiral arm region, that was selected to not be dominated heavily by ongoing star formation, i.e. a minimal number of bright H ii regions, in order to maximize the detection rate of faint PNe. We find an appreciable number of O star candidates, including a significant fraction of hot emission line stars (11 out of 38, i.e. 29%), albeit the problem of our current approach using the MIUSCAT and GLIB libraries that do not allow for an accurate classification of hot stars due to the lack of O template stars in MIUSCAT, and the lack of non-LTE models for GLIB -a shortcoming that we are planning to resolve with an improved approach in the future. For the remaining spectral types there is good coverage such that we expect to having obtained a complete sample of spectral types B...M brighter than or equal to luminosity class III. 13 Fig. 10. Radial velocity uncertaintites ∆v rad plotted against signal-tonoise ratio. The sample is segregated into hot and cool stars as in Fig. 7. The insert shows a histogram of the ∆v rad distribution, the 50th percentile at 20 km s −1 .
The radial velocities determined by the MIUSCAT and GLIB fitting procedures were found in the majority of cases to be both plausible and in good agreement with each other, even for spectra with S/N as low as 3, where sometimes important diagnostic absorption lines such as e.g. the calcium triplet were only marginally discernible from the noise. The measure- ments merely failed in cases where excessive residuals from very strong nebular emission line background, that was too bright to be removed by PampelMUSE, prevented a reliable determination. We attribute this result to the robustness of our fitting procedures that are essentially based on a multitude of features over a free spectral range as large as an entire octave, rather than looking only into a few selected absorption lines. The values of v rad all seem to be realistic in that they are close to the systemic velocity of 144 km s −1 , with a mean of 169 km s −1 and a dispersion of 23 km s −1 in the case of field (i). Preliminary results of radial velocities from the remaining pointings indicate that the increase of median radial velocity per field with growing galactocentric distance in the range of 140...200 km s −1 is indeed sampling the rotation curve of NGC 300. This finding is further discussed in §4.7 and Fig. 25 below. Fig. 10 illustrates the scatter of formal radial velocity errors ∆v rad determined by ULySS versus S/N, plotted with red and blue symbols for hot and cool stars, respectively, as in Fig. 7. The cumulative histogram for the blue and red subsample shows that half of the radial velocity errors are below 20 km s −1 , and the 78th percentile at 40 km s −1 . The plot also illustrates that even faint red giant and supergiant spectra down to 24. . .25 mag with S/N≈3 or less have yielded useful radial velocity information.
An unexpected finding was the presence of spectra of cool stars that were not fitted at all, neither with MIUSCAT, nor with GLIB, despite a reasonable S/N. By comparison with spectra from the literature, e.g. van Loon et al. (2005), these stars turned out to be carbon stars that are as yet not covered by the libraries we use. We find 23 such stars in field (i) with a high level of confidence, rendering the C to M star ratio C/M=0.13, which is reasonably in line with the study of Hamren et al. (2015) in M31, who find C/M≈ 0.15 . . . 0.32 at [O/H]=8.6. In NGC 300, this compares to a metallicity for field (i) at a galactocentric distance of 1.5 kpc (R/R25=0.3) of [O/H]=8.5±0.1 dex, following Bresolin et al. (2009).
We have also discovered a number of emission line stars from the visual inspection after automatic processing for MIUS-CAT and GLIB fits. Their spectra are typically showing broad Hα and Hβ emission lines that are characteristic of hot stars with strong stellar winds. As our current fitting procedure is not sensitive enough to distinguish well enough between different classes of such stars, we have assigned, for the time being, merely a classification OBem. An alternative, more complete approach to discover and measure emission line stars and the discovery of a WR star is discussed in § 4.3.
For resolved stellar population studies one has to correct for contamination of the sample under study by foreground stars which are difficult to identify merely on the basis of photometry. The spectrum of star A in Fig. 26, ID 32224 in Table 6 with F606W=24.25, is classified a K5 dwarf, translating to an apparent R magnitude of 33.3 for the distance modulus of NGC 300, which is obviously by far too faint to be observable. The radial velocities determined from the MIUSCAT and GLIB fits are 66 ± 6 km s −1 and 69 ± 3 km s −1 , respectively, i.e. almost identical, each with a very small uncertainty. This velocity is far away from the systemic velocity of NGC 300 and therefore, together with the photometric evidence, indicative of a foreground star. In this particular case, the situation is somewhat more complicated as it turns out that in the HST image star A is resolved into two distinct stars. From the ULySS fit, we conclude that there is a weak blend from an M6III star, which explains the less than perfect fit of the spectrum. Nevertheless, the facts remain that for the main sequence K spectrum component, such a star is not visible at the distance of NGC 300, and that the well-constrained radial A&A proofs: manuscript no. NGC300_MUSE_2column_arXiv Notes. Excerpt of catalogue for field (i) with spectral type classification for spectra extracted with PampelMUSE, covering ranges of brighter and fainter F606W magnitudes (full table available online). Column 1, ANGST catalogue ID; Col. 2, F606W magnitude; Col. 3, right ascension (J2000); Col. 4, declination (J2000); Col. 5, signal-to-noise ratio estimate from PampelMUSE; Col. 6, flag "blending": 2=not obvious, 1=minor, 0=significant; Col. 7, effective temperature of best fit [K]; Col. 8, gravity of best fit; Col. 9, metallicity of best fit; Col 10, radial velocity [km s −1 ]; Col. 11, radial velocity uncertainty [km s −1 ]; Col 12, spectral type of best fit, or probable range of spectral type for ambiguous cases; Col. 13: flag "plausibility of fit": 2=very plausible, 1=marginal, 0=uncertain. "999" in Col. 10 and 11 indicate uncertain results. Hot stars assigned OB, T eff is set to 30000. PER: peculiar red star visible in HST image, not detected in MUSE cube.
velocity is incompatible with stars within this galaxy. Using the new Besancon Galaxy model (Czekaj et al. 2014) we have estimated the number of stars that one would expect in the direction towards NGC 300 down to V = 20 to 0.5 stars per 1×1 arcmin 2 . Extrapolation to V=24.25 predicts approximately 4500 star per square degree, i.e. 1.25 stars per MUSE pointing, which seems to be in accord with our observation.
Discussion: Our first results from a complete analysis of field (i) have shown that the application of PSF-fitting crowded field 3D spectroscopy as previously exercised in globular clusters could successfully be established in the more challenging case of nearby galaxies, where, owing to a much larger distance, we can no longer expect to reach down to the main sequence (except for the most massive O stars). However, we have demonstrated that we are able to sample the population of giants, subgiants, and supergiants, from which we obtain spectra with sufficient S/N to make a spectral type classification, as well as to measure radial velocities with accuracies around 20 km s −1 . We have also demonstrated that the deblending technique works well, even for stars as faint as m F606W ≈ 23, provided we can rely on an input catalogue of HST stars.
We note that, thanks to the high efficiency of MUSE, the spectra of blue, yellow, and red supergiants in our sample have good enough S/N, despite the relatively short total exposure time of 1.5 h of the current data set, to allow for quantitative spectroscopy along the lines of, e.g. Bresolin et al. (2002), Kudritzki et al. (2008), considering numerous absorption lines of the Balmer and Paschen series of hydrogen, of Mg, Fe, Na, Ca, etc., that are immediately apparent in the spectra. We have one Article number, page 13 of 37 A&A proofs: manuscript no. NGC300_MUSE_2column_arXiv star in common with the Bresolin et al. (2002) sample that is suitable to make a direct comparison: their object B11, classified as A5 supergiant (marked in Fig. 2). This star was observed with FORS1 at the VLT with approximately 5 Å resolution, a total exposure time of 3.75 h, and seeing conditions ranging from 0.4" to 0.7". The MUSE parameters for comparison are: 2.5 Å resolution, 1.5 h exposure time, 1.2" seeing. The star is listed as ID25 from a preliminary analysis of field (b). The MIUSCATbased classification of this star is A1Iab...A5II, in agreement with Bresolin et al. (2002), despite a lower S/N due to poorer observing conditions and a shorter exposure time. Fig. 11 shows our spectrum of this star, along with the best ULySS fit (with an offset for clarity). Unfortunately, the overlap between the FORS and MUSE spectra is restricted to a narrow interval in the Hβ region of 4600. . . 4900 Å. However, the comparison is good enough to qualitatively demonstrate that MUSE performs at least as well as FORS, considering the less favourable observing conditions for the MUSE spectrum (exposure time, seeing), and despite of the fact that in this wavelength region the efficiency curve of MUSE experiences a steep drop towards the blue. The importance of the statistics of red and blue supergiants was stressed by Massey et al. (2017). From a review on massive stars by Massey (2003), and a more recent paper by Massey et al. (2016) describing a survey in M31 and M33, the amount of work necessary to safely detect and classify massive stars is quite apparent. However, the exercise to create a census of massive stars is a prerequisite to understand stellar evolution at the high mass end. The method is based, firstly, on broadband and narrowband photometry, to identify candidate stars. As photometry in the optical merely samples the Rayleigh-Jeans tail of the spectral en-ergy distribution of hot stars, it is, secondly, necessary to perform follow-up spectroscopy, which is a tedious task with conventional multi-object spectrographs. As we have shown with a proof-of-principle in field (i), such a two-step procedure is no longer necessary with IFS. MUSE allows, for the first time over a reasonable FoV, to combine the different survey strategies into one: imaging and spectroscopy combined in a single datacube under identical observing conditions. Concerning the imaging capabilities, it is worthwhile noting that the discovery of evolved massive stars with heavy mass loss and fast stellar winds such as LBVs or WR stars is much facilitated by narrowband images that are readily extracted from a datacube at the wavelengths of interest, e.g. Hα, or He ii λ4686, see for comparison Massey et al. (2015). Although this capability was already demonstrated with PMAS by Relaño et al. (2010) and Monreal-Ibero et al. (2011), who found WR stars in M33, it is only now with the advent of MUSE that we are able to analyse such data over a wide FoV.
The argument also holds for stars with spectral types other than OB. Hamren et al. (2016), for example, have studied carbon stars in the satellites and halo of M31 using 10 years worth of data from the SPLASH survey (Guhathakurta et al. 2006;Tollerud et al. 2012). The ratio (C/M) of carbon-rich to oxygenrich AGB stars has been used to study the evolution of AGB stars, concerning observational constraints on the thermally pulsating AGB phase, dredge-up processes, opacities, and mass loss during this phase. C/M ratios have also been used to study the galactic environment in which the stars have formed. The data of Hamren et al. (2016) include 14143 stellar spectra taken with DEIMOS at the Keck-II 10m telescope. The spectra were obtained from 151 individual DEIMOS masks for 60 separate fields, with a typical exposure time of 3600 s per mask. The observations were targeting the environment of M31, i.e. dwarf spheroidals, dwarf ellipticals, the smooth virialized halo, halo substructure, and M32. From this entire dataset, 41 unambiguous carbon stars were identified. It is interesting to note that spectroscopy is apparently essential for obtaining a reliable identification. We take as an example the dwarf elliptical galaxy NGC 147: based on NIR photometry, Davidge (2005) found 65 carbon stars, Sohn et al. (2006) reported 91 carbon stars using the same technique, while Hamren et al. (2016) confirmed merely 12 stars on the basis of photometry and spectroscopy in the visual. For comparison, our 1.5 h exposure in field (i) of NGC 300 has yielded immediately 23 detections of carbon stars. While this is not really a fair comparison to the sample of Hamren et al. (2016) owing to effects of metallicity and the size of the underlying stellar population, we note that a similar study in the disk of M31 by Hamren et al. (2015), which would be more comparable to NGC 300, has yielded a total of 103 carbon star identifications from a sample of 10,619 spectra, again taken with DEIMOS and the circumstances described above.
Surveys of stars in external galaxies like the ones cited above must be corrected for foreground contaminant stars. This is usually accomplished spectroscopically, e.g. by constraining spectral type and luminosity class, and thereafter comparison with photometry, but also by comparing the measured radial velocity with the systemic velocity of the galaxy under study. We emphasize that through our automated procedure described in §3.3, spectral type and radial velocity are directly reported as part of the process, so such contaminant stars are immediately identified without further effort. The example in §4.1, Fig. 26, illustrates how the foreground star detection is an integral part of the data analysis pipeline. As shown in Fig.7 and Table 5, the extracted spectra for field (i) span a range of S/N with a total of 265 that are flagged "good". Of those latter spectra several tens have a S/N>10. We have already referred to the blue supergiant B11 whose spectrum in Fig. 11 is a less than ideal example due to the mediocre seeing of the observation in field (b). The S/N estimate for this spectrum is 17, decreasing somewhat towards the blue. Despite the modest S/N, a number of absorption lines can be seen. A higher quality spectrum is presented for the A1 supergiant ID151 of field (i) in Fig. 12 that exhibits a S/N of 22 -despite the fact that this star is about a magnitude fainter than B11, illustrating the importance of high image quality for this kind of observation. In contrast to the blue supergiants with rel-atively few identifiable lines, our example of the late K to early M supergiant ID301 in field (i) with m F606W =21.56 and S/N=14 shows a multitude of absorption lines, including TiO bands, that are reproduced by the MIUSCAT fit (Fig. 13).
It is clear that with longer exposure times, the quality of such spectra will be improved. We stress again, however, that excellent seeing is a prerequisite. Insofar the installation of the adaptive optics (AO) facility at VLT-UT4 and the GALACSI module in front of MUSE (Stuik et al. 2006) will even further improve the deblending performance of our technique in crowded fields. For the time being, we have focused our effort on the automatic processing of hundreds of spectra with a pipeline for the purpose of spectral classification and the determination of radial velocities. The results from a provisional analysis for all fields (a). . . (i) are discussed later in §4.7.
Obviously, the comparison with model atmospheres and the determination of abundances is the next logical step. However, we have as yet only begun to develop the technique of crowded field 3D spectroscopy at 2 Mpc distance. In the first instance, we required to visually inspect all spectra one by one in order to gain confidence in the procedure. This turned out to be a very time-consuming procedure. For the future, the goal is to develop a fully automated, robust pipeline with very little, if any, human interaction. Also, as a weakness of the current scheme, the MIUSCAT and GLIB libraries do not support well the analysis of hot stars. This is an issue that we are planning to address in the near future.
In summary, we believe that the first analysis of spectra in field (i) presented in this paper has validated crowded field 3D spectroscopy with MUSE as a powerful tool for quantitative spectroscopy of individual stars in nearby galaxies.

Planetary nebulae and compact H ii regions
Results: The analysis of our datacubes for the discovery and measurement of point source and extended emission line objects began with the set of narrowband images for the most important wavelengths as listed in Table 3. An example for the wavelength of Hα is shown in Fig. 14 as a montage of fields (a). . . (j) over the Hα image from Bresolin et al. (2009), obtained with the ESO/MPI 2.2m WFI. Not surprisingly, the MUSE data obtained at the VLT show much finer detail and go significantly deeper than the WFI image, which can be appreciated well e.g. from the giant H ii region, marked De100, that is cut between WFI and MUSE at the eastern edge of field (a), or the variety of giant shells and filaments in fields (a). . . (j) that can be found nowhere at this level of detail in the WFI image.
One of the immediate scientific goals motivating our NGC 300 observations was to measure the [O iii] planetary nebula luminosity function (PNLF) down to unprecedented faint magnitudes in order to infer its diagnostic value for stellar population studies. Jacoby (1989) introduced the m 5007 magnitude, based on the line flux of [O iii] at 5007 Å that is conventionally used to represent the PNLF: m 5007 = −2.5 log F 5007 − 13.74 (2) The bright end of the PNLF was initially proposed as a standard candle for extragalactic distance determinations by Jacoby (1989) and collaborators, and subsequently utilized to measure distances to more than 60 galaxies; see review by Mendez (2016). However, it was also suggested that the exact shape of the PNLF at fainter magnitudes than the bright cutoff could potentially unravel useful information on intermediate age stars of (a) 0.7'' De42 S16 S10 S14 S8  Deharveng et al. (1988) are marked with red circles. Yellow circles denote the supernova remnants recorded by Blair & Long (1997).
the underlying stellar population (Ciardullo 2010), and eventually help to understand the apparent invariance of the cutoff in the first place (Kwitter et al. 2014).
Beyond the classical narrowband filter techniques that have been employed at 4m class telescopes to typically sample the brightest ∼1.5 mag of the PNLF, which is sufficient for distance determinations, MUSE offers the advantage of providing filter maps with extremely narrow width at many wavelengths in parallel, combined with the VLT 8.2m light collecting power and excellent seeing, thus extremely high sensitivity. Moreover, we were able to apply the PampelMUSE PSF-fitting technique to extract PN spectra, deblended them from their crowded stellar and nebular environments, and then measure accurate emission line fluxes and radial velocities.
In order to assess the sensitivity that we can reach with MUSE, we created sequences of artificial PN point sources in [O iii] 5007 Å with decreasing brightness that were inserted in our observed datacubes (Fig. 15). We were able to demonstrate for all fields (a),. . . ,(j), representing different seeing conditions of 0.6 ". . . 1.2", that we indeed detect these objects using FIND, and in the best case reach completeness down to m 5007 =28.0, which is about 6 magnitudes fainter than the PNLF cutoff. The diagram in Fig. 15 shows that the completeness limit, as expected, depends strongly on the image quality, ranging from the best case at 0.6" and m 5007 =28.0 for field (i) to 1.2" and m 5007 =26.0 for field (b). The faintest PNe that we detected in the experiment have m 5007 =30.0.
After having convinced ourselves that the data quality of our datacubes would be sufficient to search for faint PNe, we performed the search and classification procedure as described above.  (1) aligned with M star; (2) bright, variable nebular background; (3) high extinction; (4) visual detection, not detected by FIND, also flagged as ":" in Col. 1; (5) retrograde velocity; (6) edge of field; (7) blend with compact H ii region d97, located 0.5 arcsec to the SW; (8) blend with bright star at same position.
from the analysis is presented in Table 7 where we list a total of 45 PNe, out of which 34 are new discoveries, and 11 were already known from the work of Soffner et al. (1996), and Peña et al. (2012). Nine extremely faint objects (m 5007 ≈ 27 . . . 29) are considered PN candidates on the basis that we could measure the [O iii] line at the expected wavelength, however other lines fell below the detection limit. All of the known objects from the literature were easily rediscovered in our data, whereas almost all of the new detections are as expected significantly fainter than the objects found with the NTT by Soffner et al. (1996)), and VLT-FORS by Peña et al. (2012). Our faintest PN candidate e08 was measured to have m 5007 =29.67. The analysis of results also lead to the unexpected observation that the visual detection method is surprisingly sensitive: all of the FIND detections were recovered, whereas a total of 19 PN where found only by visual examination, without detections from FIND -which is in line with practical experience from the early PNLF distance determination work (R. Ciardullo, priv. communication). Table 7 also lists the radial velocity estimates that were derived from the emission line fits with P3D. The overall distri- bution matches well the one of stars, thus giving confidence to indeed observing the rotation of the galaxy -however with one exception: PN i84 exhibits a velocity of 108 ± 13 km s −1 which is far off the median in field (i), and also blueshifted from the systemic velocity of 144 km s −1 , thus either indicating an extreme retrograde motion and the fact that the central star of this PN may be a runaway or a halo star, or that the object does not belong to NGC 300.
Apart from the fact that cHII present a challenge to identify a clean sample of PNe, they are interesting objects on their own right. We have therefore compiled a listing in Table 8 with essentially the same entries as for the PNe in Table 7. We note that we have, again, recovered all of the objects found by Soffner et al. (1996), and Peña et al. (2012). However, upon inspection of the spectra we concluded that as many as 7 of their PN candidates must be classified cHII, as they show a [S ii] line ratio > 1 as discussed above. Similar to the PNe, the cHII kinematics in the seven fields (a),...,(j) match well the stellar kinematics in those fields. Discussion: One of the major objectives of our MUSE observations of NGC 300 was the study of the PNLF down to very faint magnitudes. We have detected a total of 45 objects, reaching seeing-dependent completeness limits between m 5007 =28.0 and m 5007 =26.0 for a seeing FWHM of 0.6" and 1.2", respectively. PN candidates that were identified in field (e) reach magnitudes even fainter than m 5007 ≈ 29.0. The latter may seem to be surprising in view of the fact that pointing (e) is not complete in the sense that only one third of the nominal exposure time of 1.5 h could initially be secured. However, with 9 firm PN detections and 4 more PN candidates, field (e) has the highest PN number density of all pointings. Even though our completeness limit estimates are supported by this result, another surprising finding is, that pointings (i) and (j) have yielded significantly fewer detections, although the excellent to good image quality would have suggested otherwise.
We have compared our detections with results from the literature, namely Soffner et al. (1996) who used narrowband imaging at the 3.5m ESO-NTT, and Peña et al. (2012) with data from narrowband imaging and spectroscopy using FORS at the VLT. Firstly, the MUSE observations have recovered all of the previous detections. We have classified, however, from the 18 objects in common with Peña et al. (2012) as many as 5, i.e. almost one third, as cHII candidates. The same is true for 2 out of 8 PN candidates in common with Soffner et al. (1996). As these objects are all fainter than m 5007 =24.5, our reclassification should however have no major impact on the bright cutoff of the PNLF constructed from their samples. Our new detections that have no counterpart in the literature are typically fainter than m 5007 =27.
We have also compared our m 5007 magnitudes with the samples from Soffner et al. (1996) and Peña et al. (2012) and plot the outcome in Fig. 16. We have 17 objects in common, of which 5 are cHII. The results obtained by Soffner et al. (1996) are within the errors in agreement with our data, except for the very faintest magnitudes. However, there is a systematic offset with a median of ∼0.67 mag between the FORS (Peña et al. 2012) and our MUSE data. The radial velocities determined for our sample of PNe are discussed in §4.7 below.
In an attempt to understand the unequal detection rate over our pointings that cover the nucleus, parts of the central spiral arm to the NW, versus parts of the leading and trailing interarm regions around this spiral arm, we have complemented our detections with the PN candidates from Peña et al. (2012) and plotted all objects for orientation over an Hα map of NGC 300, see Fig. 25. The coloured plot symbols refer to MUSE planetaries, while the filled white circles refer to the Peña et al. (2012) objects.
Although the latter do not reach the same faint magnitudes as our sample, there seems to be a hint of the surface density distribution to favour a higher concentration of PNe near the nucleus and along the spiral arms. Owing to the lack of deep coverage at the level of our MUSE data, this observation is not entirely compelling, however supported by the sparse population in field (j) and the interarm extension to the south-east (SE). Whether the latter is a generic feature of this region, or perhaps an effect of extinction owing to a dust lane extending along the interarm region, remains to be clarified. At any rate, a concentration along star forming regions of a spiral arm would not at all be expected from a theoretical point of view. According to Renzini & Buzzoni (1986), and more recently Buzzoni et al. (2006), the number A&A proofs: manuscript no. NGC300_MUSE_2column_arXiv It remains unclear, however, why there is an extreme number of detections in field (e). More than half of the objects extend from the leading edge of the spiral arm into the interarm region to the NW. Two of the PNe in this area exhibit extremely low radial velocities (see below). The PN surface density pattern does not seem to continue into the adjacent field (d). Perhaps we are observing just a spurious density peak that would average out on larger scales. As the exposures in (d) and (e) are as yet incomplete, we consider the current census as preliminary.
We note that this is not the first work to discover extragalactic PNe in MUSE datacubes. Recently, Kreckel et al. (2017) published a paper on the PNLF of the face-on grand design spiral galaxy NGC 628. The faintest PNe in their sample have a magnitude of ∼27.9, which is consistent with PNe in our field (e), where the exposure time was 30 min and the seeing 0.75". For NGC 628, Kreckel et al. (2017) estimate PN surface densities of 1.6. . . 3.1 PN/kpc 2 for two regions with different galactocentric distances. These numbers are significantly lower than the preliminary values that we find in our fields: 23.2 and 33.5 PN/kpc 2 for fields (a) and (e), and 16.6 PN/kpc 2 as an average of the total of 7 fields, respectively, -while NGC 628 is ∼3 mag more luminous in the H band than NGC 300, suggesting a larger parent stellar population for the observable PNe.

Emission line stars
Results: The visual examination of the Hα maps yielded point source candidates whose spectra are not compatible with PNe or cHII, however showing hints of a continuum. Blinking the Hα map with the corresponding VRI image sometimes suggested coincidence with a star. Only a few of these stars could be associated with an object from the catalogue of PampelMUSEextracted spectra, probably because the emission lines in Hα and Hβ, filling in the absorption lines at the corresponding wavelengths, have prevented ULySS from converging to a reasonable fit, as emission line stars are not part of our empiricial library, nor would be supported by GLIB. In such cases, we have employed P3D to perform aperture spectrophotometry on the stars in question, with the caveat of potentially subtle background subtraction issues at the wavelengths of emission lines. However, it was possible to extract spectra that in most cases presented broad, non-Gaussian emission line profiles with extended wings, distinct from the narrow nebular background emission, that are characteristic for hot massive stars with strong stellar winds, see e.g. Klein & Castor (1978), Leitherer (1988), Lamers & Leitherer (1993), Puls et al. (1996), Kudritzki et al. (1999). An example is shown in the bottom panel (i) of Fig. 20. For this object, a Gaussian fit to the Hα emission line profile yields a FWHM of 6.6 Å, however with a significant residual owing to extended wings to the red and the blue (the MUSE line-spreadfunction has a FWHM of 2.4. . . 3Å). The linewidth and the extended wings point to stellar wind velocities of several hundreds up to 1000 km s −1 , typical for hot massive stars.
Discussion: In many cases we could immediately associate a blue star from our VRI maps with the emission line object. However, there were also cases where we could not find an obvious bright optical counterpart, which could possibly be attributed to strong dust extinction: given that there is an immediately visible distinct dust lane at the trailing edge of the north-western (NW) spiral arm in NGC 300, and another one at the leading edge of this spiral arm, and possibily significant intrinsic dust extinction in star forming regions as well, the broad line emission may act as a beacon to detect the associated massive stars even in the event of high extinction. In some cases we also found objects with a narrow line width. At the distance of NGC 300 it is impossible to disentangle compact nebulosities from intrinsic stellar emission on the basis of the MUSE data. Again, we visually inspected all of the spectra and decided which objects would  Table 9.
Direct comparison with VRI colour maps has revealed emission line objects that most likely are associated with red giants rather than with hot massive stars as discussed above. As an example, i54 shows a broad Hα line on top of a pronounced M star spectrum. Unfortunately, the object was found in a region where there is no coverage with HST images. Although we cannot rule out a chance alignment of the AGB star with a hot star that would be responsible for the emission, there is no indication from the MIUSCAT fit that this might be the case. There is, however, a striking similarity to the spectrum of the prototypical object R Aquarii, which is a well-studied case of a symbiotic M star, showing a hot companion, Roche lobe overflow, accretion onto the hot star, and jet activity (Zhao-Geisler et al. (2012)). The latter authors estimate the distance to R Aqr as 250 pc. Schmid et al. (2017) used observations with HST and SPHERE ZIMPOL at the VLT to measure an Hα flux of 4.6 × 10 −11 erg/cm 2 /s. Pushing this object to the distance of NGC 300, this would translate to a flux of roughly 10 −18 erg/cm 2 /s, i.e. just detectable with MUSE. We therefore flagged red giant emission line stars as symbiotic star candidates.
A special case is the Wolf-Rayet star (WR) that was readily detected from the He ii narrowband image created for field (c). Inspection of the PampelMUSE-extracted spectrum confirmed that this object is indeed a hot star with the characteristic blue WR bump (see Fig. 17).
Article number, page 21 of 37 A&A proofs: manuscript no. NGC300_MUSE_2column_arXiv  [Å], "w" indicates extended line wings; Col. 11, quality: confirmed, "?" uncertain, "-" no hint of a star, probably other type of emission line point source; Col. 11, apparent colour of star, "-f": faint, "-s": symbiotic star candidate, "-c": carbon star.  Fig. 14 bears some similarity with the 30 Dor region in the LMC. There is quite a number of prominent shells of diameters up to 100 pc visible, e.g. quite prominently in field (d). Besides the point sources that were catalogued as cHII, we also found many compact, only slightly extended H ii regions. As one can appreciate from Fig. 14, the earlier compilation of H ii regions from Deharveng et al. (1988) mainly comprise bright objects, while our new data reveal a wealth of newly discovered, fainter objects.
We have made an attempt to distinguish between classical photoionized, optically thick H ii regions, optically thin H ii regions, and diffuse ISM using the technique of ionization parameter mapping (IPM) that is based on the extent of Strömgren spheres in different ions (Koeppen 1979a,b), as excercised extensively by Pellegrini et al. (2012) in the Magellanic Clouds, and also recently applied for NGC 628 with MUSE data by Kreckel et al. (2016). Fig. 18 shows two examples of optically thick and thin H ii regions, catalogued as i97 and i81, respectively. i97 has a compact appearance with a peak Hα surface brightness of 9.72 × 10 −16 erg cm −2 s −1 arcsec −2 (left). The IPM greyscale image, with low values of [S ii]/[O iii] in black and high values in white, shows the typical picture of an optically thick region with a high excitation core and an enhanced line ratio as a ring, as presented by Pellegrini et al. (2012) in their Fig. 3 for DEM S38 in the LMC. Conversely, i81 has a peak Hα surface brightness of 1.54 × 10 −16 erg cm −2 s −1 arcsec −2 and shows qualitatively the same appearance as DEM S159 (their Fig. 4.) with areas of both high and low optical depths, namely a distinct rim to the E, complemented by an open lobe to the W, and again highly ionized gas in the central part.
We inspected the IPM images for all fields (a). . .(j) and marked the discovered H ii regions accordingly in the catalogue. The [S ii]/[O iii] line ratio map is also sensitive to reveal other interesting features. Besides the brighter PNe appearing as black blotches that were already discovered otherwise as described above, the map brings out high excitation emission, e.g. a peculiar structure seen as extending radially from the rim of i97 (arrow in Fig. 18, top row, middle). It it not visible in any other emission line rather than [O III] and is roughly, but not quite aligned with a blue star. The nature of this strange feature is unclear.
The catalogue of H ii regions is given in Table 10. Emission line fluxes, radial velocities, positions, as well as sizes and morphology were determined as described in Section 3.4.
Also the compact H ii region a32 that is embedded in a high surface brightness area of the H ii region a65, impossible to discern on the basis of an Hα image, was discovered this way. Altogether, we discovered 114 H ii regions, measured their fluxes, sizes, radial velocities, and recorded any luminous blue stars visible with the nebulae that could be responsible for photoionization.
For the ISM, Chu (2008) distinguishes interstellar bubbles with diameters up to 30 pc that are powered by the stellar winds of individual massive stars, superbubbles with sizes of ∼100 pc, dynamical ages of ∼ 10 6 years that require only one episode of star formation, and supergiant shells that have sizes of ∼1 kpc, dynamical ages of ∼ 10 7 years and require multiple episodes of star formation. Furthermore, there are SNR, that can be distinguished from the latter class of objects on the basis of the line intensity ratio F([S II)]/F(Hα) > 0.4 (Mathewson & Clarke 1973). A detailed investigation of the nature of all of these nebulae is beyond the scope of this paper, so we have restricted ourselves to list the objects by characterizing them as compact H ii regions, classical optically thick/thin H ii regions, or giant H ii regions with sizes larger than 100 pc, following the scheme of Franco et al. (2003), furthermore as SNR candidates, shells, or diffuse/filamentary ISM. SNR and the DIG are discussed further below.
Discussion: The most comprehensive catalogue of H ii regions in NGC 300 is still the one of Deharveng et al. (1988). It was complemented with more recent data from Bresolin et al. (2009), whose 28 H ii regions, however, do not overlap with our MUSE observations. The red circles in Fig. 14 indicate 8 objects from Deharveng et al. (1988), that coincide partially or fully with fields (a). . . (j). It is not surprising that with a total of 61 H ii regions our MUSE observations go significantly deeper than their photographic data from 30 years ago. Since we were able to sample low emission line intensities down to levels below 10 −17 erg cm −2 s −1 arcsec −2 , we could use the IPM technique to categorize H ii regions as optically thick or thin, and to identify superbubbles and giant shells. The fluxes measured by Deharveng et al. (1988) are in accord with our measurements within the errors, as far as their objects are entirely contained in our fields. As a convenient feature of the MUSE datacubes, we are able to immediately identify blue stars that are likely to be the ionizing sources of the H ii regions. As a caveat, it is important to realize that the spatial resolution of MUSE per se is no guarantee for having resolved individual stars. With a sampling of 1.8 pc per spaxel, even in the event of seeing at the Nyquist limit, the image of the PSF projected into NGC 300 corresponds to 3.6 pc, which is a size too large to safely resolve e.g. OB associations. Therefore, the use of HST images and catalogue stars as input to PampelMUSE is an important prerequisite for the identification and possibly deblending of the ionizing stars. Even with HST at a sampling of 0.048"/pixel, i.e. an image scale of 0.44 pc/pixel, very compact clusters or binaries may not be resolved.
The decomposition of blended spectra with our library fitting technique may provide a solution to this problem. In order to make this work, we are planning to increase the scope of our stellar library with more O and B stars. The goal of this exercise would be to determine the energy budget of H ii regions and obtained a unified picture concerning leaking Lyman continuum radiation and the ionization of the DIG (see below). The wellknown fact that photometry in the optical cannot constrain the number of ionizing photons of hot stars, because the Rayleigh-Jeans tail of the spectral energy distribution is insensitive to the effective temperature (Hummer et al. 1988), was realized in the outcome of an experiment by Niederhofer et al. (2016), that consisted in the attempt of predicting the content of O stars in clusters within H ii regions from isochrone fitting to HST photometry. In a follow-up project, similar to the work of, e.g. Bresolin et al. (2002), or Kudritzki et al. (2008), we are planning to experiment with MUSE spectra to constrain the luminosity of the ionizing stars. In addition, we are also planning to further investigate the chemical (oxygen) content of our data in order to better constrain the distribution of the metal within this galaxy using different metallicity indicators as well combining our data with the literature ones.

Supernova remnants
Results: We inspected the [O iii] maps for extended emission to find SNR candidates and subsequently used the spectra to select only those objects that show evidence for shock excitation on the basis of the canonical line ratio criterion [S ii]/Hα > 0.4 (Mathewson & Clarke 1973). We found prominent, high surface brightness examples, e.g. Fig. 21, as well as low surface brightness objects that barely stand out from the noise. In our dataset, we have reached limiting surface brightness levels of 2×10 −17 erg cm −2 s −1 arcsec −2 in [O III]. In total, we detected 38 SNR and measured their fluxes, sizes, and radial velocities. Nebulae that have a similar appearance as SNR and spectra pointing to shock excitation rather than photoionization, however formally failing the [S II]/Hα line ratio criterion, were classified as shells to distinguish them from H ii regions. Discussion: SNR in NGC 300 were first recorded by Dodorico et al. (1980) who made photographic observations with narrowband filters at the UK Schmidt and found 8 SNR candidates. Blair & Long (1997) performed a systematic narrowband filter search in Hα, [S II], and continuum using the 2.5m telescope at Las Campanas Observatory with a CCD imager and identified 28 SNR candidates in NGC 300. With follow-up spectroscopy at the same telescope they investigated in detail the empirical classification parameter of the [SII]/Hα line ratio that is required to be larger than or equal to 0.4 to qualify for a SNR. The extra effort of spectroscopy was necessary to eliminate the additional flux contribution from the [N II] lines adjacent to Hα that would normally fall within the transmission band of a narrowband filter. Their survey reached a limiting Hα surface brightness of 1×10 −16 erg cm −2 s −1 arcsec −2 . For comparison, the Hα surface brightness limit of our data is below 1×10 −17 erg/cm2/s/arcsec 2 , i.e. an order of magnitude fainter, see [S II] map in Fig. 23 for comparison. The total number of SNR detected in our MUSE datacubes is 38, with diameters in the range of 20. . .110 pc. The size distribution of SNR in M33, a galaxy similar to NGC 300, was modeled and inferred observationally by Asvarov (2014). Assuming a similar distribution for NGC 300, we expect that the median of the distribution is found at SNR diameters near 40 pc, with maximum diameters between 150 and 200 pc, depending on different values of the filling factor of the warm phase of the ISM. However, as pointed out by Asvarov (2014), the predicted size distribution is subject to observational selection effects, diminishing the detection rate for extended low surface brightness SNR in particular. We found one example of an extremely nitrogen-rich SNR (a71) similar to the one reported for the galactic SNR Pup A by Dopita et al. (1977), or several such objects in M81, M101, and NGC 6946 Matonick & Fesen (1997).
The dichotomy of filter spectrophotometry and the associated problem of [N II] contamination of Hα fluxes vs. slit spectroscopy, as discussed by the latter authors, is lifted when using IFS, which is an advantage in terms of efficiency, because no extra observations are necessary. Moreover, the fact that all derived quantities were obtained under identical observing conditions within a given datacube is an asset.
From the observational point of view in the optical, there exists no sharp distinction between SNR, superbubbles, or giant shells as all of these objects are characterized by shock excitation rather than photoionization and a high [SII]/Hα line ratio. SNR are defined by a line ratio larger than 0.4, altough Blair & Long (1997) found that this criterion is washed out in some galax- Notes. Column 1, name; Col. 2, Objects in common with Blair & Long (1997): S16 (1) , S14 (2) , S10 ( ies, e.g. NGC7793. We have chosen to assign a classification "shell" if blue stars can be seen within the nebula whose stellar winds are likely shaping the structure, and "SNR" if this is not the case. Fig. 21 shows examples of a SNR and several superbubbles/shells. Although we recorded the presence of luminous blue stars within the nebulae whose stellar winds might contribute to the expansion of the shell, e.g. the WR star d38 in shell d95 (Fig. 17), we note that we are not able at this stage on the basis of the optical data alone to assess with certainty whether the shells are driven by stellar winds, by expanding SNR, or a combination of the two. As an interesting feature, we show three objects that exhibit a double-ring structure whose nature is not exactly clear -perhaps only a by-chance alignment. Detailed studies of extragalactic SNR, including Xray data, were for a long time restricted to the Magellanic Clouds. Recent advances, e.g. the study of Long et al. (2018) in M33, show that more distant galaxies have come into reach. Considering the data quality and sensitivity obtained with MUSE, as well as the advantages of imaging spectroscopy, we are now in the position to collect larger samples within different environments for the purpose of testing numerical simulations, e.g. the SNR models from Li et al. (2015), wind-driven superbubble models from Ntormousi et al. (2011), or the more recent results from the SILCC project (Peters et al. 2017).

Diffuse ionized gas
Results: Diffuse ionized gas, that has been thought for some time to be due to photoionization from field stars and Lyman continuum radiation leaking from H ii regions (Oey & Kennicutt 1997), is ubiquitously found in our dataset. Fig. 19 illustrates the distribution of DIG emission, whose visibility is enhanced by choosing the [S II] emission line map as opposed to Hα, as well as a suitable color table. Apart from bright H ii regions, SNR, and shells that appear in pink and white hues, the greenish structure can mostly be attributed to the ISM. Typical [S II] surface brightness values corresponding to this color coding are from a few up to 50×10 −18 erg cm −2 s −1 arcsec −2 . The presence of DIG is very prominent along the NW spiral arm, and even near the nucleus of the galaxy, with a particularly enhanced level over a roughly 500 pc wide strip extending SW-NW in the lower left quadrant of field (b) to the upper right quadrant of field (a). The interarm regions in the northern parts of field (d12) and in field (j) are comparatively quiet. Still the minimum [S II] surface brightness levels in these areas amount to 2×10 −18 erg cm −2 s −1 arcsec −2 and 2×10 −17 erg cm −2 s −1 arcsec −2 in fields (d) and (j), respectively.
We have investigated in more detail field (a) with the expectation that the old stellar population in the center should exhibit little star formation, hence providing a more undisturbed view on the diffuse emission of the ISM. Within a radius of roughly 200 pc around the nucleus we find indeed no major star forming regions. From Figs. 19 and 23 it is apparent that the distribution of [S II] emission is not smooth and homogenous, but shows a distinct fuzzy and filamentary structure that has as yet not been reported in the literature. In particular, there are several narrow ridges of enhanced surface brightness up to 2 . . . 5 × 10 −17 erg cm −2 s −1 arcsec −2 that extend over lengths of 100. . . 200 pc. In order to obtain alternative surface brightness estimates independently from any unknown systematic errors of our continuum-subtracted narrowband images, and also with the goal to measure radial velocities of the gas, we performed spectrophotometry with P3D over apertures of 5 × 5 spaxels, i.e. 1 arcsec 2 , sampling a total of 20 positions aD01. . .aD20 along filaments, some other high surface brightness areas, and three areas of the lowest apparent surface brightness (Fig. 23).
A correction for the underlying stellar continuum was applied using a ULySS fit to the local unresolved stellar population, that takes into account the LOSV and velocity dispersion. This was subtracted from the flux measured within the aperture defined with the P3D tool. The background spectrum needed to be  since ULySS applies a polynomial fit to accomodate variations of the overall spectral energy distribution, that however does not necessarily match the exact shape of the background continuum at the percent level of accuracy. Since the scaling was iteratively performed under visual control, it was also reassuring to see that the model background spectrum reproduced not only the strong absorption lines, but also numerous weak absorption features. The background correction turned out to be instrumental for the elimination of systematic errors that would otherwise arise strongly at stellar absorption lines. The results from these measurements are listed in Table 12. Pellegrini et al. (2012) have stressed that ionization parameter mapping is effective at highlighting large, extremely faint structures. We have applied the IPM technique exemplarily in field (a). As reported in §4.4, the [S II] surface brightness recorded in the narrowband images agrees generally well with spectrophotometry that was performed in selected areas. We have chosen to measure with P3D specifically (a) prominent filamentary DIG features, and (b) local minima of low surface brightness regions, in order to shed light on the physical nature of those regions. The results of these distinct surface brightness measurements are tabulated in Table 12, presenting the flux in [O iii], Hα, and [S ii] across an aperture of 1 arcsec 2 , and the LOSV of the measured spots. The radial velocities for the filament SF1 north of the nucleus, sampled by aD03. . .aD06 are 145±2 km s −1 , which is within the error identical with the systemic velocity of the galaxy, i.e. quite reasonably close to the nucleus. Interestingly, the branch of the filament SF2 sampled by aD07. . .aD09 exhibits a slightly more positive value of 151±2 km s −1 , which can be understood if SF1 and SF2 are aligned in projection, however are not physically associated with each other, each having their own kinematic history. The [S II]/Hα line ratio map in Fig. 22 shows that filaments SF1 and SF2 stand out from their environment with a value of ∼0.6. Discussion: Following Mathewson & Clarke (1973), we interpret the enhanced line ratio as shock fronts from ancient supernova remnants, rather than due to photoionization from leaking H ii regions. There are more features suggesting the same interpretation, namely filaments SF4 and SF5. SF3 is especially prominent, marking the transition from the extremely low surface brightness area around aD14 towards the bright giant H ii region De100. For comparison, the SNR a64 (marked "S16" in Fig. 14,"SNR" in Fig. 22) appears prominently at the northern edge of the field. The [S II]/Hα map further supports this interpretation, e.g. delineating a region of high excitation gas by an almost perfectly circular arc (SF2), as indicated in Fig. 22 with a circle. The same argument can be made for the arc of SF1, as well as for the pair of SF4/SF5, and SF3. The four examples are somehow similar to the enigmatic large ionized bubble discovered by Pellegrini et al. (2012) in the LMC (their Fig. 10), although the latter is lacking a rim. Ogden & Reynolds (1985) have studied a similar, kinematically distinct filament in the Milky Way, concluding that most likely the object is not photoionized by a nearby O star, but rather shock ionized.
The emission line intensity measurements of aD01. . .aD020, that were chosen specifically along filaments bright in Hα, and then also in regions of local surface brightness minima, were used to determine line ratios and place the resulting values in diagnostic diagrams, such as the one of Sabbadin et al. (1977) that distinguishes shock-excitation gas in SNR from photoionized gas in H ii regions.  Therefore the question arises, whether, perhaps, the emission of DIG is not dominated by shock excitation owing to SNR and stellar winds of massive stars, rather than photoionization. We felt compelled to study this more closely, because already Reynolds (1985aReynolds ( , 1985b who had observed the DIG in the Milky Way using the Wisconsin Fabry-Perot Spectrometer, pointed out a high [S II]/Hα ratio, typical for shock excitation, and concluded that shock ionization was likely to be at work, or perhaps a combination of photo-and shock ionization. The model of Wood et al. (2010), however, demonstrated, that also  image of diffuse ionized gas in field (a). Contours are set to surface brightness levels of 0.5, 1.0, and 1.5×10 −17 erg cm −2 s −1 arcsec −2 . Right: Diagnostic diagram to distinguish photoionization from shock excitation according to Sabbadin et al. (1977). Red plot symbols correspond to shock ionized filaments aD03...aD09, aD15...aD17, etc., and blue symbols to photoionized local surface brightness mimima aD02, aD14, aD18 (see text).
photoionization alone could possibly be the relevant mechanism, owing to the patchy and turbulent nature of the ISM.
As a first test towards a more thorough future quantitative analysis, we have used the region-of-interest feature of the P3D tool to measure the flux emerging from filaments likely to emit shock-excited features in proportion to flux from the remaining areas in field (a) outside of H ii regions to obtain an order of magnitude estimate of the fraction of DIG emission that might be due to Lyman continuum photons leaking from H ii regions, rather than due to shock excitation. As a result, we find a total Hα flux of 1.1×10 −14 erg/cm 2 /s for the filaments (5342 spaxels coadded), and a total of 8.8×10 −15 erg/cm 2 /s for the remaining area (9934 spaxels). In essence, the two contributions are of the same order of magnitude, while it is worth mentioning that one cannot exclude that even the latter may be, to some fraction, also due to shocked gas.
We note in passing that bright blue field stars without a noticeable H ii region, also emission line stars that possibly represent massive luminous stars, show no prominent enhancement of DIG emission in their vicinity that would be comparable to the intensity of the filaments. An example is shown in Fig. 24, where the LBV candidate star i60 is seen in Hα in comparison to H ii regions and SNR within a radius of ∼200 pc. This star emits only in Hα a flux of 6×10 −16 erg/cm 2 /s, which is more than half of the total Hα a flux of the entire optically thick H ii region i62 (to the right in Fig. 24). If the large width of the Hα line (5.3 Å FWHM, extended wings) is indeed an indication of a strong stellar wind, the star is likely extremely hot and luminous, probably outshining any of the O stars that are ionizing the H ii regions in this field. Since there is no indication of an H ii region around i60, one would expect significant leakage of Lyman continuum photons that should ionize the neutral hydrogen known to exist in NGC 300 (Westmeier et al. 2011). The contour levels in Fig. 24 show no evidence that this would be the case (the object i59 to the north is a shock-ionized SNR).

Kinematics
As highlighted in §4.1, the analysis with ULySS and the P3D tool has provided us with LOSV information for stars and emission line sources. The results from our measurements are plotted in Fig. 25. The panel (a) shows the LOSV of individual stars that were extracted with PampelMUSE with a colour code that is explained in the legend, stretching from below 130 km s −1 to above 190 km s −1 in bins of 10 km s −1 . For better orientation, the symbols are plotted over the WFI Hα image shown already in Fig. 14. One can immediately see that on average the stellar LOSV increases steadily with galactocentric distance from a median at 144 km s −1 near the nucleus to around 200 km s −1 at the extreme north-west. We have compared our results with the literature and reproduced to this end the velocity field from Hlavacek-Larrondo et al. (2011) as projected onto our fields (panel (d) in Fig. 25). This data was obtained from observations with a scanning Fabry-Perot instrument that has sampled the Hα line in emission and a post-processing that produced a smoothly varying velocity field without gaps. They found their results to agree with H I data from Puche et al. (1990) and Westmeier et al. (2011) concerning the systemic velocity of 144 ± 2km s −1 , and the distribution of the velocity field. As the sampling of the H I maps is rather coarse (10 arcsec and 30 arcsec grids, respectively), we have preferred a direct comparison to the Hα velocity map. The stellar LOSV follows well the Hα velocity field, with values around the systemic velocity close to the nucleus in field (a), out to 200 km s −1 in fields (d) and (e). This result should be taken with some caution because only the stars in field (i) have been fully analyzed including visual inspection and quality parameter assignment at this stage, such that the output of ULySS fits were taken at face value for the remaining fields, only excluding obvious outliers from failed fits. From our experience with the exercise in field (i) that has shown a robust behaviour with regard to radial velocity determinations, we are however expecting no fundamental change of the picture after full analysis of the remaining fields.
Since we have obtained ULySS fits also for the local background spectra as an output from PampelMUSE for each datacube, there are also radial velocity maps for the unresolved stellar population. For the parameter setting of our PampelMUSE runs, we have obtained unresolved background spectra for tiles with sizes of 20 × 20 spaxels, i.e. areas of 4"×4" on the sky. However, these measurements are hampered by contamination from nebular lines for tiles at or close to the location of bright H ii regions, superbubbles, shells, etc. that mostly rendered a fit impossible. For the scope of this paper, we have not attempted to remove the disturbing emission, e.g. following the discussion of Falcón- Barroso et al. (2004) for the measurement of emissionfree stellar kinematics. We reserve this step to the future analysis of the complete dataset with the pPXF tool (Cappellari & Emsellem 2004), that would also address a stellar population analysis of the background component. The present results are plotted in panel (b) of Fig. 25 with the same colour-code as the Hα velocity map, where black areas indicate that no ULySS fit was obtained due to bright nebular emission lines. Despite the fact that significant areas are not covered with LOSV information, for the tiles that do carry information the velocities are almost indistinguishable from the Hα map. The predominance of G type main sequence star spectra of the ULySS fits for the background in field (a), changing to predominantly B type spectra for the fields associated with the spiral arm is suggestive of a population age ≥ 10 Gyr near the nucleus (Bertelli et al. 2008) and ongoing star formation elsewhere, prompting us to perform a detailed analysis in a forthcoming paper.
Finally, the radial velocities measured for PNe are plotted in panel (c) of Fig. 25 using the same colour code as for the stars. We have also added the PN candidates from Peña et al. (2012) as white circles, for which we have, however, no velocity information. The first striking feature has already been addressed above, namely the unequal distribution of PN detections across our seven fields, despite the fact that we have estimated completeness limits of m 5007 = 27 . . . 28, depending on the seeing attained in each field. The next observation is the apparent scatter of velocities, which for most objects track the general distribution of the Hα velocity map. However, there is a number of remarkable outliers: two high velocity PNe in the nuclear region: (a03) with v rad =197±9 km s −1 , (b06) with v rad =244±3 km s −1 , and at larger galactocentric distances four low velocity PNe: v rad =139±11 km s −1 (d63), 144±7 km s −1 (d79), v rad =145±4 km s −1 (e02), and 143±3 km s −1 (e20), respectively, as well as the extreme case of PN i84 with v rad =108±8 km s −1 . It is going to be interesting to see after the final analysis of our full dataset whether a similar behavior can be found for AGB stars as the progenitors of PNe. We also speculate that we may be able to find evidence for the existence of migration along the leading and trailing edges of spiral arms as suggested by the predictions from numerical simulations by Grand et al. (2016), although it would be premature at this stage to draw any conclusions from our as yet provisional dataset.

Serendipitous discoveries
As mentioned in §3, the process of data analysis has involved significant human interaction and visual control. As an asset, the detailed inspection of our data has enabled the discovery of objects that were not initially targeted, or expected. One such example is the WR star found in field (d) as reported in §4.3 that immediately became apparent upon inspection of the He II image of that region. This star was not specifically classified or ear-marked by ULySS. Other examples are the detection of blue emission line stars representing the rare class of hot massive stars whose discovery conventionally would have involved the twostep procedure of imaging and follow-up spectroscopy, the discovery of symbiotic star candidates, and the detection of carbon stars.  Ta ble 6), that in the HST image appears even 1.4 mag brighter than A, is not detectable in our MUSE image. The less than perfect fit of the K5V spectrum from A is understood as an effect of blending.

Variable sources
Although the layout of our observations was not directly addressing photometric or spectroscopic variability, multi-epoch MUSE surveys inherently offer the capability of detecting variable sources. Incidentally, by comparison of MUSE data with HST images, we have serendipitously made such a discovery: upon the investigation of the apparent foreground star ID 32224 as described in §4.1, the visual inspection of the stellar images shown in Fig. 26 revealed that the star denoted B in the HST image (epoch 2005) is in fact absent from the MUSE image (epoch 2015). Apparently, we have picked up a transient or variable star (perhaps a Mira?), whose brightness has faded over a decade by at least 2 orders of magnitude, reminiscent of a similar case of a red supergiant that was described as a failed supernova by Adams et al. (2017). A possible alternative is that star B has moved closer to the position of star A. In that case B would be the actual foreground star with a rather high proper motion. Unlike an ongoing MUSE GTO program on globular clusters (Giesers et al. 2018), we have as yet not attempted to obtain multi-epoch observations in NGC 300 with the purpose of detecting radial velocity or spectrophotometric variability. However, such an observing strategy of future MUSE surveys in nearby galaxies would open unique opportunities, e.g., to discover O star binaries.

Nuclear star cluster
Like most Scd galaxies, NGC 300 has a nuclear star cluster at its center (Seth et al. 2008), CL1 in Fig. 27. Nuclear star clusters are thought to grow dissipatively from in situ star formation (Milosavljević 2004) and/or by dissipationless accretion of migrating star cluster (Tremaine et al. 1975;Antonini 2013). The at first glance point-like object CL2 in field (a) attracted our attention, in that the spectrum is similar to the spectrum of the nuclear star cluster CL1 that is known to harbour a young stellar population (Walcher et al. 2005;Carson et al. 2015). Inspection of the corresponding HST F606W image revealed that the image of CL2 has a FWHM of 4.0 pixels, that is significantly wider than the PSF with a FWHM ∼2.2 pixels. This confirms that we have found another, less luminous cluster at a projected distance of ∼71 pc to the SW of the nucleus.
The decomposition of the spectrum of the nuclear cluster CL1 with ULySS yields major fractional contributions from spectral types G8V (66.1%), A2IV (10.8%), K2IIIb (9.9%), and CL2 CL3 CL1 v rad =139±2km s −1 , whereas the decomposition of the new cluster CL2 is reporting G8V (64%), B0.5IV (19%), B9III (12%), and AIa (30%) as the major components, indicating an even younger population than the one of the nuclear cluster. The LOSV is determined as v rad =119±3km s −1 . There is another bright object with a distinct blue color in the vicinity of the nucleus, whose HST F606W image is resolved with a FWHM of 3.4 pixels. This new cluster CL3 is located to the NE of CL1 with a projected distance of 35 pc. The spectral decomposition in this case yielded B0.5IV (45%), G8V (32%), B3Ib (22%). The LOSV is measured to v rad =114±10 km s −1 , a value that is significantly lower than the average of 144 km s −1 for the stellar population in this region, which is also the case for CL2. The presence of of a blue supergiant component, as well as a broad Hα emission line in the composite spectrum, points to a cluster age of only a few million years. With a LOSV differential of 20 km s −1 , the cluster would have travelled within a timespan of 1 Myr over a distance of no more than 20 pc on an inclined trajectory through the disk of NGC 300. There is also no evidence for a recent merger or an interacting companion (Bland-Hawthorn et al. 2005). Clusters C1, C2, and C3 are also visible in HST I band images of NGC 300 by Böker et al. (2002). The latest compilation of nuclear star clusters in late type galaxies is published in Georgiev & Böker (2014). The stellar clusters orbiting the nucleus together with the presence of young stars motivates a more detailed study of the assembly of the nuclear star cluster in NGC 300.

Background galaxies
During the extensive use of the P3D visualization tool, we incidentally noticed emission line signatures in the display of colour coded stacked spectra, that did not seem to match any familiar pattern. As it turned out, we discovered several background galaxies in each of our fields, with redshifts of z = 0.13 . . . 1.33, and even two candidates for Lyα emitting galaxies at redshifts Table 13. Background galaxies (redshifts with a colon are uncertain).

ID
X Y z > 4. These objects are listed in Table 13. IDs with an index a,b,... indicate systems with kinematically distinct components, that could either correspond to a rotating disk, or otherwise to a group of galaxies. Although the study of background galaxies was not considered an immediate objective of this work, we note that such galaxies, in particular AGN, may serve as astrometric references for proper motion studies with future instruments delivering very high astrometric resolution and accuracy, such as MICADO for the E-ELT (Davies et al. 2016).

Summary and conclusions
This paper is presenting the first results obtained from MUSE observations of the nearby galaxy NGC 300 in the wide field mode, utilizing the novel technique of crowded field 3D spectroscopy. The dataset used at this stage comprises five complete pointings at the nucleus of the galaxy and adjacent regions to the west with a total exposure time of 1.5 h per field, two of which are compromised by poor seeing. Two more pointings were as yet incomplete with exposure times of 0.5 and 1 h, respectively. The major purpose of this first paper is to demonstrate the feasibility of crowded field 3D spectroscopy in nearby galaxies, and to inform the astronomical community about the legacy value of MUSE datacubes, in particular concerning the multitude of scientific results over a variety of different objects one can expect from such data. Table 14 summarizes the different classes of objects for which individual spectra and images were obtained. The total number of objects, such as stars, PNe, H ii regions, SNR, etc., including background galaxies, amounts to 2187. At this early stage of exploration of the new method of crowded field 3D spectroscopy, we have decided to visually inspect each and every object, assigning quality flags, and reassuring ourselves to not be trapped by immature or not fully reliable automatic algorithms. Concerning the classification of stars, this task has so far fully been accomplished for field (i). Numbers flagged with a colon for the remaining fields in Table 14 indicate that the classification is preliminary. However, significant parts of the analysis have already been automized: the extraction of stellar spectra with PampelMUSE, and spectral type classification as well as radial velocity measurements with ULySS and GLIB fits. The significant level of human interaction at this stage is timeconsuming and cumbersome. We are actively working on the improvement of the currently available tools and expect that future developments will enable fully automatic procedures, which is a prerequisite for the analysis of surveys larger than the present study.
As the major outcome of this pilot study, that can be considered a proof-of-principle for MUSE crowded-field 3D spectroscopy in nearby galaxies, we summarize the following results: 1. It is possible to decompose nearby galaxies into individual giant stars and an unresolved component of faint background stars. From the analysis of field (i), which is the best one in our dataset in terms of seeing, we present spectra for a total of 517 stars whose S/N was sufficient for spectral type classification and the determination of radial velocities ( S/N ≥ 3). As inferred from simulations, the limiting S/N for a reliable spectral type classification was found to be in the range of 3. . . 5, depending on spectral type. From the remaining, as yet uninspected fields, 1329 stars have yielded spectra, part of which were used to measure preliminary radial velocities, subject to quality control and confirmation in the near future. Of the immediate findings from this study, we point out the discovery and classification of ∼100 blue supergiants, tens of yellow giants/supergiants, and more than 300 K and M giants/supergiants in field (i). The radial velocities measured for individual stars follow the rotation curve of the galaxy in accord with published velocity fields for Hα and H I. Likewise, the velocity field determined from ULySS fits to the unresolved background stellar population is consistent with this data.
2. From the single field (i) with an exposure time of 1.5 h, we identified a total of 23 carbon stars. This number compares with a total of 115 carbon stars found in the disk of M31 (Hamren et al. 2015) and the surrounding halo and satellite dwarf galaxies (Hamren et al. 2016), that were obtained from a total of 24762 stellar spectra secured over 10 years worth of Keck DEIMOS observations as part of the SPLASH survey.
3. In addition to the stars selected and extracted on the basis of the ANGST catalogue, we found 82 blue emission lines stars, most of which show broad Hα emission lines, that are characteristic for evolved massive stars with strong stellar winds, e.g. LBV stars. Such stars are very rare and difficult to find (Massey et al. 2015). For comparison, Massey et al. (2016) reported the discovery of 4/5 LBVs and 0/1 WR star in M31/M33, respectively. Point sources that show a continuum spectrum and unresolved Balmer lines in emission may represent circumstellar shells around stars that are small enough to remain spatially unresolved for ground-based observations.
4. Another unexpected discovery amongst the detected emission line stars is the finding that some of the emission line point sources are apparently associated with AGB stars, some of which are interpreted as the discovery of candidate symbiotic stars. We identified a total of 4 symbiotic star candidates that would match such a spectral pattern. It is interesting to note that according to the NASA HEASARC archive merely 188 such stars are known in the Milky Way (Belczyński et al. 2000). According to Mikołajewska (2004), 6 such objects are known in the SMC, and 8 objects in the LMC.
5. Both the analysis of stellar spectra through automized fits to the MIUSCAT and GLIB libraries, and the interactive measurement of emission line objects using P3D, have yielded LOSV estimates for hundreds of individual objects. While the stars match very well the velocity field of H II, which should be characteristic for the kinematics of the disk, amongst the sample of PNe, that generally follow the same trend, there are several distinct outliers with extreme radial velocities that may be attributed to halo objects. The velocity field, that was obtained from local background estimate spectra as output from the PampelMUSE code, although as yet uncorrected for gaseous emission systematics, is also indistinguishable from the H II velocity field.
6. From the complete analysis in field (i), so far only one foreground star was found on the basis of luminosity class and radial velocity, in accord with the expectations.
7. Amongst emission line point sources we have registered 36 bona fide PNe (out of which 13 objects were known from previous studies), and 9 more uncertain candidates at faint magnitudes. Extragalactic PNe can easily be confused with unresolved compact H ii regions that show emission line spectra similar to the ones of PNe. From previous PN surveys in NGC 300, we found several cases where objects were misclassified as PN and better match the criteria of cHII. In total, we registered 61 cHII.
8. We discovered a total of 53 new extended H ii regions and confirmed 8 that were already catalogued previously. On the basis of the IPM technique (Pellegrini et al. 2012), we assigned categories of optically thick or thin H ii regions, and shells.
9. Supernova remnants were identified on the basis of the [S II]/Hα line ratio. We registered a total of 38 SNR candidates, with sizes ranging between a few and more than 100 pc. Some of the SNR candidates exhibit very low surface brightness, and therefore a possible confusion with the DIG. 10. We measured DIG emission lines in the nuclear region of NGC 300 down to surface brightness levels as low as a few 10 −18 erg cm −2 s −1 arcsec −2 and identified filamentary structures that show diagnostic emission line ratios that are characteristic for shock excitation, whereas only a few areas with very low surface brightness were found whose spectra are characteristic for photoionization.  11. Two young clusters within a projected distance of less than 75 pc from the nucleus of NGC 300 were discovered.
12. Of serendipitous discoveries made in the process of data analysis, we mention that 28 background galaxies with redshifts z=0.127. . .1.33 were identified, that otherwise would likely have remained undetected from imaging surveys.
It is obvious that these early results have merely scratched the surface of a potentially very powerful new tool to probe resolved stellar populations in nearby galaxies, a research topic that we are planning to expand in the future.
The data products from this work will be made publicly available through the Strasbourg astronomical Data Center (CDS) and the MUSE website 5 .

Outlook
As a technical shortcoming at the present stage, we have found that the PampelMUSE technique to account for contributions from the unresolved background as a constant average over a finite area (the "tiles"), which works well for globular clusters, is a major limitation in regions with strong nebular emission lines. Of the stars with acceptable S/N, for which however no trustworthy classification could be obtained, the majority was found to be suffering from strong nebular contamination. We concluded that mainly massive stars or star clusters within bright H ii regions are affected by this approximation. In previous experiments with the PMAS instrument published by Roth et al. (2004), we have separately modeled the 2-dimensional background surface brightness distribution in the continuum and in emission. For an instrument like MUSE that has a factor of ∼350 times more spaxels than PMAS, such an approach is computationally very expensive. The further development of PampelMUSE will specifically address this issue.
Another limitation was imposed by the incomplete coverage of the HRD by the MIUSCAT and GLIB libraries, especially for hot stars, emission line stars, and also for carbon stars. We are planning to resolve this problem by complementing the libraries with simulations and new observational data. The MaStar library currently being planned as part of the SDSS-IV collaboration (Yan et al. 2017) may become a promising solution to this end.
Future work shall address the completed dataset of stellar spectra of all fields with secured quality checks, including two new pointings at larger galactocentric distance that are not shown in this paper. We are also planning new observations of the fields that have suffered from poor seeing with support of the adaptive optics facility, which should substantially increase the number of useful stellar spectra. Moreover, we are planning to improve 5 http://muse-vlt.eu/science/ our methods towards automatic data analysis to more efficiently process the expected number of thousands of spectra.
From the entire data set, we intend to study the stellar population from the nuclear region to the north-western spiral arm, including the leading and trailing interarm regions covered by our pointings, with the goal to infer the star formation history on the basis of individual stars and their relation to H ii regions, SNR, and PNe, complemented by spectral synthesis modelling of the background of unresolved stars. Furthermore, from a detailed kinematic study we expect insight into phenomena like migration, mergers, and runaway stars. Another challenging goal will be the determination of chemical abundances for individual stars, with interesting prospects concerning young and old stellar populations, e.g. blue and red supergiants vs. RGB and AGB stars.
Again on the basis of the complete dataset, we shall study the origin of the PNLF and influences on its shape owing to the underlying stellar population, in particular M and C stars as their progenitors. This goal has become more accute with the revised picture of post-AGB evolution from the latest models of Miller Bertolami (2016), that are three to ten times faster and ∼0.1-0.3 dex brighter than the previous models by Vassiliadis & Wood (1994) and Bloecker (1995), that were being used in many studies.
With the hypothesis that DIG emission may be dominated by shock excitation, we are planning a follow-up paper to systematically investigate all of our pointings in more detail in order to obtain a more global view of the characteristics and origin of DIG in different parts of the galaxy.
For the modest investment in observing time (9 h) the yield as summarized in Table 14 is impressive. It has been a fundamental goal of the development of the MUSE instrument to break the conventional sequence of photometry -follow-up spectroscopy, and replace it with a single-step observation, that provides spectra and images from a single homogenous dataset. Beyond this already fundamental achievement, it is worthwhile to stress that not just quantitatively, but also qualitatively it has been only the format and performance of MUSE that has enabled the efficient use of crowded field spectroscopy, both in terms of multiplex advantage, as well as in utilizing the PSF fitting technique to deblend overlapping stellar images. Combined with high instrumental throughput, the light collecting power of an 8m telescope, and the superb image quality of VLT-UT4, we conclude that MUSE is an ideal instrument to study (partially) resolved stellar populations, gas, and dust in nearby galaxies. We expect that this capability will be even more enhanced with the new adaptive optics facility for VLT-UT4 that has become available recently. Our findings also strongly support new survey telescope concepts suitable for massively-multiplexed spectroscopy as presented by Pasquini et al. (2016).