© ESO 2018

1. Introduction

The quest to understand 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 fiber, or a slit (see Sanchez et al. 2012). The CALIFA survey is the first effort toward 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) such as 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 IFU, it is not entirely clear whether the technique of stellar population synthesis for the resulting spectra is capable of uniquely determining stellar populations and their star formation history. The study of resolved stellar populations in nearby galaxies, where stars can be measured accurately on an individual basis, would be able to lift any degeneracies and allow the calibration of the integral light measurements of more distant galaxies. Such work has been attempted with Hubble Space Telescope (HST) photometry, for example, 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, for example, age-metallicity degeneracy, sensitivity to extinction, or the inability to distinguish objects such as Wolf–Rayet (WR) stars from normal O stars, carbon stars from M stars would rather call for spatially resolved spectroscopy. Except for conventional slit spectroscopy of individual bright giants or supergiants in local volume galaxies, for example, Bresolin et al. (2001) and Smartt et al. (2001) in the optical, or Urbaneja et al. (2005) and Gazak et al. (2015) in the near-infrared (NIR), no comprehensive spectroscopic dataset exists to date for stars in heavily crowded fields of nearby galaxies beyond the Magellanic Clouds. Pioneering attempts to use the technique of IFS for crowded field spectroscopy, i.e., fitting the point spread function (PSF) to point sources through the layers of a data cube, analogous to CCD photometry in a single image, for example, 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; see, for example, 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, for example, 8 × 8 arcsec² at 0.5 arcsec sampling in the case of PMAS (Roth et al. 2005), the light collecting power of 4 m 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² 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.5 m 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 ≈60 s snapshot exposures, resulting in a spectacular dataset of spectra for more than 12 000 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, and 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) again studied the radial velocities and velocity dispersion of individual stars, however with a sample size of more than 12 000 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⁻¹.

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: asymptotic giant branch (AGB) and red giant branch (RGB) stars, the most massive O stars, blue and red supergiants, luminous blue variable (LBV) and WR stars. As opposed to the Milky Way, where severe selection effects are at work, for example, dust extinction, small numbers of known objects, lack of accurate parallaxes (see Crowther 2007), we 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- and nar-rowband filter imaging, and then follow-up spectroscopy (see Massey et al. 2015), IFS allows us 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 us 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; Arnaboldi 2015a,b). Modelling internal dust extinction, the discovery of supernova remnants (SNR) on the basis of Hα/[S II] and [O I]/Hβ line ratios (Fesen et al. 1985); the distribution and physical properties of H II regions, and the discovery of compact H II regions, in particular; and diffuse and filamentary interstellar medium (ISM) are further pieces of information that can be retrieved from the data cube.

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 to demonstrate the feasibility of the proposed approach. Basic parameters of this galaxy are listed in Table 1. We report the first results of our GTO observations in ESO periods P93, P94, P95, also demonstrating new legacy value opportunities provided by MUSE. We stress that any data cube hosts copious amounts of information in addition to the objectives of the science case proper that have led 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 dataset with a focus on different science cases (see Sect. 6).

The paper is organized as follows: Sect. 2 presents details of the observations and data reduction; Sect. 3 describes the data analysis; and Sect. 4 summarizes the results with catalogs of stars, PNe, 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.

2. 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.2 m 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⁻¹. 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 basic 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 3D 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 and standard star exposures. The average atmospheric extinction curve of Patat et al. (2011; included in the pipeline) was used when creating the sensitivity functions. Sensitivity functions and 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 90 000 spectra of the sky cube, which was subtracted from the extracted dataset of the respective OB. The two to six extracted datasets 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 for this purpose. The resulting cubes were created using the standard sampling of 0″̣2 × 0″̣2 × 1.25 Å.

Fig. 1. Footprint of MUSE pointings (a), (b), …, (j) in NGC 300, covering the nucleus (a), different fractions of the spiral arm that extends to the NW of the center (b), (c), (d), (e), (i), and the inter-arm region (j). Orientation: N up, E left. (Image credit: ESO).

3. Analysis

The data analysis was carried out 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 performed 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 data cubes using DS9, respectively. P3D is an open source IFS software package developed at the AIP and publicly available under GPLv3 from SourceForge¹ (Sandin et al. 2010). Previous applications of P3D for the VIMOS 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.

3.1. 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 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 chose to define three broadband transmission curves V, R, I, which 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. Figure 2 shows color composite images of our fields, overlaid onto a wide-field color composite image obtained with the ESO 2.2 m telescope.

Furthermore, we defined narrowband filters for important emission lines such as Hα, Hβ, [O III], [N II], [S II], [O I] with central wavelengths adjusted to the redshift of NGC 300, and wide enough to cover the instrumental profile of MUSE, respectively, in order to not lose any flux. The filter width, however, was chosen to be small enough to suppress continuum background light as much as possible. With a typical narrowband filter width of five 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.

Fig. 2. VRI images reconstructed from MUSE data cubes 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 color, 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.2 m WFI (credit: ESO).

3.2. 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²). 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 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 manually 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 PNe 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, kurtosis β, ellipticity e, and 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 brightnesses in the reference catalog, a mock MUSE image is created for a predefined 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 by estimating the S/N of each source based on its magnitude in the reference catalog, PSF, and 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 catalog 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 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. Afterward, 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 behavior, 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 repeated with 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 catalog number per field as listed in Col. (1) of Table 6.

3.3. 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 classifi-cations, and to measure radial velocities, even with moderate to low signal-to-noise ratios (S/N). To this end, we 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). The MIUSCAT catalog 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), which 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³, which 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 ten) 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 us to make a good guess of the spectral type of the observed star. However, any less than plausible spectra can be ruled out immediately; for example, 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)₀ = 26.37 (Gieren et al. 2004, 2005) for comparison with the apparent magnitudes from the ANGST catalog, 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 VRI map from the MUSE data cube as a post-stamp-like image of a size of 4 × 4 arcsec², accompanied with an HST ACS image of the same region, color 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), which 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), the most probable radial velocity, and a global quality flag were determined for the final result. The adopted outcome of the fits and visual inspection were recorded in individual log files and summarized in a catalog 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 Sect. 4.1. An excerpt of the catalog is presented in Table 6.

Fig. 3. MUSE vs. HST images of 4 × 4 arcsec2 size from field (i), as examples of the star in the center being affected by different severity of crowding. a/b: ID1154, A1 Ib star, unblended. c/d: ID1179, B3 III star, negligible blending. e/f: ID995, central star cluster in bright H II region, spectrum heavily affected by blending and nebular emission, no reliable classification possible.

3.4. Extracting nonstellar sources

As expected from the outset, visual inspection has shown indeed that our MUSE data cubes 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 such as H II regions, SNR, supershells, diffuse ionized interstellar gas (DIG), and PNe. Also faint background galaxies reveal their presence through redshifted emission lines that stand out from the foreground galaxy continuum. The search and classification of such objects was carried out by visually inspecting the emission line maps described in Sect. 3.1 with DS9. H II regions were readily identified on the basis of spatial extension and brightness. The PNe were detected by blinking [O III] versus Hα and He II to find high excitation objects. Supernova remnants, superbubbles, and supershells were distinguished from H II regions by the technique of ionization parameter mapping (IPM) 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 data cube with Pam-pelMUSE. 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 Col. (1) of Tables 7–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 IPM 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 P3 D 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 to accommodate 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 such as 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 Sect. 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 catalog, without consulting the PampelMUSE catalog to avoid any subjective bias in the detection process. By blinking against the [O III] imagesf we identified high excitation objects such as PN candidates, also finding 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 catalog 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 us 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 line-of-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 as follows:

$${\sigma }_l={[ {\sigma }_{\text{cont}}^2 N (1+EW/N\Delta )+{(0.05\times F_l)}^2 ]}^{1/2}$$(1)

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 double check 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.

4. Results and discussion

4.1. 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 northwestern spiral arm at a galactocentric distance of ≈1.5 kpc. We adopted both of the two procedures introduced in Sect. 3.2, i.e., the input of stellar centroid priors from high spatial resolution HST images, and alternatively from a search in the data cube using DAOPHOT FIND. PampelMUSE was run with both catalogs, resulting in a total of 3540 and 552 extracted spectra, respectively. It should be noted that the HST coverage is only two-thirds 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 catalogs, 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 data cube, 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 Sect. 3.2, we created the final catalogs for the two sets of spectra. An excerpt of the catalog for field (i) with HST input is presented in Table 6. The complete catalog for field (i) is 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 top panel star numbered ID119 has a S/N = 19 and is classified A1Ib with v_rad = 149 ± 18 km s⁻¹, the spectrum in the middle from ID260 has S/N = 14 and is classified K3Iab with v_rad = 171 ± 17 km s⁻¹, the spectrum in the bottom panel from ID1379 has S/N = 6 and is classified M6III with v_rad = 174 ± 7 km s⁻¹. 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 lines to the best MIUSCAT fits, which are for the most part almost indistinguishable from the measured data.

Figure 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, which is clearly identifiable through its prominent TiO absorption bands, whereas 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, which are 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; for example, less then 20% contamination is from neighboring stars (Massey et al. 2016). The HST input catalog stars that have delivered useful spectra cover a magnitude range in the F606W filter of 20 … 25 mag.

Figure 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). Because of 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 broad with FWHM(Hα) = 6.5 Å, i.e., not of nebular origin. The Hα flux was measured to 4.0 × 10⁻¹⁸ erg cm⁻² s⁻¹, and the radial velocity as 144 km s⁻¹. The collapsed broadband image from the MUSE cube shows merely a vague blue hue and illustrates the fact that for ground-based observations only IFS opens the opportunity 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 plausi-blity 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 cases are found mostly at S/N levels below 4. Secondly, marginally plausible fits show a similar behavior and only a handful of cases are 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 bottom panel shows the distribution broken down by effective temperature; the blue symbols correspond 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, for example, 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, for example, measured as FWHM = 0.60″ at 460 nm versus 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 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 already works with a catalog 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 catalog 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⁴. 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 colors (F606W–F814W) ≥ 0.4; blue symbols indicate hot stars with (F606W–F814W) < 0.4.

Both datasets show a very similar distribution that has a tight correlation at bright magnitudes and a rapid degradation toward large errors at faint magnitudes. Residuals against the 1:1 relation in magnitude bins of 0.2 mag are shown in the bottom right part of the graph (with an offset of 20.5 mag). The error bars indicate the standard deviation in each bin and have 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 toward bright magnitudes as the object photon shot noise dominated regime with a slope of –2, as expected. On the other hand, the branch toward 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 ten 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. Figure 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.

We 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 strongly depends on the spectral type. For example, K0 and K4 subgiants still show 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, which 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, for example, 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 stars show a different behavior, 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^b4 is the number of blue stars with S/N > 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 there 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 clas sified by the procedure explained above, comprising roughly two-thirds of the spectra extracted from the HST input catalog, and complemented for the remaining one-third 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, which was selected to not be dominated heavily by ongoing star formation, i.e., a minimal number of bright H II regions, to maxi-mize the detection rate of faint PNe. We find an appreciable number of O star candidates, including a significant percentage of hot emission line stars (11 out of 38, i.e., 29%). This is the case despite the problem of our current approach using the MIUSCAT and GLIB libraries, which do not allow for an accurate classification of hot stars because of the lack of O template stars in MIUS-CAT and the lack of non-LTE models for GLIB; this is 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 have obtained a complete sample of spectral types B...M brighter than or equal to luminosity class III.

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 the calcium triplet, were only marginally discernible from the noise. The measurements merely failed in cases where excessive residuals from very strong nebular emission line background, which 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⁻¹ and have a mean of 169 km s⁻¹ and a dispersion of 23 km s⁻¹ 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⁻¹ is indeed sampling the rotation curve of NGC 300. This finding is further discussed in Sect. 4.7 and Fig. 25.

Figure 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⁻¹, and the 78th percentile at 40 km s⁻¹. 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, for example, 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 also discovered a number of emission line stars from the visual inspection after automatic processing for MIUSCAT and GLIB fits. Their spectra typically show 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 Sect. 4.3.

For resolved stellar population studies it is necessary to correctforcontaminationofthesampleunderstudybyforeground stars that 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⁻¹ and 69 ± 3 km s⁻¹, 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 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 to be expected in the direction toward NGC 300 down to V = 20–0.5 stars per 1 × 1 arcmin². 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⁻¹. 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 catalog 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 dataset, to allow for quantitative spectroscopy along the lines of, for example, Bresolin et al. (2002) and 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 star in common with the Bresolin et al. (2002) sample that is suitable to make a direct comparison: their object B11, classified as A5 su-pergiant (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, and 1.2″ seeing. The star is listed as ID25 from a preliminary analysis of field (b). The MIUSCAT-based classification of this star is A1Iab … A5II, in agreement with Bresolin et al. (2002), despite a lower S/N owing to poorer observing conditions and a shorter exposure time. Figure 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 favorable observing conditions for the MUSE spectrum (exposure time, seeing), and despite the fact that in this wavelength region the efficiency curve of MUSE experiences a steep drop toward the blue.

The importance of the statistics of red and blue supergiants was stressed by Massey et al. (2017). From a review of massive stars by Massey (2003), and a more recent paper by Massey et al. (2016) describing a survey in M31 and M33, it follows that the amount of work necessary to safely detect and classify massive stars is significant. 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 energy 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. The MUSE instrument allows, for the first time over a reasonable FoV, to combine the different survey strategies into one: imaging and spectroscopy combined in a single data cube 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 data cube at the wavelengths of interest, for example, 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 analyze 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 ten years worth of data from the SPLASH survey (Guhathakurta et al. 2006; Tollerud et al. 2012). The ratio (C/M) of carbon-rich to oxygen-rich 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. These 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 14 143 stellar spectra taken with DEIMOS at the Keck-II 10 m telescope. The spectra were obtained from 151 individual DEIMOS masks for 60 separate fields and have 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 to obtain 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 such as those cited above must be corrected for foreground contaminant stars. This is usually accomplished spectroscopically, for example, 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 Sect. 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 Sect. 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 because of the mediocre seeing of the observation in field (b). The S/N estimate for this spectrum is 17, decreasing somewhat toward 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 even though 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 relatively 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 as 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 Sect. 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 inspect visually all spectra one by one 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.

Fig. 4. Example spectra from field (i). Images in the two columns to the left are analogous to Fig. 3. The spectra cover the full extended MUSE wavelength range of 4600 … 9300 Å. Flux is plotted in units of 10−20 erg cm−2 s−1 Å−1. Black curves: observed spectra; orange: ULySS fits.

Fig. 5. Deblending of two overlapping stellar images in field (i). The MUSE and HST maps to the left are analogous to Fig. 3, for units of spectra see Fig. 4.

Fig. 6. Deblending of the extremely faint emission line star (A) from an unresolved group of K stars (C). For an explanation, see Fig. 5 and text below. The insert spectrum covers the wavelength range 6200 … 7000 Å, flux is plotted in units of erg cm−2 s−1 Å−1 (continuum background not subtracted).

Fig. 7. Distribution of S/N as a function of F606W magnitude in field (i), selectively by quality flag (top), and Teff (bottom).

Fig. 8. MUSE F606W photometry versus HST ACS magnitudes, broken down into cool stars (red dots) and hot stars (blue dots). Mean residuals against the 1:1 relation in 0.2 dex mag bins are shown as black dots. The error bars indicate the standard deviation while the gray horizontal line indicates zero.

Fig. 9. Recovered Teff as a function of S/N for a total of 630 simulated spectra, as determined by ULySS. For details, see text.

Fig. 10. Radial velocity uncertaintites ∆vrad plotted against S/N. The sample is segregated into hot and cool stars as in Fig. 7. The insert shows a histogram of the ∆vrad distribution, the 50th percentile at 20 km s−1.

Fig. 11. MUSE spectrum of A5 supergiant B11 from the sample of Bresolin et al. (2002), seen in field (b), ID25. This star has a V magnitude of 19.95. The MIUSCAT fit is plotted in orange with an offset for clarity. Wavelengths in units of Å, flux in units of 10−20 erg cm−2 s−1 Å−1. The reference wavelengths of a number of absorption lines (hydrogen Balmer and Paschen series, helium lines) are indicated.

Fig. 12. MUSE spectrum of AIa star of magnitude mF606W = 20.78 in field (i), ID151. The MIUSCAT fit is plotted in orange with an offset for clarity. Wavelengths in units of Å, flux in units of 10−20 erg cm−2 s−1 Å−1.

Fig. 13. MUSE spectrum of M1Ia star of magnitude mF606W = 21.56 in field (i), ID301. The MIUSCAT fit is plotted in orange with an offset for clarity. Wavelengths in units of Å, flux in units of 10−20 erg cm−2 s−1 Å−1. The reference wavelengths of absorption lines of calcium, magnesium, sodium, Hα, and Hβ are indicated.

4.2. Planetary nebulae and compact H II regions

Results. The analysis of our data cubes 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 of 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.2 m 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, for example, from the giant H II region, denoted 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), which cannot be found 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 to infer its diagnostic value for stellar population studies. Jacoby (1989) introduced the m₅₀₀₇ magnitude, based on the line flux of [O III] at 5007 Å that is conventionally used to represent the PNLF, i.e.,

$$m_{5007}=-2.5 \log F_{5007}-\mathrm{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 used to measure distances to more than 60 galaxies; see review by Mendez (2017). 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 the underlying stellar population (Ciardullo 2010) and eventually help us 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 4 m 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.2 m 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, deblend the spectra 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 data cubes (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₅₀₀₇ = 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₅₀₀₇ = 28.0 for field (i) to 1.2″ and m₅₀₀₇ = 26.0 for field (b). The faintest PNe that we detected in the experiment have m₅₀₀₇ = 30.0.

After having convinced ourselves that the data quality of our data cubes would be sufficient to search for faint PNe, we performed the search and classification procedure as described above. Typically, extragalactic PN candidates appear as point sources that are brighter in [O III] than in Hα. This fact used to be the canonical criterion for the selection of PN candidates for PNLF distance determinations of galaxies, practically applied by blinking narrowband images in [O III] versus Hα. However, a problem with reliably classifying PNe consists in the possible confusion with compact H II region (cHII), which may appear as [O III] bright point sources at the distance of NGC 300. In adopting the argument of Stasińska et al. (2013), we therefore excluded all PN candidates whose spectra show a [S II] line ratio F(6716)/F(6731) > 1 as indicative of the low density limit, and in such cases classified them cHII. Also, Hβ brighter than [O III] λ 4959 is a good indicator for cHII, which is not affected by extinction in contrast to the [O III]/Hα line ratio. The results from the analysis is presented in Table 7, in which we list a total of 45 PNe 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₅₀₀₇ ≈ 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₅₀₀₇ = 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 were 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. comm.).

Table 7 also lists the radial velocity estimates that were derived from the emission line fits with P3D. The overall distribution matches well that of stars, thus giving confidence to indeed observing the rotation of the galax. However, there is one one exception: PN i84 exhibits a velocity of 108 ± 13 km s⁻¹, which is far off the median in field (i) and is also blueshifted from the systemic velocity of 144 km s⁻¹. Thus, this either indicates 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₅₀₀₇ = 28.0 and m₅₀₀₇ = 26.0 for a seeing FWHM of 0.6″ and 1.2″, respectively. The PN candidates that were identified in field (e) reach magnitudes even fainter than m₅₀₀₇ ≈ 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 compared our detections with results from the literature, namely Soffner et al. (1996), who used narrowband imaging at the 3.5 m ESO-NTT, and Peña et al. (2012), who used data from narrowband imaging and spectroscopy using FORS at the VLT. Firstly, the MUSE observations have recovered all of the previous detections. We 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₅₀₀₇ = 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₅₀₀₇ = 27.

We also compared our m₅₀₀₇ 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 Sect. 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 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 colored 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 favor 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 it is supported by the sparse population in field (j) and the interarm extension to the southeast (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 of observable PNe at any time is given by the evolutionary flux that is tied to the underlying stellar population. The PNe that we observe stem from an intermediate age stellar population of ∼1 Gyr, so we would expect the PN surface density to correlate with the stellar mass surface density of the old stellar population, however not with young O stars and H II regions. As we found evidence for contamination of their sample from cHII in 30% of the cases of PN candidates from Peña et al. (2012), it is not excluded that the provisional picture of PNe aligned with spiral arms is simply an artifact.

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 data cubes. 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) estimated PN surface densities of 1.6 … 3.1 PN kpc⁻² 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⁻² for fields (a) and (e), and 16.6 PN kpc⁻² 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.

Fig. 14. Hα images reconstructed from MUSE data cubes for pointings (a), (b), …, (j). The background image is from the ESO/MPI 2.2 m WFI (credit: F. Bresolin). H II regions in common with the catalog of Deharveng et al. (1988) are indicated with red circles. Yellow circles denote the SNR recorded by Blair & Long (1997).

Fig. 15. Left: examples for artificial PN point sources inserted in data cubes of fields (a) and (b), with seeing of 0.7″ and 1.2″, respectively. Right: fraction of artificial PN detections as a function of m5007, obtained with DAOPHOT FIND for fields (a), (b), (c), (i), (j), and seeing of 0.7″, 1.2″, 1.0″, 0.6″, 0.85″, respectively.

Fig. 16. Comparison of PN and cHII magnitudes obtained with MUSE versus results from Soffner et al. (1996) in red and Peña et al. (2012) with blue plot symbols.

4.3. 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 catalog of PampelMUSE-extracted spectra. This is 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 because emission line stars are not part of our empiricial library nor would they be supported by GLIB. In such cases, we 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, for example, Klein & Castor (1978), Leitherer (1988), Lamers & Leitherer (1993), Puls et al. (1996), and 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 spread function 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⁻¹, 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 in which 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 northwestern (NW) spiral arm in NGC 300, another at the leading edge of this spiral arm, and possibly 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 neb-ulosities from intrinsic stellar emission on the basis of the MUSE data. Again, we visually inspected all of the spectra and decided which objects would qualify as emission line star candidates. The results of this analysis are listed in Table 9.

Direct comparison with VRI color 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⁻¹¹ erg cm⁻² s⁻¹. Pushing this object to the distance of NGC 300, this would translate to a flux of roughly 10⁻¹⁸ erg cm⁻² s⁻¹, i.e., just detectable with MUSE. We therefore flagged red giant emission line stars as symbiotic star candidates.

A special case is the WR star 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).

Fig. 20. Example spectra of PNe, cHII, SNR, and an emission line star in units of erg cm−2 s−1 per 1.25 Å spectral bin in the Hβ-[O III] and Hα regions, plotted against wavelength in Å. The spectra are offset by some arbitrary bias to better show the noise level. Panels a–d: PNe over the full range from the brightest to the faintest objects. The very faint object a03 with m5007 = 28.88 is an uncertain PN candidate, given that the [SII] line ratio, although critically dependent on diffuse nebular background subtraction, could also be indicative of a cHII, or even a SNR. Panels e, f: a bright and a faint cHII, and panels g, h: a bright and a faint SNR. Panel i: an emission line star is shown, presenting an Hα emission line width of 6.6 Å.

Fig. 17. WR star in field (d). Top left: VRI map, candidate star indicated as object d38. Top right: Hα map, showing surrounding shells. Bottom: spectrum showing the He II feature (arrow), in arbitrary units.

4.4. H II regions

Results. In addition to emission line point sources, we catalogued extended objects from the visual examination of emission line maps, appearing bright in Hα and/or in [O III] 5007 Å, with a variety of morphologies such as contiguous high surface brightness regions with circular or elliptical geometry, bright regions with irregular geometries, contiguous patches of low surface brightness, ring-like shells, diffuse and filamentary low surface brightness emission, and large scale arcs or closed shells.

For example, the H II region labeled De100 in Fig. 14 bears some similarity with the 30 Dor region in the LMC. There are a number of prominent shells of diameters up to 100 pc visible, for example, prominently in field (d). Aside from 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 made an attempt to distinguish between classical pho-toionized, 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). Figure 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⁻¹⁶ erg cm⁻² s⁻¹ arcsec⁻² (left). The IPM grayscale image, which have 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⁻¹⁶ erg cm⁻² s⁻¹ arcsec⁻² 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 catalog. The [S II]/[O III] line ratio map is also sensitive to reveal other interesting features. Aside from the brighter PNe appearing as black blotches that were already discovered otherwise as described above, the map brings out high excitation emission, for example, a peculiar structure seen as extending radially from the rim of i97 (arrow in Fig. 18, top row, middle). This structure not visible in any other emission line rather than [O III] and is roughly, but not exactly aligned with a blue star. The nature of this strange feature is unclear.

The catalog of H II regions is given in Table 10. Emission line fluxes, radial velocities, positions, sizes, and morphology were determined as described in Sect. 3.4. Also the compact H II region a32 that is embedded in a high surface brightness area of the H II region a65, which is impossible to discern on the basis of an Hα image, was discovered this way. Altogether, we discovered 112 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⁶ yr that require only one episode of star formation, and supergiant shells that have sizes of ∼1 kpc, dynamical ages of ∼10⁷ yr, 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), and furthermore as SNR candidates, shells, or diffuse/filamentary ISM. The SNR and DIG are discussed further below.

Discussion. The most comprehensive catalog of H II regions in NGC 300 is still that 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 yr ago. Since we were able to sample low emission line intensities down to levels below 10⁻¹⁷ erg cm⁻² s⁻¹ arcsec⁻², 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 insofar as their objects are entirely contained in our fields. As a convenient feature of the MUSE data cubes, 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 too large in size to safely resolve, for example, OB associations.

Therefore, the use of HST images and catalog 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 obtain a unified picture concerning leaking Lyman continuum radiation and the ionization of the DIG (see below). The well-known 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), which 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, for example, 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 to better constrain the distribution of the metal within this galaxy using different metallicity indicators and combine our data with those of the literature.

Fig. 18. Ionization parameter mapping of two H II regions, from left to right: Hα, [S II]/[O III], VRI. The grayscale for the Hα maps shows high intensity in black, whereas the IPM maps feature high ratios in white. The region i97 (top row, 4.8″ × 5.8″) is optically thick, while i81 (bottom, 10.4″ × 11.8″) is apparently leaking ionizing radiation, a “blister” in the nomenclature of Pellegrini et al. (2012). A strange, [O III]-bright feature near the rim of i97 is indicated by an arrow (see text). The stars that are likely responsible for the ionization of the nebulae are visible as blue objects in the VRI images.

4.5. 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) and low surface brightness objects that barely stand out from the noise. In our dataset, we reached limiting surface brightness levels of 2 × 10⁻¹⁷ erg cm⁻² s⁻¹ arcsec⁻² 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 such nebulae from H II regions.

Discussion. The 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.5 m 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⁻¹⁶ erg cm⁻² s⁻¹ arcsec⁻². For comparison, the Hα surface brightness limit of our data is below 1 × 10⁻¹⁷ erg cm⁻² s⁻¹ arcsec⁻², 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 data cubes is 38 and have 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 and have 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 that reported for the galactic SNR Pup A by Dopita et al. (1977), or several such objects in M 81, M 101, and NGC 6946 Matonick & Fesen (1997).

The dichotomy of filter spectrophotometry and the associated problem of [N II] contamination of Hα fluxes versus 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 data cube is an asset.

From the observational point of view in the optical, no sharp distinction exists 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. Supernova remnants are defined by a line ratio larger than 0.4, although Blair & Long (1997) found that this criterion is washed out in some galaxies, for example, NGC 7793. 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. Figure 21 shows examples of a SNR and several superbubbles and shells. Although we recorded the presence of luminous blue stars within the nebulae whose stellar winds might contribute to the expansion of the shell, for example, 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 (Fig. 21).

Detailed studies of extragalactic SNR, including X-ray data, were for a long time restricted to the Magellanic Clouds. Recent advances, for example, 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, for example, 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).

Fig. 21. [S II] images of a supernova remnant and shells. Upper left panel: SNR a64, truncated at the northern edge of field (a); the object to the SW from the center is the optically thick H II region a10. Upper right panel: double-lobe shell d95/d96, hosting WR star. Lower left panel: double-ring shell d87/d88, reminiscent of the rings of SN1987A, however with a diameter of ∼10 pc. Lower right panel: shell b23, suspiciously aligned with the somewhat fainter shell B22.

4.6. Diffuse ionized gas

Results. Diffuse ionized gas, which 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. Figure 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⁻¹⁸ erg cm⁻² s⁻¹ arcsec⁻². 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 bottom left quadrant of field (b) to the top right quadrant of field (a). The interarm regions in the northern parts of field (e) and in field (j) are comparatively quiet. Still the minimum [S II] surface brightness levels in these areas amount to 2 × 10⁻¹⁸ erg cm⁻² s⁻¹ arcsec⁻² and 2 × 10⁻¹⁷ erg cm⁻² s⁻¹ arcsec⁻² in fields (d) and (j), respectively.

We 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 indeed find 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⁻¹⁷ erg cm⁻² s⁻¹ arcsec⁻² 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 of measuring radial velocities of the gas, we performed spectrophotometry with P3D over apertures of 5 × 5 spaxels, i.e., 1 arcsec², 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, which takes into account the line-of-sight velocity (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 scaled individually to the actual continuum level near each line for the wavelength regions around Hβ/[O III] and Hα/[N II]/[S II] since ULySS applies a polynomial fit to accommodate variations of the overall spectral energy distribution; this 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) stressed that ionization parameter mapping is effective at highlighting large, extremely faint structures. We applied the IPM technique exemplarily in field (a). As reported in Sect. 4.4, the [S II] surface brightness recorded in the narrowband images agrees generally well with spectrophotometry that was performed in selected areas. We chose to measure with P3D specifically (a) prominent filamentary DIG features and (b) local minima of low surface brightness regions 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², 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⁻¹, which is within the error identical with the systemic velocity of the galaxy, i.e., 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⁻¹, 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 SNR, 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 toward 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, for example, 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 and 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) 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, which 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, similar to Sabbadin et al. (1977), who distinguished shock-excitation gas in SNR from photoionized gas in H II regions. The resulting values of log(Hα/[SII]) versus log(Hα/[NII]) are plotted in Fig. 23, along with the envelopes for SNR and H II regions that were determined empirically by the latter authors. The red plot symbols refer to samples aD03 … aD09, aD15 … aD17, and aD19 … aD20, which are associated with filaments of enhanced Hα and [S II] emission, whereas the blue squares refer to the local minima of emission line surface brightness aD02, aD14, aD18; see Fig. 23. There is a clear separation between shock excitation in filaments (red circles) and photoionization in the undisturbed, low surface brightness regions (blue squares). It is interesting to note that aD14 is located within 20 pc of the giant H II region a66, i.e., relatively close to the photoion-izing star cluster.

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 (1985a,b), 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 photoionization alone could possibly be the relevant mechanism, owing to the patchy and turbulent nature of the ISM.

As a first test toward a more thorough future quantitative analysis, we 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 from shock excitation. As a result, we find a total Hα flux of 1.1 × 10⁻¹⁴ erg cm⁻² s⁻¹ for the filaments (5342 spaxels coadded) and a total of 8.8 × 10⁻¹⁵ erg cm⁻² 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⁻¹⁶ erg cm⁻² 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. Given the interest in understanding the nature of DIG, we are planning a future paper to extend this analysis to a study of the entire set of pointings in NGC 300.

Fig. 19. [S II] images reconstructed from MUSE data cubes for pointings (a), (b), …, (j). The color scale brings out the DIG more prominently than the corresponding Hα map, also suppressing many point sources.

Fig. 22. Ionization parameter mapping used to analyze the DIG in field (a). Left: [S II]/Hα, grayscale 0.15 … 0.8 (black–white). Right: [S II]/[O III], grayscale 1.5 … 10 (black–white).

Fig. 23. Left: continuum-subtracted [S II] 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).

Fig. 24. Hα image around the LBV candidate i60 (blue circle to the left) in field (i), covering 245 × 320 pc2. Red and yellow circles indicate H II regions and SNR, respectively. Contours at surface brightness levels of 1.5 and 2 × 10−17 erg cm−2 s−1 arcsec−2.

4.7. Kinematics

As highlighted in Sect. 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 color code that is explained in the legend, stretching from below 130 km s⁻¹ to above 190 km s⁻¹ in bins of 10 km s⁻¹. 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⁻¹ near the nucleus to around 200 km s⁻¹ at the extreme northwest. We 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. The latter authors 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 ± 2 km s⁻¹ 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 preferred a direct comparison to the Hα velocity map. The stellar LOSV follows well the Hα velocity field and has values around the systemic velocity close to the nucleus in field (a), out to 200 km s⁻¹ 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 behavior 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 also obtained ULySS fits for the local background spectra as an output from PampelMUSE for each data cube, there are also radial velocity maps for the unresolved stellar population. For the parameter setting of our PampelMUSE runs, we 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, for example, following the discussion of Falcón-Barroso et al. (2004) for the measurement of emission-free stellar kinematics. We reserve this step to the future analysis of the complete dataset with the pPXF tool (Cappellari & Emsellem 2004), which 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 color code as the Hα velocity map, where black areas indicate that no ULySS fit was obtained owing 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 color code as for the stars. We 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, even though we estimated completeness limits of m₅₀₀₇ = 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 are a number of remarkable outliers: two high velocity PNe in the nuclear region, i.e., (a03) with v_rad = 197 ± 9 km s⁻¹ and (b06) with v_rad = 244 ± 3 km s⁻¹, and at larger galactocentric distances four low velocity PNe, i.e., v_rad = 139 ± 11 km s⁻¹ (d63), 144 ± 7 km s⁻¹ (d79), v_rad = 145 ± 4 km s⁻¹ (e02), and 143 ± 3 km s⁻¹ (e20), respectively, and the extreme case of PN i84 with v_rad = 108 ± 8 km s⁻¹. 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 of 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.

Fig. 25. Kinematics in NGC 300. Panel a: radial velocities determined for individual stars; panel b: radial velocities obtained from unresolved background tiles produced by PampelMUSE; panel c: radial velocities determined for PNe, errors coded with plot symbols (full circles <10 km s−1, thick rings <20 km s−1, thin rings >20 km s−1); panel d: reference velocity field reproduced from Hlavacek-Larrondo et al. (2011). The velocity field color coding is similar to, but differs slightly from, the color code for point sources in panels a and c.

4.8. Serendipitous discoveries

As mentioned in Sect. 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 Sect. 4.3 that immediately became apparent upon inspection of the He II image of that region. This star was not specifically classified or earmarked 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 two-step procedure of imaging and follow-up spectroscopy, the discovery of symbiotic star candidates, and the detection of carbon stars.

Variable sources. Although the layout of our observations did not directly address 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 serendipitously made such a discovery. Upon the investigation of the apparent foreground star ID 32224, as described in Sect. 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, for example, 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), i.e., CL1 in Fig. 27. Nuclear star clusters are thought to grow dissipatively from in situ star formation (Milosavljević 2004) and/or by dis-sipationless 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 harbor 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 v_rad = 139 ± 2 km s⁻¹, 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 that of the nuclear cluster. The LOSV is determined as v_rad =119 ± 3 km s⁻¹.

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⁻¹, a value that is significantly lower than the average of 144 km s⁻¹ for the stellar population in this region, which is also the case for CL2. The presence of of a blue supergiant component and 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⁻¹, the cluster would have traveled 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 of 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 visu-alization tool, we incidentally noticed emission line signatures in the display of color-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 that have redshifts of z = 0.13 … 1.33, and even two candidates for Lyα emitting galaxies at redshifts z > 4. These objects are listed in Table 13. The IDs with an index a, b, … indicate systems with kinematically distinct components, which 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).

Fig. 26. Foreground star ID 32224 denoted A, from field (i). The nearby peculiar red star B (ID 1003, F475W = 24.45, F606W = 22.82, see Table 6), which 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.

Fig. 27. Reconstructed VRI image from MUSE data cube, showing the nuclear cluster CL1 and two newly discovered clusters in its vicinity. The FoV measures 245 × 245 pc2 (in projection).

5. 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, using 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 data cubes, in particular concerning the multitude of scientific results over a variety of different objects one can expect from such data.

Table 14 summarizes the various classes of objects for which individual spectra and images were obtained. The total number of objects, such as stars, PNe, H II regions, and SNR, including background galaxies, amounts to 2189. At this early stage of exploration of the new method of crowded field 3D spectroscopy, we decided to visually inspect each and every object, assigning quality flags and reassuring ourselves not to 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 Pampel-MUSE, 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 time consuming 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), which were obtained from a total of 24762 stellar spectra secured over ten 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 catalog, we found 82 blue emission lines stars, most of which show broad Hα emission lines, which are characteristic for evolved massive stars with strong stellar winds, for example, 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, which 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 accordance with the expectations.

7. Amongst emission line point sources we registered 36 bona fide PNe, 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 those of PNe. From previous PN surveys in NGC 300, we found several cases in which objects were misclassified as PN and better match the criteria of cHII. In total, we registered 51 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 that range between a few and more than 100 pc in size. Some of the SNR candidates exhibit very low surface brightness, and are therefore possibly confused 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⁻¹⁸ erg cm⁻² s⁻¹ arcsec⁻² 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 30 background galaxies with redshifts z = 0.13 … 1.33 were identified, which 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 are publicly available through the CDS and the MUSE website⁵.

6. Outlook

As a technical shortcoming at the present stage, we 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 were 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 separately modeled the 2D background surface brightness distribution in the continuum and in emission. For an instrument such as 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 our methods toward automatic data analysis to process the expected number of thousands of spectra more efficiently.

From the entire dataset, we intend to study the stellar population from the nuclear region to the northwestern spiral arm, including the leading and trailing interarm regions covered by our pointings, with the goal of inferring 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 such as migration, mergers, and runaway stars. Another challenging goal will be the determination of chemical abundances for individual stars; interesting prospects concern young and old stellar populations, for example, blue and red supergiants versus 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 acute with the revised picture of post-AGB evolution from the latest models of Miller Bertolami (2016), which 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 used in many studies.

With the hypothesis that DIG emission may be dominated by shock excitation, we are planning a follow-up paper to investigate systematically all of our pointings in more detail 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, which 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 only the format and performance of MUSE has enabled the efficient use of crowded field spectroscopy, both in terms of multiplex advantage and in using the PSF fitting technique to deblend overlapping stellar images. Combined with high instrumental throughput, the light collecting power of an 8 m 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).

Acknowledgments

The authors would like to thank the anonymous referee for helping to improve the quality of this paper. SK, PMW, and CS received funding through BMBF Verbundforschung (grants 05A14BAC, 05A14MGA, ERASMUS-F, grant 05A08BA1, and ELT-MOS, grant 05A15BA1). SD acknowledges BMBF support through grant 05A14MGA, and from DFG through DR281/35-1. AMI acknowledges support from the Spanish MINECO through project AYA2015-68217-P. MMR acknowledges fruitful discussions with PER and TJR. The authors are thankful to A. Vazdekis for providing the (unpublished) MIUSCAT library. Based on observations made with the NASA/ESA Hubble Space Telescope, and obtained from the Hubble Legacy Archive, which is a collaboration between the Space TelescopevScience Institute (STScI/NASA), the Space Telescope European Coordinating Facility (ST-ECF/ESA), and the Canadian Astronomy Data Centre (CADC/NRC/CSA). This research has made use of the Spanish Virtual Observatory (http://svo.cab.inta-csic.es) supported from the Spanish MINECO/FEDER through grant AyA2014-55216.

References

All Tables

All Figures

	Fig. 1. Footprint of MUSE pointings (a), (b), …, (j) in NGC 300, covering the nucleus (a), different fractions of the spiral arm that extends to the NW of the center (b), (c), (d), (e), (i), and the inter-arm region (j). Orientation: N up, E left. (Image credit: ESO).
In the text

Fig. 2. VRI images reconstructed from MUSE data cubes 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 color, 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.2 m WFI (credit: ESO).

In the text

	Fig. 3. MUSE vs. HST images of 4 × 4 arcsec2 size from field (i), as examples of the star in the center being affected by different severity of crowding. a/b: ID1154, A1 Ib star, unblended. c/d: ID1179, B3 III star, negligible blending. e/f: ID995, central star cluster in bright H II region, spectrum heavily affected by blending and nebular emission, no reliable classification possible.
In the text

	Fig. 4. Example spectra from field (i). Images in the two columns to the left are analogous to Fig. 3. The spectra cover the full extended MUSE wavelength range of 4600 … 9300 Å. Flux is plotted in units of 10−20 erg cm−2 s−1 Å−1. Black curves: observed spectra; orange: ULySS fits.
In the text

	Fig. 5. Deblending of two overlapping stellar images in field (i). The MUSE and HST maps to the left are analogous to Fig. 3, for units of spectra see Fig. 4.
In the text

	Fig. 6. Deblending of the extremely faint emission line star (A) from an unresolved group of K stars (C). For an explanation, see Fig. 5 and text below. The insert spectrum covers the wavelength range 6200 … 7000 Å, flux is plotted in units of erg cm−2 s−1 Å−1 (continuum background not subtracted).
In the text

	Fig. 7. Distribution of S/N as a function of F606W magnitude in field (i), selectively by quality flag (top), and Teff (bottom).
In the text

	Fig. 8. MUSE F606W photometry versus HST ACS magnitudes, broken down into cool stars (red dots) and hot stars (blue dots). Mean residuals against the 1:1 relation in 0.2 dex mag bins are shown as black dots. The error bars indicate the standard deviation while the gray horizontal line indicates zero.
In the text

	Fig. 9. Recovered Teff as a function of S/N for a total of 630 simulated spectra, as determined by ULySS. For details, see text.
In the text

	Fig. 10. Radial velocity uncertaintites ∆vrad plotted against S/N. The sample is segregated into hot and cool stars as in Fig. 7. The insert shows a histogram of the ∆vrad distribution, the 50th percentile at 20 km s−1.
In the text

Fig. 11. MUSE spectrum of A5 supergiant B11 from the sample of Bresolin et al. (2002), seen in field (b), ID25. This star has a V magnitude of 19.95. The MIUSCAT fit is plotted in orange with an offset for clarity. Wavelengths in units of Å, flux in units of 10−20 erg cm−2 s−1 Å−1. The reference wavelengths of a number of absorption lines (hydrogen Balmer and Paschen series, helium lines) are indicated.

In the text

	Fig. 12. MUSE spectrum of AIa star of magnitude mF606W = 20.78 in field (i), ID151. The MIUSCAT fit is plotted in orange with an offset for clarity. Wavelengths in units of Å, flux in units of 10−20 erg cm−2 s−1 Å−1.
In the text

	Fig. 13. MUSE spectrum of M1Ia star of magnitude mF606W = 21.56 in field (i), ID301. The MIUSCAT fit is plotted in orange with an offset for clarity. Wavelengths in units of Å, flux in units of 10−20 erg cm−2 s−1 Å−1. The reference wavelengths of absorption lines of calcium, magnesium, sodium, Hα, and Hβ are indicated.
In the text

	Fig. 14. Hα images reconstructed from MUSE data cubes for pointings (a), (b), …, (j). The background image is from the ESO/MPI 2.2 m WFI (credit: F. Bresolin). H II regions in common with the catalog of Deharveng et al. (1988) are indicated with red circles. Yellow circles denote the SNR recorded by Blair & Long (1997).
In the text

	Fig. 15. Left: examples for artificial PN point sources inserted in data cubes of fields (a) and (b), with seeing of 0.7″ and 1.2″, respectively. Right: fraction of artificial PN detections as a function of m5007, obtained with DAOPHOT FIND for fields (a), (b), (c), (i), (j), and seeing of 0.7″, 1.2″, 1.0″, 0.6″, 0.85″, respectively.
In the text

	Fig. 16. Comparison of PN and cHII magnitudes obtained with MUSE versus results from Soffner et al. (1996) in red and Peña et al. (2012) with blue plot symbols.
In the text

Fig. 20. Example spectra of PNe, cHII, SNR, and an emission line star in units of erg cm−2 s−1 per 1.25 Å spectral bin in the Hβ-[O III] and Hα regions, plotted against wavelength in Å. The spectra are offset by some arbitrary bias to better show the noise level. Panels a–d: PNe over the full range from the brightest to the faintest objects. The very faint object a03 with m5007 = 28.88 is an uncertain PN candidate, given that the [SII] line ratio, although critically dependent on diffuse nebular background subtraction, could also be indicative of a cHII, or even a SNR. Panels e, f: a bright and a faint cHII, and panels g, h: a bright and a faint SNR. Panel i: an emission line star is shown, presenting an Hα emission line width of 6.6 Å.

In the text

	Fig. 17. WR star in field (d). Top left: VRI map, candidate star indicated as object d38. Top right: Hα map, showing surrounding shells. Bottom: spectrum showing the He II feature (arrow), in arbitrary units.
In the text

Fig. 18. Ionization parameter mapping of two H II regions, from left to right: Hα, [S II]/[O III], VRI. The grayscale for the Hα maps shows high intensity in black, whereas the IPM maps feature high ratios in white. The region i97 (top row, 4.8″ × 5.8″) is optically thick, while i81 (bottom, 10.4″ × 11.8″) is apparently leaking ionizing radiation, a “blister” in the nomenclature of Pellegrini et al. (2012). A strange, [O III]-bright feature near the rim of i97 is indicated by an arrow (see text). The stars that are likely responsible for the ionization of the nebulae are visible as blue objects in the VRI images.

In the text

Fig. 21. [S II] images of a supernova remnant and shells. Upper left panel: SNR a64, truncated at the northern edge of field (a); the object to the SW from the center is the optically thick H II region a10. Upper right panel: double-lobe shell d95/d96, hosting WR star. Lower left panel: double-ring shell d87/d88, reminiscent of the rings of SN1987A, however with a diameter of ∼10 pc. Lower right panel: shell b23, suspiciously aligned with the somewhat fainter shell B22.

In the text

	Fig. 19. [S II] images reconstructed from MUSE data cubes for pointings (a), (b), …, (j). The color scale brings out the DIG more prominently than the corresponding Hα map, also suppressing many point sources.
In the text

	Fig. 22. Ionization parameter mapping used to analyze the DIG in field (a). Left: [S II]/Hα, grayscale 0.15 … 0.8 (black–white). Right: [S II]/[O III], grayscale 1.5 … 10 (black–white).
In the text

Fig. 23. Left: continuum-subtracted [S II] 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).

In the text

	Fig. 24. Hα image around the LBV candidate i60 (blue circle to the left) in field (i), covering 245 × 320 pc2. Red and yellow circles indicate H II regions and SNR, respectively. Contours at surface brightness levels of 1.5 and 2 × 10−17 erg cm−2 s−1 arcsec−2.
In the text

Fig. 25. Kinematics in NGC 300. Panel a: radial velocities determined for individual stars; panel b: radial velocities obtained from unresolved background tiles produced by PampelMUSE; panel c: radial velocities determined for PNe, errors coded with plot symbols (full circles <10 km s−1, thick rings <20 km s−1, thin rings >20 km s−1); panel d: reference velocity field reproduced from Hlavacek-Larrondo et al. (2011). The velocity field color coding is similar to, but differs slightly from, the color code for point sources in panels a and c.

In the text

	Fig. 26. Foreground star ID 32224 denoted A, from field (i). The nearby peculiar red star B (ID 1003, F475W = 24.45, F606W = 22.82, see Table 6), which 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.
In the text

	Fig. 27. Reconstructed VRI image from MUSE data cube, showing the nuclear cluster CL1 and two newly discovered clusters in its vicinity. The FoV measures 245 × 245 pc2 (in projection).
In the text