© ESO 2022

1. Introduction

The ‘baryon cycle’ – the collapse of gas to form stars and the subsequent re-injection of matter, energy, and momentum into the interstellar medium (ISM) – is an intrinsically multi-phase and multi-scale process. Flows of gas in and out of galaxies, as well as internal gas dynamics, connect the small-scale cycle of baryons with the larger galactic and circum-galactic scales. These processes drive both the evolution of galaxies in a large-scale cosmological context, and the still elusive small-scale physics involved in the collapse of gas cores, and the feedback from massive stars (Scannapieco et al. 2012; Haas et al. 2013; Hopkins et al. 2013; Agertz & Kravtsov 2015; Fujimoto et al. 2019).

Fig. 1. From kpc to ∼100 pc scale, red-green-blue (RGB) colour images (channels derived from reconstructed MUSE mosaic using the SDSS i, r, and g bands) of NGC 4303 (D = 17 ± 3 Mpc) at various spatial resolutions. The four columns, from left to right: at 64 pc (the homogenised resolution for our MUSE dataset) and then convolved to 250, 500, and 1000 pc. The size of the beam (FWHM) is provided as a hatched circle at the bottom-left corner of each top row panel. The bottom row shows zoomed-in views of the top row images (area shown as a white rectangle).

A key challenge for both observations and theoretical models is to connect the population of galaxies observed across cosmic time with the sub-parsec-scale physics of the star-formation process. Observationally, the cosmological context is addressed by large redshift surveys, while the small-scale physics can be most easily accessed within the Local Group or our own Milky Way. The physical processes associated with star formation are, however, also affected by varying local conditions (Kawamura et al. 2009; Colombo et al. 2014; Hughes et al. 2016; Egusa et al. 2018; Hirota et al. 2018; Sun et al. 2020b; Querejeta et al. 2021).

Nearby star-forming galaxies therefore offer a unique viewpoint at the interface of the cosmological and Galactic scales. Several classical studies have been dedicated to detailed multi-wavelength mapping of individual nearby targets (e.g., M51, M33, and M31; Kennicutt et al. 2003; Gil de Paz et al. 2007; Calzetti et al. 2005; Boquien et al. 2011; Viaene et al. 2014; Corbelli et al. 2017; Williams et al. 2018). We have, however, so far lacked a comprehensive multi-tracer campaign covering the entirety of the discs of a representative set of star-forming main-sequence galaxies (where the bulk of today’s stars are being formed; e.g., Brinchmann et al. 2004) down to their individual star-forming regions and probing the different phases of their ISM. The target scale is the typical inter-cloud distance scale of a few tens of parsec up to about 100 parsec: it represents the intermediate ‘structuring size’ or ‘cloud scale’ of star-forming galaxies, relating to individual gas clouds, clusters of young stellar objects, and discrete star-forming regions. Already at a few hundred parsec resolution, most of the structures associated with the gas clouds, stellar clusters, and dusty features are lost (see Fig. 1). When exploiting the nearby volume of galaxies up to about 20 Mpc, a scale of 100 pc typically requires arcsecond or sub-arcsecond full width at half maximum (FWHM) beams and, hence, relatively high spatial resolution supported by spectroscopic data over fields of view (FoVs) of a few arcminutes on disc galaxies.

Optical spectroscopy, in particular, is a powerful tool for placing the star-formation process in the context of its galactic host as it can probe the ionised ISM and the stellar backbone, which characterises the local disc environment. It more specifically provides key information pertaining to chemical abundances, stellar mass, star-formation histories, gas, and stellar kinematics. Spectroscopic observations at such scales over the large FoV required to map nearby galaxies have been very challenging due to the lack of suitable instrumentation.

The capabilities offered by integral field spectroscopy (IFS) have, however, greatly evolved over the last 30 years. Different hardware solutions, including fibres, micro-lenses, or advanced slicers, currently allow for a varied set of spatial and spectral samplings, filling factors, FoVs, and overall performance (see e.g., Bacon & Monnet 2017, and references therein). These and other spectroscopic mapping techniques have already begun to be applied to samples of nearby galaxies (e.g., PINGS, Rosales-Ortega et al. 2010; VENGA, Blanc et al. 2013; TYPHOON, Seibert et al., in prep.; SIGNALS, Rousseau-Nepton et al. 2019, CALIFA Sánchez et al. 2012, SAMI, Croom et al. 2012, MaNGA, Bundy et al. 2015). The Multi Unit Spectroscopic Explorer (MUSE) at the Very Large Telescope (VLT), with its coverage of most of the optical wavelength range and a FoV of 1 arcmin², properly sampling the seeing disc, has recently opened a significant new area of parameter space. In particular, the success and versatility of MUSE lies in a combination of factors, including its high overall performance (∼35% peak efficiency around 7000 Å; see e.g., MUSE/VLT User’s Manual), its multiplexing capabilities (90 000 spaxels, each with about ∼4000 spectral pixels), and, most importantly, its robust optical setup (based on advanced slicers and an industrial approach for the building of its 24 spectrographic units) and its dedicated advanced data reduction pipeline (Bacon et al. 2016; Weilbacher et al. 2020, see also Sect. 4.1). MUSE is one among just a few integral field units (IFUs) that can deliver a spectro-photometric view of the objects it targets, a key requirement for being able to robustly derive physical parameters from IFS data.

These combined characteristics make MUSE the ideal instrument for providing, for the first time, extensive mapping of nearby (D ≲ 20 Mpc) star-forming galaxy discs and resolving the mean inter-cloud distance in the ionised ISM (i.e. accessing the cloud scale). This paper presents the realisation of this ambitious goal in the form of the PHANGS-MUSE survey, an ESO Large Programme built on a VLT MUSE pilot project that mapped NGC0628 (Programme IDs: 1100.B-0651/PI: E. Schinnerer; 095.C-0473/PI: G. Blanc; and 094.C-0623/PI: K. Kreckel) to obtain spectrophotometric maps of the ionised gas and stars for 19 nearby galaxies.

The PHANGS-MUSE survey is a key part of the Physics at High Angular Resolution in Nearby Galaxies¹ (PHANGS) project. PHANGS aims at obtaining, for the first time, a comprehensive view of the star-formation process across different ISM phases in the range of environments present within a representative sample of nearby, massive, star-forming galaxies (see the detailed discussion about the sample in Leroy et al. 2021a). Key goals for PHANGS follow the following scientific threads: (a) infer the timescales of the star-formation process (i.e. molecular cloud lifetimes, feedback timescales, feedback outflow velocities, star-formation efficiencies, and mass loading factors); (b) quantify the importance of the various stellar feedback processes (i.e. ionising radiation, stellar winds, supernova explosions, etc.) in galactic discs, (c) resolve the chemical enrichment and mixing across galaxy discs (in both radial and azimuthal directions), and (d) establish how the clustering of young stars is seeded by and disrupts the structure of the ambient ISM.

To achieve these goals, observations must map and resolve the individual structures of the star-formation process (with sizes of a few pc to ∼100 pc; e.g., Sanders et al. 1985; Oey et al. 2003), such as (giant) molecular clouds, H II regions, and the resulting star clusters. Sampling the variety of environments (related to e.g., bars, spirals, centres, mass, and dynamical structures) present in nearby galaxies is needed to assess the environmental impact within and among different galaxy discs. In order to link our results to the larger-scale studies of galaxy populations, we have prepared a selection in accordance with the main sequence of star-forming galaxies (e.g., Brinchmann et al. 2004).

Fig. 2. Overview of large spectroscopic surveys of nearby galaxies. Left: large spectroscopic surveys in the plane defined by their spatial resolution (in physical units) and the number of spatial resolution elements surveyed. IFU surveys are shown with diamond symbols (green for those achieving cloud-scale or better resolution and blue if probing at kiloparsec scales). VENGA (Blanc et al. 2013) is not shown because it features fewer than 105 resolution elements. For reference, the single fibre SDSS survey (Abazajian et al. 2009, red triangle) is added. PHANGS-MUSE sits in the top-left of this space, ranking highly on both metrics. Right: large spectroscopic surveys in the plane defined by their spatial resolution (in physical units) and the number of galaxies surveyed. PHANGS-MUSE lies on the overall trend line of other IFU surveys.

Focusing on the optical spectroscopy aspect, PHANGS-MUSE covers a largely unexplored area of parameter space in terms of number of spectra versus spatial resolution compared to other state-of-the-art spectroscopic studies of nearby galaxies. In Fig. 2 (left) we compare PHANGS-MUSE with several other IFU surveys of local galaxies and with the Legacy (single 3″ fibre) Sloan Digital Sky Survey (SDSS Strauss et al. 2002; Abazajian et al. 2009). For each IFU survey, we estimate the number of independent spatial resolution elements as the ratio between the total area surveyed and the area of the point spread function (PSF) FWHM. In this parameter space, PHANGS-MUSE occupies a unique region, combining high spatial resolution with the largest number of spectral elements. PHANGS-MUSE resolves the galactic discs about 1.5 orders of magnitude better than large IFU surveys of the nearby Universe, such as CALIFA, SAMI, and MaNGA (Sánchez et al. 2012; Croom et al. 2012; Bundy et al. 2015), while delivering a factor of ∼2 more independent spatial resolution elements than MaNGA, the largest of these surveys. This comparison highlights the impressive information-gathering power of the MUSE instrument. It also helps contextualise the challenges associated with the processing of the PHANGS-MUSE dataset. PHANGS-MUSE complements the large kiloparsec-scale surveys, such CALIFA, SAMI, and MaNGA, which have observed ∼10³ − 10⁴ galaxies (Fig. 2, right), and accesses new physics by trading sample size for spatial resolution.

Three other MUSE surveys, the MUSE Atlas of Discs (MAD; Erroz-Ferrer et al. 2019), the Time Inference with MUSE Extragalactic Rings (TIMER; Gadotti et al. 2019), and GAs Stripping Phenomena in galaxies with MUSE (GASP; Poggianti et al. 2017), target nearby star-forming galaxies. TIMER focuses on the impact of bars and active galactic nuclei (AGN) on galaxy evolution, while GASP studies ongoing and past ram pressure stripping events: their samples are therefore focused on addressing specific science questions and are not representative of the population of star-formation main-sequence (SFMS) galaxies. The MAD survey, on the other hand, focused on main-sequence galaxies with log(M_⋆/M_⊙)> 8.5, selected to be nearby (z <  0.013, D <  55 Mpc) and moderately inclined (i <  70°), but obtained only one central MUSE pointing per object (and two-pointing mosaics in a few exceptional cases). MAD therefore probes only the inner regions (∼2 kpc) of nearby galaxies. For more distant targets, it samples a larger fraction of the galactic disc but at coarser spatial resolution (> 200 pc), starting to blend structures at the cloud scale (see Fig. 1). PHANGS-MUSE is complementary to all these surveys, covering the galactic discs of typical star-forming galaxies at 100 pc or better resolution.

In this paper we present the PHANGS-MUSE survey, providing both the global context for this campaign and information pertaining to the associated first public data release (DR1.0). We note that all released PHANGS-MUSE data are available via the ESO Archives (i.e. accessible programmatically using the ‘PHANGS’ data collection flag) as well as from the Canadian Astronomy Data Centre (CADC). We start by reviewing the top-level scientific goals of the PHANGS-MUSE programme and present the galaxy sample (Sect. 2). The MUSE observation strategy is described in Sect. 3. In Sects. 4 and 5 we detail the data reduction and data analysis pipelines we have developed specifically for this survey and provide further quality assessment in Sect. 6. In Sect. 7 we briefly describe the relevant internal and public data releases associated with the PHANGS-MUSE dataset. Finally, in Sect. 8 we provide a set of data-demonstration figures to illustrate the potential of the survey, and we present our conclusions in Sect. 9.

2. PHANGS-MUSE survey: Observational context, sample, and science goals

2.1. The PHANGS-MUSE galaxy sample

Fig. 3. PHANGS-MUSE sample in the M⋆ − SFR plane. Left: PHANGS sample compared with the population of local galaxies from z0MGS (Leroy et al. 2019, small grey dots). The large red circles represent the PHANGS-MUSE galaxies. We show the overlap with the ALMA (blue dots) and HST (black empty squares) components of the PHANGS project. The dashed line is the best fit to the SFMS from Leroy et al. (2019). Right: PHANGS-MUSE sample compared to two complementary projects, EDGE-CALIFA (Bolatto et al. 2017) and ALMaQUEST (Lin et al. 2020), which also target local galaxies with optical IFS and CO interferometric mapping. The dashed line is the best fit to the SFMS from Leroy et al. (2019) with associated scatter (grey shaded area).

The parent sample of the overall PHANGS programme was originally constructed according to the following four criteria (see details in Leroy et al. 2021a): (i) southern sky accessible, in order to be observable by ALMA and MUSE, with −75° ≤δ ≤ +25°; (ii) nearby (5 Mpc ≤ D ≤ 17 Mpc), in order to probe star-forming regions out to at least an effective radius in a reasonable timescale and simultaneously provide < 100 pc resolution; (iii) low to moderate inclination (i <  75°), to limit the effects of extinction and line-of-sight confusion and facilitate the identification of individual star-forming sites; and (iv) massive star-forming galaxies with log(M_⋆/M_⊙)≳9.75 and log(sSFR/yr⁻¹)≳ − 11. The mass cut is driven by the desired overlap with ALMA observations, which become increasingly expensive in the low-mass, low-metallicity regime.

These criteria ensure that individual molecular clouds and star-forming regions can be isolated without confusion while the selected galaxies are representative for galaxies where most of the star formation in the local Universe occurs (e.g., Brinchmann et al. 2004). The cuts do not strictly apply to the PHANGS sample described in Leroy et al. (2021a) because of subsequent revisions of the fiducial distance estimates (Anand et al. 2021), improved derivation of the mass-to-light ratios (M/Ls) and star-formation rates (SFRs; Leroy et al. 2019) and the incorporation of additional galaxies extending the sample in important directions. A subset of 90 galaxies in the PHANGS parent sample has been observed by the PHANGS-ALMA survey, to create wide-field (covering the actively star-forming disc, roughly 1−2 R_e), high-resolution (FWHM ∼ 1″) CO(2–1) maps resolving the molecular phase into individual molecular clouds (Leroy et al. 2021a).

As the MUSE effort started at the same time as the ALMA Large Programme, the target selection focused on the 19 targets that were already observed as part of the ALMA pilot projects, or had ALMA archival data of similar characteristics. More specifically, a MUSE pilot programme (PIs K. Kreckel and G. Blanc; see Kreckel et al. 2016, 2017, 2018) focused on a close nearby face-on target, namely NGC 628, which was then followed by 16 galaxies observed as part of the ALMA campaign to probe the SFMS, further complemented with targets present in the ESO archive: this led to a sample of 19 nearby systems. Key properties of the PHANGS-MUSE targets are summarised in Table 1, and the distribution (PHANGS-MUSE objects in red) in the SFR versus stellar mass, M_⋆, plane is shown in Fig. 3 (left), both in relation to the other PHANGS surveys – ALMA and the Hubble Space Telescope (HST) – and relative to the main sequence of local (D <  50 Mpc) star-forming galaxies, as derived by Leroy et al. (2019) via a joint UV+IR analysis. Our sample covers a wide stellar mass range (9.4 <  log(M_⋆/M_⊙)< 11.0), but is, biased towards high masses (median stellar mass is log(M_⋆/M_⊙)=10.52), and towards the upper envelope of the SFMS – the median SFMS offset is +0.21 dex with respect to the z = 0 Multi-wavelength Galaxy Synthesis (z0MGS) relation (Leroy et al. 2019) – due to the need to adopt early, but sometimes uncertain, measurements of distances and SFR for the full PHANGS sample of star-forming main-sequence galaxies (see details in Leroy et al. 2021a). The PHANGS-MUSE sample does not include any of the Green Valley targets from PHANGS-ALMA and PHANGS-HST. Passive galaxies are excluded from the PHANGS sample by design, although a few passive galaxies have been targeted by PHANGS-ALMA as part of an ancillary programme.

It is useful to compare PHANGS-MUSE with other programmes aiming at obtaining both molecular gas and optical IFU spectroscopy of local galaxies. EDGE-CALIFA (Bolatto et al. 2017) and ALMaQUEST (Lin et al. 2020), consisting of 125 and 46 galaxies, respectively, are the only comparable efforts in this category. These surveys provide a more uniform sampling of the main sequence, and extend to the Green Valley, as shown in Fig. 3 (right). Unlike PHANGS, however, they both observe galaxies at approximately kiloparsec resolution, insufficient to resolve the physics of star formation on the scale of individual clouds. The MUSE targets within the PHANGS sample further have a large set of ancillary data on resolved scales, as emphasised in Sect. 2.2.

2.2. PHANGS-MUSE in the multi-wavelength context

The PHANGS programme is built on three main pillars. In addition to the PHANGS-MUSE programme, PHANGS leverages a Large Programme imaging the cold molecular phase in CO(2–1) with ALMA (Leroy et al. 2021a, 2021b, PI: E. Schinnerer), and a high-resolution Legacy survey of star clusters and stellar populations using five-band NUV-U-B-V-I imaging with the HST (Lee et al. 2022, PI: J. Lee). A fourth pillar is expected in the coming years, as we have also been awarded a James Webb Space Telescope (JWST) Treasury programme (PI: J. Lee) to image the PHANGS-MUSE sample. This practically means that the PHANGS-MUSE sample of 19 galaxies will ultimately have the full suite of MUSE, ALMA, HST, and JWST data. In Fig. 4 we give an example of the combined MUSE, ALMA and HST coverage for one of the galaxies in our sample, NGC 3351. In Fig. 5 we provide an illustrative view on the synergy between PHANGS-MUSE, PHANGS-ALMA and PHANGS-HST datasets (top left panels), with the specific power of optical spectroscopy allowed by MUSE, leading superb constraints on stellar and gas kinematics, the distribution and properties of the ionised gas and stellar populations.

Fig. 4. Synoptic view of the PHANGS multi-wavelength data using NGC 3351 for illustration. Left: blue and red contours showing the footprints of the PHANGS-ALMA and PHANGS-MUSE data. Within the respective footprints, we show flux maps for Hα from MUSE (light red) and CO(2–1) from ALMA (light blue). The footprint of the HST observations included in the PHANGS-HST data release is also shown as an orange contour. Right: HST imaging (F814W filter) shown in colour. The image covers a larger FoV with respect to the left panel in order to show the entire area imaged by HST. We also indicate several radial metrics: the disc scale length (Rd, as derived in Leroy et al. 2021a) and R25. The ALMA and MUSE footprints are also shown, same as in the left panel. We note that all MUSE, ALMA, and HST footprints of the 19 galaxies can be found at https://archive.stsci.edu/hlsp/phangs-hst.

Fig. 5. Multi-wavelength, multi-phase view of NGC 4535. The top panels present the stellar and gas distribution and kinematics. Top, second panel from the left: multi-emission line view (Hα in red, [O III] in blue, [S II] in green) tracing the sites of massive star formation along the spiral pattern, with differences in the combined colours highlighting the changes in local physical condition (e.g., abundance, ionisation parameter) and ionising source. Second panel from the left: extracted spectra (marked as white circles) demonstrating the typical characteristics of four (numbered) regions, namely: 1- dominated by stellar continuum (red); 2- AGN (orange); 3- H II regions (light purple), and 4- supernova remnants (SNR; purple). AGN line emission is superimposed on the strong stellar continuum absorption, with distinctive strong [O III] and [N II] emission. In the outer disc, H II regions and SNRs have less contribution from the stellar continuum. Broadened line shapes are apparent in the expanding SNR, in contrast to the narrower H II region line emission, and show strong [S II] and [N II] relative to Hα. Zooming into one section of the spiral arm (white box), spatial offsets between the HST star clusters and ionised gas (left top) and between the ionised gas and ALMA molecular gas (left bottom) demonstrate a time sequence evolution across the spiral pattern. Centre right: the gas and stellar velocity fields, mapped through the MUSE spectroscopy, highlight deviations from regular rotation and indicate dynamically driven radial flows along the spiral and bar structures. Right: the stellar mass surface density (ΣM⋆; see Sect. 5.2.4) highlights the location of these dynamical spiral and bar structures, and provides crucial constraints on the underlying gravitational potential affecting all stellar and gaseous processes in the disc.

This core observational effort is supplemented by a suite of complementary data, including ground-based narrow-band imaging (PHANGS-Hα; Razza et al., in prep., PIs G. Blanc and I-T. Ho), Keck Cosmic Web Imager (KCWI) spectroscopy (PI: K. Sandstrom), Canada-France-Hawaii Telescope SITELLE [O II] imaging (PI: A. Hughes), Russian 6m Fabry-Perot Interferometre spectroscopy (PI: E. Egorov), atomic hydrogen 21 cm mapping (Utomo et al., in prep., PI: D. Utomo), far-UV imaging with AstroSAT (Rosolowsky et al., in prep., PI: E. Rosolowsky), stellar mass maps with corresponding environmental masks (Sheth et al. 2010; Querejeta et al. 2015, 2021), dust maps obtained from archival Spitzer and Herschel imaging (Kennicutt et al. 2003, 2011; Clark et al. 2018) as reprocessed by Chastenet et al. (in prep.), and maps of molecular dense-gas tracers (Jiménez-Donaire et al. 2019). Tailored numerical work aims to provide the required reference simulations (e.g., Jeffreson et al. 2020; Utreras et al. 2020).

2.3. PHANGS-MUSE survey science goals

The MUSE observations of our selected nearby galaxies provide observational constraints on the structures (i.e. spirals, bars, centres) that make up the galaxy discs through various measurements by both covering various hosts and spatially resolving those structures. That includes the identification and spectroscopy of individual H II regions, their spatial distributions, brightnesses, metallicities, ionisation properties, and for a subset of regions even measurement of weaker temperature-sensitive lines. Those observations also provide detailed optical coverage of the stellar populations, constraining a two-dimensional view of the star-formation history and stellar mass distribution. They represent a unique probe of the gaseous and stellar dynamics via resolved kinematics; and additional components such as dust (extinction via the stellar continuum and Balmer decrement) or AGN (broad lines or high-ionisation emission lines). The PHANGS-MUSE dataset will more specifically inform several key science goals, which we now briefly review in turn.

A local perspective on scaling relations. Scaling relations have often guided observational and theoretical work towards our understanding of star-formation-related processes in galaxies (Kennicutt 1998; Brinchmann et al. 2004; Salim et al. 2007; Bigiel et al. 2008; Blanc et al. 2009; Leroy et al. 2013; Cano-Díaz et al. 2016; Hsieh et al. 2017; Medling et al. 2018; Lin et al. 2020; Sánchez et al. 2021; Ellison et al. 2021; Pessa et al. 2021; Querejeta et al. 2021). Basic properties such as stellar, molecular gas, or total gas surface density and SFR surface density have been probed to anchor representative timescales, such as how long it takes to deplete available gas reservoirs on galactic scales (Schmidt 1959; Kennicutt 1989, 1998; Saintonge et al. 2011, 2017; Cicone et al. 2017), as well as more locally on sub-kiloparsec scales (Wong & Blitz 2002; Bigiel et al. 2008; Leroy et al. 2008, 2013; Genzel et al. 2010; Schruba et al. 2011; Momose et al. 2013; Bolatto et al. 2017; Lin et al. 2020; Sorai et al. 2019).

The PHANGS-MUSE dataset provides robust constraints on the SFR, through direct Hα flux measurements corrected for extinction (via the Balmer decrement) and contribution from the diffuse gas component. The MUSE spectral coverage also probes gas-phase metallicity tracers, and key age and metallicity-sensitive stellar continuum spectral features, therefore allowing accurate derivations of M/Ls and stellar mass surface densities. The PHANGS dataset allows us to probe scaling relations at the global, kiloparsec-scale or cloud-scale levels, and to investigate the effect of different galactic environments (pressure budget, morphological, dynamical; see Sect. 8.1). Understanding how scaling relations vary as a function of scale and galactic environment (e.g., Pessa et al. 2021) will shine new light on the driving mechanisms for the observed trends, as well as the source of the associated scatter.

The impact of stellar feedback, in relation to local and global environments. All theories and simulations now agree that stellar feedback plays a central role in the self-regulation of the star-formation process across a wide range of galactic environments (e.g., Mac Low & Klessen 2004; McKee & Ostriker 2007; Ostriker et al. 2010; Hopkins et al. 2014; Agertz & Kravtsov 2015, 2016; Grisdale et al. 2017; Semenov et al. 2018, 2021; Fujimoto et al. 2019). Yet stellar feedback comes in many forms: radiative ionisation and heating, radiation pressure, stellar winds, and supernova explosions. These processes affect not just the local (< 100 pc) surroundings but can impact on kiloparsec scales and contribute to the pervasive diffuse ionised gas (DIG) observed throughout spiral galaxies (Zurita et al. 2000; Haffner et al. 2009; Zhang et al. 2017). While simulations have made progress considering the combined effects of these feedback processes (e.g., Hopkins et al. 2014; Rathjen et al. 2021), observations have lagged behind. Only a few nearby targets, including our Milky Way (e.g., Barnes et al. 2020; Olivier et al. 2021), the Magellanic Clouds (Pellegrini et al. 2011; Lopez et al. 2011, 2014) or NGC 300 (McLeod et al. 2020), have seen a careful inventory of the relative strength and location of sources of stellar feedback (see also Chevance et al. 2022; Barnes et al. 2021).

The PHANGS-MUSE survey of nearby galaxies will enable the quantitative study of the different forms of stellar feedback (radiative and mechanical) across galactic discs. In combination with the ALMA CO maps, the MUSE data can probe the interactions (such as localisation, dynamical and pressure differences) between the warm (10⁴ K) and the cold (< 100 K) gas reservoir on local and global scales, and potential variations with key galactic parameters. The data can provide measurements of the local balance of input momentum and energy from radiation, stellar winds, and supernovae constrained via the HST-derived massive stars and cluster catalogues (Turner et al. 2021, Whitmore et al., in prep., Larson et al., in prep.), MUSE-derived H II region properties, and the MUSE-modelled star-formation histories against the restoring forces of gas self-gravity and stellar gravity (derived from the ALMA molecular gas maps and the MUSE stellar mass maps) at a succession of spatial scales (see e.g., Sun et al. 2020a; Barnes et al. 2021). We can furthermore model the escape of radiation from individual H II regions and quantify its contribution to the ionisation of the kiloparsec-scale DIG, in combination with hot evolved low-mass stars and other ionising sources (Belfiore et al. 2022). This will give us a local assessment of the impact of individual feedback processes from the scale of individual regions to large parts of galaxies.

Quantifying the chemical enrichment and mixing in galactic discs. Radial metallicity trends have been observed in galaxy discs for decades, and more recently quantified in the overall population of local galaxies by large IFU surveys (i.e. CALIFA, Sánchez et al. 2014; SAMI, Croom et al. 2021; MaNGA, Bundy et al. 2015). Going beyond the radial trends, measurements of azimuthal variations and small-scale patterns remain poorly constrained, while being crucial to understand the key processes driving the chemical evolution of the ISM (e.g., Zaritsky et al. 1994; Sánchez et al. 2014; Belfiore et al. 2017), flows of gas (pristine or enriched), and the redistribution of metals from their birth sites to kilo-parsec scales. Tantalising evidence of azimuthal variations in gas-phase oxygen abundance have been obtained by high spatial resolution IFU studies (Sánchez-Menguiano et al. 2016; Vogt et al. 2017; Kreckel et al. 2019) and multi-slit spectroscopy (Berg et al. 2015; Croxall et al. 2016). Ho et al. (2017, 2018), for example, observed clear azimuthal metallicity variations associated with the spiral arms of NGC 1365 and NGC 2997 in ∼100 pc resolution pseudo-IFU long-slit data. Yet, it the origin of these azimuthal metallicity variations remains unclear; it is unclear if they are driven by localised self-enrichment and spiral-arm-induced mixing (Ho et al. 2017) or by radial flows in the disc (Sánchez-Menguiano et al. 2016).

Establishing the physical meaning of such variations requires high spatial resolution IFS observations of nearby galaxies, reaching out beyond the brightest H II regions, minimising the contamination by the DIG. It also requires probing various emission lines, isolating individual H II regions and inventorying them. The PHANGS-MUSE data are able to derive measurements (both from strong-line calibrations, as in Kreckel et al. 2019, and from direct electron temperature T_e determinations, as in Ho et al. 2019) of such patterns and their relationships to bars and spiral arms, the key drivers of radial flows in galaxies. By linking such variations with a detailed analysis of bar and spiral arm pattern speeds (Williams et al. 2021) and gas flows through our galaxies based on combined MUSE and ALMA data, we will be able to determine the key driver of mixing within galactic discs.

The role of dynamical regimes on the triggering, boosting or inhibiting of star formation. Dynamical environments play a key role in the redistribution of the gas reservoir and in setting the efficiency of star formation within discs. Bars, spirals, rings, resonances and central regions (e.g., Verley et al. 2007; Sanchez-Blazquez et al. 2011; Meidt et al. 2013; Renaud et al. 2015; Sun et al. 2020b; Kretschmer & Teyssier 2020; Gensior et al. 2020; Henshaw et al. 2020b) are characterised by different regimes associated with, for example, shear, torques, instabilities, gas flows, compression, or shocks. The PHANGS-MUSE dataset will help us constrain the local star-formation history via spectral fitting techniques, to thus characterise the stellar mass contribution. It will also provide unique leverage on the gravitational potential via the mapping of stellar and gas kinematics, and its various tracers (ionised gas, stellar populations; Kalinova et al. 2017; Leung et al. 2018; Bryant et al. 2019; Shetty et al. 2020). The determination of the star-formation history of the stellar disc and its link with the underlying dynamical orbital structure will provide a key constraint to understand the assembly and evolution of stellar discs. At the same time, we will compare gas and stellar surface densities and kinematics to predictions from equilibrium disc models to assess the scale at which vertical equilibrium (e.g., Ostriker et al. 2010; Ostriker & Shetty 2011) and radial disc stability (e.g., Hunter et al. 1998; Martin & Kennicutt 2001; Krumholz et al. 2018; Romeo 2020) emerge, balancing stellar feedback and gravity. In such a context, the synergy with numerical (hydro-dynamical) simulations will be paramount to further probe the relevant processes and their respective timescales (Utreras et al. 2020; Fujimoto et al. 2019; Jeffreson et al. 2020).

A multi-purpose legacy dataset. The sensitivity and physical resolution of the PHANGS-MUSE data will additionally allow for further investigation in a number of areas, including, for example, precise distance determination via planetary nebula luminosity functions (Kreckel et al. 2017, Scheuermann et al., in prep.), and identification of supernova remnants (SNRs) via line ratio diagnostics and line broadening (see e.g., Kopsacheili et al. 2020, and references therein). Our science goals are highly complementary to existing kiloparsec resolution IFU studies of hundreds (CALIFA Sánchez et al. 2012, CALIFA) or thousands of galaxies (SAMI, Croom et al. 2012; MaNGA, Bundy et al. 2015) and to future imaging spectroscopy of high-redshift targets, for instance using JWST and the Extremely Large Telescope (ELT). The PHANGS-MUSE dataset will serve calibration purposes and act as a training sample to enable accurate physical parameter estimation from lower-resolution observations of massive main-sequence star-forming galaxies. Our overall goal is that PHANGS-MUSE becomes a long-standing legacy dataset, providing the community with a reference for future observational, theoretical and simulation works.

3. Observations

3.1. Observing strategy

The PHANGS-MUSE survey covers a substantial fraction of the star formation within the galactic discs for a sample of 19 nearby star-forming spiral galaxies (see Sect. 2.1). The footprint of the PHANGS-MUSE survey is shown in Fig. 6: it was designed to overlap with the area of sky imaged in CO(2–1) by PHANGS-ALMA, which was itself aimed at encompassing all regions of active star formation inside the disc, including on average 70% of the WISE3 luminosity across the PHANGS-ALMA sample (see Leroy et al. 2021a). We designed the mosaics with a preference for a north-south orientation for the individual MUSE pointings. In a few cases, we rotated the position angles of the pointings to best cover the area of the galaxy we wished to map (e.g., NGC 1672 or NGC 7496).

Fig. 6. Footprints for the MUSE observations of PHANGS galaxies. Each panel represents one target of the PHANGS-MUSE sample, with a 5 × 5 arcmin2 FoV from the WFI Rc-band images (r-band du Pont for NGC 7496), and the footprints of the MUSE exposures are overlaid in red. Pointings marked with the ⌀ symbol (in NGC 1365, NGC 1512, NGC 1566, NGC 2835, and NGC 3351) and outlined in blue correspond to observations acquired outside of the main PHANGS campaign but reduced following the same data flow and released as part of PHANGS-MUSE. The vertical white bar on the left side of each panel indicates a scale of 5 kpc.

Individual galaxies are covered by a variable number of pointings depending on their angular size, ranging from 3 (e.g., NGC 7496) to 15 (e.g., NGC 1433), for a total of 168 individual pointings (5 of which were obtained from the ESO archive and had been observed by other programmes). Pointings were placed to have an overlap of 2″ (∼two resolution elements, and ∼10 MUSE spaxels) between adjacent fields to assist in the alignment process for the final mosaics. We further optimised the pointings’ gridding, to be as regular as possible to achieve our planned coverage, also for those galaxies that included archival pointings.

With the PHANGS-MUSE campaign, we aimed at detecting the Hβ line (at an S/N >  5) on top of the galaxy’s stellar continuuum (at S/N >  10). The existing ALMA coverage typically extends to I-band surface brightnesses of ∼22 mag arcsec⁻². This initially translated to a total of 43 min exposure on target (for a dark night, 7 days from the moon, and a typical airmass of 1.25), assuming a 5σ detection for Hβ on individual pointings, a 10σ detection on the stellar continuum, and about 20σ detection on Hα in the interarm regions (assuming a minimum binning of 5 × 5 spaxels, or 1″²). Considering the significant spatial variation in stellar populations, extinction, ionised gas content and regimes, this was intended to be a simplified and broad observational strategy. As detailed in Sect. 6.2.2, we typically detect Hβ in 50–80% of the individual spaxels throughout our sample (see Fig. 20).

Fig. 20. Fraction of 0.2″ spaxels inside of 0.5 R25 that have 3σ detections above a given surface brightness threshold (SB) for a representative sample of emission lines. Hα is detected in upwards of 95% of all pixels in most galaxies, with a lower ∼80% fraction in some strongly barred systems (e.g., NGC 1300, NGC 1433, and NGC 1512). Typical 3σ flux sensitivity in Hα is 4−7 × 1037 erg s−1 kpc−2 (3−7 × 10−19 erg s−1 cm−2 per 0.2″ spaxel). Hβ is typically detected in 50−80% of spaxels. Low-ionisation lines ([N II]6584, [S II]6717, and [S II]6731) are detected in 60−95% of spaxels, while the high-ionisation [O III]5007 line emission is less common (50−70% of pixels). For contrast, the faint auroral [N II]5754 emission line is detected at a 3σ level in ∼5% of spaxels (though fewer than 1% are detected at 5σ).

For our Large Programme, each pointing was observed at four different orientations separated by 90 degrees to mitigate the impact of individual MUSE slicers and the effects of the instrumental line spread function (LSF). These exposures were intertwined with two dedicated offset sky pointings, following an Object(O)-Sky(S)-O-O-S-O pattern. For each pointing the total on-source and sky integration times are 43 and 4 min, respectively. The entire PHANGS-MUSE Large Programme alone consists of a total telescope time of ∼172 hours. Observations were carried out using MUSE wide field mode (WFM), using the nominal (non-extended) wavelength range, either in seeing-limited (WFM-noAO) for most of the early acquired data or ground-layer adaptive optics (WFM-AO) mode. The observing campaign was designed for and started with the noAO (without the use of ground-layer adaptive optics) mode, the only mode available at the time. The WFM-AO mode for MUSE was offered in ESO Period 101, one semester after the start of the PHANGS observing runs. After some testing on early PHANGS datasets in terms of the extraction of stellar populations and emission lines information, and considering the small overheads incurred by the AO mode, we decided to implement the systematic use of the WFM-AO mode in 2018 for all targets that had not then yet been started (thus securing a single setup per mosaic). The 168 pointings result in about 15 Million spectra, covering a wavelength range of [4750−9350 Å], with the spectral resolution (σ) going from about 80 km s⁻¹ (at the blue end) to 35 km s⁻¹ (at the red end).

A pre-existing set of archival pointings, targeting the centres of some of our target galaxies, were not re-observed, but were reduced together with all other pointings, and included in our mosaics. The central pointing of NGC 3351 was acquired in the course of the MAD survey (Erroz-Ferrer et al. 2019), the central pointing of NGC 1365 was taken from the Measuring Active Galactic Nuclei Under the MUSE Microscope Survey (MAGNUM; Venturi et al. 2018) of local AGN, while the central fields for NGC 1512, NGC 1566 and NGC 2835 were observed as part of the TIMER project (Gadotti et al. 2019). The observations of NGC 0628 were obtained by the PHANGS collaboration as part of a dedicated pilot programme (PIs K. Kreckel and G. Blanc; see Kreckel et al. 2016, 2017, 2018). Overall our mosaics include 168 MUSE pointings: twelve pointings of the pilot target NGC 0628, 151 pointings executed as part of the PHANGS-MUSE Large Programme, and five archival pointings.

Table A.1 summarises, for each target, the main properties of the pointings such as the sky coordinates, day and time of execution, number of exposures, PSF, and observing mode. Occasionally, due to weather-related or technical issues during the observations, one pointing had to be split into two different observing blocks (OBs). The observing strategy for the Large Programme consists of a nominal four exposures per pointing. For the prototype NGC 0628, data for each pointing was observed with three exposures instead. For some pointings, we had to discard one or two exposures due to high sky brightness or extreme variability. For some pointings, on the other hand, we obtained one extra exposure, mainly as a result of OB splitting, as mentioned above. Due to a technical problem during the observations of NGC 1385, one of the pointings (i.e. P03) was observed twice and, therefore, consists of eight science exposures.

Fig. 7. Schematic data flow representing the PHANGS-MUSE data reduction part (DRP) of the pipeline. pymusepipe is used as the global process wrapper. The data flow includes a data organiser that addresses the raw data files, followed by standard MUSE data reduction steps (via MUSE DRS), which leads to the first reduced ‘Pixel Tables’ (PixTables) and cubes. These are aligned using WFI or Direct CCD reference images (via reproject). Aligned individual PixTables are then either resampled on a common grid or directly mosaicked. Optional convolution is performed using the framework from mpdaf and astropy, while pypher provides the required kernel cube.

3.2. Ancillary wide-field imaging data

In several steps of the reduction and characterisation of the MUSE data (e.g., exposure alignment, PSF measurement), it was necessary to compare our products to reference wide-field imaging data. For this purpose, we used R_c-band images obtained by the PHANGS collaboration with the Wide Field Imager (WFI; Baade et al. 1999) on the La Silla’s 2.2m MPG/ESO telescope, and r-band imaging from the Direct CCD camera of the 100 inch Las Campanas du Pont telescope. This imaging data were taken as part of the PHANGS-Hα survey, aiming to obtain narrow-band continuum-subtracted Hα maps for all PHANGS-ALMA galaxies (Razza et al., in prep.). Eighteen of the PHANGS-MUSE galaxies have been observed with WFI, while NGC 7496 has been observed with Direct CCD on the du Pont telescope (see Fig. 6 for examples of these data).

Here we briefly summarise the relevant processing steps for the R-band imaging data. A more detailed description of the observations and reduction steps can be found in Razza et al. (in prep.). For the R-band imaging we observe each field with a total integration time ranging from 900 to 1200 s, and perform a standard reduction. The imaging FoV (WFI: 34′ × 33′, Direct CCD: 8.85′ × 8.85′) is much larger than the MUSE FoV, and also large enough to allow for a robust astrometric and photometric calibration via multiple bright stars.

Each exposure was astrometrised independently and re-projected into a common grid. Photometric calibration was achieved by performing aperture photometry on a sample of stars and comparing to Gaia DR2 magnitudes (Gaia Collaboration 2018; Riello et al. 2018) converted to Johnson-Cousin (WFI) and Gunn (Direct CCD) broadband magnitudes via the appropriate conversions (Evans et al. 2018). A two-dimensional background was fitted after masking all the sources in the field (galaxy included), and subtracted from all the exposures, which are then combined to obtain the final frames. Typical seeing for the R-band observations of the PHANGS-MUSE sample was 0.8″, providing a close match to the typical MUSE resolution (listed in Table 1).

4. Data reduction

MUSE itself is a monolithic instrument consisting of 24 IFUs, each IFU containing an image slicer and spectrograph. The removal of the instrument signature is efficiently addressed via the excellent and advanced MUSE data processing pipeline software (MUSE DRS, hereafter) developed within the MUSE consortium (Weilbacher et al. 2020). In order to address, organise, reduce and analyse such a large dataset, we developed a dedicated framework, taking advantage of individual pieces of software or packages, which we describe in the following sections. Figure 7 represents a schematic of the data flow for the data reduction part (DRP) of the PHANGS-MUSE campaign, also flagging the usage of specific packages, with pymusepipe² providing the overall wrapper around these recipes.

4.1. Framework and software environment

We aimed at an almost fully automated pipeline that could be easily tuned to specific needs and potential changes associated with the survey and science goals. A number of constraints were expressed early in the project to satisfy that need, as well as the ability to easily add new acquired datasets, or rerun the entire data flow several times.

This required a modular approach, based on a set of robust packages including the MUSE processing pipeline itself. We adopted Python as the main coding and scripting language, exploiting its object-oriented capabilities, its wide community support and the existence of both robust libraries and packages. We now briefly describe the main packages and framework developed or used for the PHANGS-MUSE data.

pymusepipe. The DRP is controlled by a newly developed and dedicated Python package pymusepipe, which serves as a wrapper around the main processing steps of the data reduction. pymusepipe includes a simple data organiser and prescriptions for the structure of the data files (but no database per se), a wrapper around the main functionalities of MUSE DRS, accessed via EsoRex command-line recipes, to remove the instrumental signatures. pymusepipe additionally provides a set of modules supporting the alignment, mosaicking, (two-dimensional and three-dimensional) convolution, and data quality processing. Some of the details pertaining to each module or process are described in the next sections. pymusepipe is thus at the core of the DRP we apply to the MUSE data.

MUSE DRS. The MUSE spectrograph delivers ∼90, 000 spectra, each of about 4000 pixels covering most of the optical wavelength range, over a contiguous FoV of about 1′ × 1′. The raw fits MUSE data thus reflect, via its 24 extensions, the size and shaping of information ending up on the 24 individual detectors. The goal for MUSE DRS is to address such a complex data and science format, remove the instrument signature, combine and resample exposures for further scientific usage. MUSE DRS, developed by the MUSE team (Weilbacher et al. 2020), thus represents a pillar of any data reduction dealing with MUSE data, and follows an approach that minimises the need for resampling steps, using a table-based (PixTable) representation of the data. PixTables encompass the exact origin of the signal on the charge-coupled device and its associated IFU and slice identities. These PixTables can be projected as data cubes onto a given three-dimensional (sky positions and wavelength) grid using given geometric, astrometric and calibration information as derived with MUSE DRS. We are using the latest version of the MUSE DRS available at the time of writing (v2.8.3-1).

mpdaf. The MUSE processed data are natively derived in a PixTable format, which can be projected onto regular three-dimensional data cubes. The PixTables and cubes can themselves be used to reconstruct images in specific filters or extract spectra. mpdaf (Bacon et al. 2016) is a Python package providing an efficient and user-friendly framework to address such data in a transparent way. mpdaf is thus an important component of pymusepipe where existing mpdafPython classes have been complemented with utility functions (e.g., alignment, convolution, or image reconstruction).

pypher. The final MUSE mosaics, built either directly from PixTables (via muse-scipost or muse-exp_combine in MUSE DRS) or from already aligned and resampled data cubes (see Fig. 7), have PSFs that are varying over the spatial FoV and spectral range (Sect. 4.2.6). pypher³ (Boucaud et al. 2016) provides a robust tool to derive kernel cubes feeding a fast-Fourier-transform-based convolution algorithm to homogenise the end-product MUSE data cubes. Given two arbitrary PSF images, the pypher software uses a Wiener filter with a regularisation parameter to compute the convolution kernel needed to move from the input PSF to the output one. The power of such an algorithm is its applicability to general PSFs, expressed analytically or not. We used pypher to move from the wavelength-dependent circular Moffat PSF typical of the MUSE spectrograph, to a wavelength-independent circular Gaussian. The details about the characterisation of the MUSE PSF and the convolution of the data cubes are presented in Sects. 4.2.6 and 4.2.8.

GeneralPythonframework. We make use of generic but powerful Python packages, including numpy (Harris et al. 2020), scipy (Virtanen et al. 2020), matplotlib (Hunter 2007), and most importantly astropy (Astropy Collaboration 2018) for its excellent interface with FITS files (astropy.io.fits), and some specific modules including data management tools (e.g., astropy.tables) and units (astropy.units).

4.2. The data reduction work flow

In the following, we provide additional details regarding key steps of the MUSE data reduction work flow (DRP, Fig. 7).

4.2.1. Data organiser

The pymusepipe data organiser relies on reading existing FITS files in a user-defined directory via configuration files. All (compressed or uncompressed) FITS files are scanned and searched for specific keywords (e.g., OBJECT, TYPE, DATE, MODE), which are then used to sort them in dedicated astropy Tables, and stored as attributes of a pymusepipe internal class. Each file is then classified and sorted within predefined categories (including standard stars, flats, biases, arc lamp or illumination exposures, astrometry and geometry calibrations, science and sky exposures), the file structure prescriptions and properties being defined via the pymusepipe configuration module. These are further used for all the subsequent data reduction steps, always prioritising calibration files that are close in time except if otherwise specifically requested by the user (again via a configuration file). The local folder structure is initialised according to the outcome of the file scanning, sorted by file types and categories, following a configuration-dependent pymusepipe hierarchy. The data organiser is both reflected, as mentioned, in the Python structure, but also in a set of astropy FITS tables written to disc.

4.2.2. Instrument signature

The MUSE instrument signature is built up and removed using a set of calibration files (biases, flat fields, arc lamps, illumination frames, twilights) acquired soon before or after the main science exposures. These calibrations are processed to derive Master frames including a Master Bias and Flat (via ESO Recipes Execution Tools - EsoRex - recipes, muse_bias, muse_flat), as well as a Trace Table containing the tracing solution for each individual IFU. We note that dark current levels are less than 1 electron per exposure and can be neglected: no dark current correction was applied.

The wavelength calibration solution (EsoRex recipe muse_wavecal) is then derived given a fixed line catalogue. The wavelength solutions are very stable over the 6-year period (October 2014 to December 2020) within which the PHANGS data were accumulated, the distribution of RMS residuals having a median of 0.027 Å (and an average of 0.027 ± 0.004 Å). There is also no detectable difference between the AO and non-AO WFM modes. The next step involved the derivation of the LSF (via muse_lsf) using the same arc lamps used for the wavelength calibration. Static calibrations provided with the raw data were used both for the geometry and astrometry corrections: hence we did not run EsoRex recipes muse_geometry and muse_astrometry. We should emphasise that other joint astrometry/geometry files, calculated via the MUSE-WISE system (implemented as part of the MUSE Guaranteed Time Observations; Vriend 2015), may be more closely following the time varying geometric distortions in the instrument: we will carry out such tests and address potential implications in future data releases (see Sect. 7). The full three-dimensional illumination correction and the relative through-put for each IFU was derived via an illumination correction (muse_twilight) using twilight sky flat frames. For PHANGS data taken before March 11, 2017, vignetting correction was included as recommended in the MUSE User’s Manual⁴. Each offset sky exposure was used via the EsoRex recipe muse_create_sky to produce a sky spectrum, which is then associated with proximate individual exposures minimising the time difference between the exposures on the offset sky and on target. The main steps covered by this data flow are represented on the left-hand side of Fig. 7, and powered via MUSE DRS routines using EsoRex shell piped commands.

Two key stages needed specific attention and proved challenging in the processing of PHANGS data: the sky subtraction for extended sources, and the mosaicking of final data cubes, including astrometric and absolute flux calibrations.

4.2.3. Satellite trails

A few individual exposures included satellite trails producing bright streaks partially or fully crossing the MUSE FoV. These trails were removed from individual exposures by using manually-defined masks. These masks follow a slit-like geometry, with widths between 9 and 12 MUSE spaxels and lengths covering the individual trails detected on reconstructed images before being fed into the mosaicking module. This trail correction step was applied to the first science exposure of pointing 4 for NGC 1365, the third exposure of pointing 3 for NGC 1672, and (a partial trail in) the first exposure of pointing 4 for NGC 3351. Those three trail-masks were fed into the mosaicking module, selecting out the corresponding pixels of the PixTables before proceeding.

4.2.4. Alignment and mosaicking

MUSE DRS provides basic alignment capabilities, for example via cross-correlation techniques and detection of point sources. For extended objects, like targets in the PHANGS sample, the preliminary relative offsets provided by MUSE DRSmuse_exp_align are unfortunately not robust enough to deliver sub-spaxel accuracy, due to the lack of bright point sources. We are also seeking an absolute astrometric solution for each individual exposure that optimises the mosaicking and the comparison with miscellaneous PHANGS datasets (e.g., HST, ALMA). We therefore decided to match the astrometry of each individual exposure with the R-band imaging acquired in the course of the PHANGS project (Sect. 3.2 and Razza et al., in prep.). We used the dedicated alignment module in pymusepipe to connect each individual exposure with the WFI and du Pont images. MUSE images reconstructed from individual PixTables via the WFI (or du Pont) filter curve are compared, given a first guess relative offset, via a set of associated contours and image plots. This step is manual in the sense that the user checks and either confirms or fine-tunes the offset via the Python interface. This process leads to a somewhat subjective assessment of the relative alignment, including a pixel offset (in both X and Y) and a potential rotation angle. To increase robustness, this step is performed by at least two (sometimes three) members of the PHANGS team. Visual inspections confirmed that the relative astrometry is well within one fourth of a spaxel (namely, $0\stackrel{″}{.}05$): any change larger than 0.1 spatial pixel in an individual exposure is easily identified visually and leads to strongly asymmetric residuals (especially around point-like sources) when dividing the reference and MUSE images.

We identified a few issues associated with the geometric and astrometric solutions provided via predefined MUSE calibrations. About 20% of all exposures exhibit a global small but still significant rotation between 0.1 and 0.3 degrees with respect to the R-band images, with no apparent correlation with RA, Dec or time when the target was observed. This residual rotation is also corrected simultaneously via the pymusepipe alignment module, as a free parameter set by the user: such a rotation is also confirmed by at least two different users (sometimes three).

4.2.5. Sky background subtraction

As mentioned in Sect. 4.2.2, we associate each target exposure with the sky spectrum built from the offset sky exposure closest in time. The sky subtraction itself is included in the muse_scipost processing step, thus using an associated continuum spectrum, and always re-fitting the sky emission lines (via the default skymethod=’model’ option). We used the reference R-band images to further constrain any residual (sky) background resulting from an imperfect sky subtraction, or global flux normalisation discrepancy (per exposure). Assuming that the R-band reference image has zero background and the correct absolute flux normalisation, and that the flux in the MUSE reconstructed image represents a linear function of the true flux (involving a normalisation constant plus a background), we can write

$$\begin{array}{cc}& {\mathrm{Flux}}^R(x,y)=a\times {\mathrm{MUSE}}_1^R(x,y)\hfill \\ \hfill & =a\times ({\mathrm{MUSE}}_{\mathrm{raw}}^R(x,y)+\mathrm{Sky}-{\mathrm{Sky}}_1)+b,\hfill \end{array}$$(1)

where Sky is a constant representing the true sky background for that specific exposure, Sky₁ is another constant representing the actual value removed during the initial sky subtraction process, and a and b are constants representing a linear regression representation of the Flux_R−band versus the MUSE reconstructed image. A perfect sky subtraction and normalisation would lead to a = 1 and b = 0. We then use the fitted a value as a normalisation correction, and b to fix the sky contribution by applying Sky = α × Sky₁ where α = 1 − b/(a⋅Sky₁). Hence, knowing a and b as well as Sky₁, the value of the sky continuum integrated within the reference image filter, we derive a correction for the sky normalisation that yields a linear regression where b = 0. The pymusepipe package implements this approach as an option, using the recorded linear regression a and b values. The regression itself is performed via an orthogonal distance regression (ODR) comparing the reference and MUSE reconstructed images after noise filtering and binning: we use bins of 15 × 15 spaxels (3″ × 3″) to minimise the impact of unresolved structure in the comparison.

We find that the distribution of the scaling factors a over the full set of PHANGS-MUSE exposures is well fit by a skew-normal distribution with location 0.99, scale of 0.06 and shape parameter α (related to the skewness) of 1.5 (the best Gaussian fit has a mean of 1.03 and sigma of 0.046). Only 28 (respectively, 27) exposures out of 676 have scaling factors lower than 0.9 (respectively, higher than 1.2). The distribution of background values (b) resembles a Gaussian function centred on 0 with a FWHM of about 1.4 (in units of 10⁻²⁰ erg cm⁻² s⁻¹ Å⁻¹), with a small tail towards positive values. It is important to note that the sky re-normalisation only acts within the R-band filter, assuming that the reference image is background free. Since the reference MUSE sky exposure may result in a reference sky spectrum that is not necessarily an exact representation of the actual sky on the MUSE science exposure, this could lead to a colour variation and, hence, to an over- or under-subtraction of the sky that depends on wavelength (see Sect. 6.1.4).

4.2.6. Point spread function

As emphasised in Fig. 6, our sample galaxies extend well beyond a single MUSE FoV. Covering a significant fraction of the galactic discs (typically well beyond 3 effective radii), each target has been observed using from 3 to 15 MUSE individual pointings, each of them observed at different times and with different sky conditions. The characteristics of the PSF thus naturally vary across a given mosaic. Some of the science goals of the PHANGS project require a homogeneous spatial resolution throughout each galaxy disc. Moreover, a good characterisation of the local PSF at a given location and wavelength within the mosaicked data cube is a key step for the exploitation of such a dataset.

We have assumed that the MUSE PSF can be described at all wavelengths as a circular Moffat function (Fusco et al. 2020), parameterised by its core width and power index. Both quantities are related to the seeing and general atmospheric conditions, combined with the instrument performance during the observations. The core width serves as a proxy for the width of the PSF, while the power index relates to the relative amplitude and shape of the PSF wings. The core size R and the Moffat FWHM are connected via the following relation:

$$\begin{array}{c}\hfill \mathit{FWHM}=2·R\sqrt{2^{1/n}-1},\end{array}$$(2)

where R is the core width and n is the power index of the Moffat function. For a Moffat function, a small power index corresponds to more prominent wings (i.e. as compared with a Gaussian of the same core size). As the power index increases, the Moffat function converges towards a Gaussian profile.

Previous observations with the MUSE WFM have suggested that the MUSE PSF can be considered as a constant over the FoV of the instrument for a single exposure (e.g., Serre et al. 2010; Bacon et al. 2017; Fusco et al. 2020). However, the MUSE spectrographs deliver a PSF whose width varies significantly with wavelength. In order to characterise the PSF of our data cubes, we therefore made use of a four parameter function: a reference core width (or FWHM), a power index at reference wavelength, and a first order polynomial (two parameters) to describe the rate of change with wavelength. These can be directly measured by fitting a Moffat function to a point-like source (i.e. a star) at different wavelengths for each pointing.

In practice, only a fraction of the observed pointings (∼40%) include a star bright enough to robustly recover the four parameters describing the PSF. We thus implemented an alternative method, namely an adaptation of the algorithm presented in Bacon et al. (2017). The principle relies on a minimisation of the difference between a reference image with a known PSF, and a MUSE reconstructed image, after spatial cross-convolution. Such an algorithm requires (i) a reference image with a known (or measured) PSF, which covers a significant fraction of the region of interest (i.e. the FoV of the MUSE pointing) and (ii) an image of the MUSE pointing extracted using the same pass-band as the reference image.

As mentioned in Sect. 3.2, we used the R-band images obtained either with WFI at the 2.2m MPG/ESO telescope in La Silla or with the du Pont telescope at the Las Campanas Observatory (see Razza et al., in prep.) as reference images for the PSF measurement process. The PSFs of the reference R-band images were consistently measured using a set of bright stars throughout their large FoV, and assuming a Moffat two-dimensional profile (which has been shown to be a good representation of their PSF; see Razza et al., in prep.).

Within our PHANGS-MUSE pointings we never cover galaxy-free regions that can serve to directly measure the sky background. Since the galaxy emission is highly structured, providing a robust and independent estimate of the power index of the Moffat using such an automatic algorithm has proven quite challenging. Furthermore, we only have a single reference image per galaxy. This means that we can only measure the PSF of the MUSE reconstructed image as a luminosity-weighted value corresponding to the R bandpass⁵.

We therefore restricted the minimisation process to a single-parameter optimisation, namely the FWHM of the PSF at 6483.58 Å, the reference wavelength of the R-band WFI filter (but see Sect. 7.3). As a first approximation, we followed Bacon et al. (2017) and assumed that the FWHM of the PSF linearly decreases as a function of wavelength, with a slope of −3.0 × 10⁻⁵ arcsec/Å (MUSE team, priv. comm.), while the power index is fixed at the values of n = 2.8 and 2.3 for the no-AO and AO modes, respectively (see e.g., Fusco et al. 2020). While we cannot directly confirm the robustness of such a rule as applied to our dataset, due to the lack of isolated bright stars in the MUSE pointings, preliminary checks with a few sources appear to be consistent with that assumption.

Our specific implementation for the PHANGS-MUSE data cubes consists of the following steps: (i) the reference image is converted to the same flux units as the MUSE image; (ii) the reference image is re-projected to match the MUSE gridding, and the two images are aligned using their respective World Coordinate System (WCS) coordinates (and their cross-normalisation checked); (iii) the MUSE image is convolved with a model of the measured PSF of the reference image; (iv) the reference image is convolved with a given model of the MUSE PSF, with its FWHM as a free parameter; and (v) foreground stars identified using Gaia-DR2 catalogue are masked. The last two steps (convolution of re-gridded narrowband image and masking) are used as input for a least-square optimisation process.

We note that the WFI PSF was estimated by building an effective PSF using a few bright sources in the image already re-gridded to resemble the MUSE pointing (FoV and sampling). Finally, while in Bacon et al. (2017) the stars are masked before the convolution, and the minimisation is performed in Fourier space, we decided to mask the stars after the convolution and to keep the computation in real space. Due to the complexity of the background, masking the stars before the convolution step would create potentially strong edge effects that could drive the minimisation towards incorrect values. Masking after the convolution and deriving the difference in the image space allow us to avoid this issue, at the expense of having a more time-consuming algorithm. Table A.1 summarises the properties of each pointing, including the FWHM of the PSF measured with the algorithm described in this section, and Fig. 8 provides histograms of the measured distribution for the FWHM measured on individual pointings and on completed mosaics for all targets. The minimum, maximum and median values of the measured FWHM are $0\stackrel{″}{.}38$, $1\stackrel{″}{.}18$ and $0\stackrel{″}{.}69$, respectively, with about 80% and 97% of all pointings having FWHM smaller than $0\stackrel{″}{.}8$ and 1″, respectively (only 3 out of 168 having FWHM between 1″ and $1\stackrel{″}{.}2$).⁶

Fig. 8. Distribution of the FWHM of the MUSE PSF. Top panel: FWHM per pointing as measured by the algorithm described in Sect. 4.2.6. The histograms respectively identify the pointings observed in AO mode (orange) and noAO mode (purple). The red dotted histogram represents the distribution for all existing pointings, while the dashed black histogram shows the distribution of PSFs of the copt data cubes (i.e. after the convolution with the worst observed PSF for each galaxy). Bottom panel: average FWHM per galaxy and the convolved FWHM (again, using the worst PSF in each target as the baseline), with the same colour scheme as in panel (a). The arrows indicate the median value of each sub-sample (AO, noAO, all exposures, convolved) with the associated colour scheme.

4.2.7. Line spread function

The LSF is in principle not affected by the observing conditions and thus solely determined by the instrument characteristics. MUSE is a relatively robust instrument and the LSF is stable enough that it does not need to be characterised for every exposure or pointing, as we did for the PSF (Bacon et al. 2017). The LSF is observed to change slightly over the FoV, but usually the variation is small enough (< 0.05 Å; Husser et al. 2016) that, for our purposes, it can be considered constant. However, it does change significantly as a function of wavelength, and a good knowledge of its behaviour is critical, for example, to measure the stellar and gas velocity dispersion (via absorption and emission lines, respectively). The MUSE LSF can be roughly approximated by a Gaussian profile (Bacon et al. 2017) whose FWHM changes as a function of wavelength. The FWHM varies between about 3 Å towards the blue end of the spectrum (4800 Å) and 2.4 Å at ∼7500 Å (Fig. 2 from Bacon et al. 2017). Over the whole range, this variation can be described by a second order polynomial. In this work, we will make use of Eq. (8) of Bacon et al. (2017):

$$\begin{array}{c}\hfill \mathit{FWHM}(\lambda [\mathrm{\AA }])=5.866\times {10}^{-8}{\lambda }^2-9.187\times {10}^{-4}\lambda +6.040.\end{array}$$(3)

This was shown to represent a fair approximation of the variation in the LSF with wavelength: those variations were measured to have a scatter of about 1 to 3%, representing an average of 0.05 Å over the full MUSE spectral range (see e.g., Bacon et al. 2017; Emsellem et al. 2019).

4.2.8. Post processing

At this stage of the data reduction, mosaicked data cubes whose astrometry and background have been calibrated to match those of the reference R-band images were computed. We refer to these mosaics as ‘native’ (for native spatial resolution): the variation in the PSF over the field and as a function of wavelength is not corrected (see Sect. 4.2.6). The native data cubes have the advantage of having the highest spatial resolution possible with the given observations, while the PSF variation may impair robust measurements throughout the FoV or depending on wavelength. In addition to the set of native resolution mosaics, we produce data cubes with homogenised PSFs, labelled as the ‘copt’ (for convolved, optimised) dataset.

The homogenisation procedure first requires a measure of the input MUSE PSF, as described in Sect. 4.2.6. Our target PSF is a circular two-dimensional Gaussian whose FWHM is constant as a function of wavelength and position within each individual mosaic. A Gaussian target PSF was selected to simplify further post-processing, including convolution to coarser spatial resolutions. We make use of a direct convolution scheme with a three-dimensional kernel (‘kernel cubes’) representing the transfer function from the original to the target PSF. Each individual MUSE exposure is addressed independently, while all exposures from a given pointing adopt the same common pointing-dependent PSF. Measuring the PSF at the pointing level leads to a more robust outcome, while the PSF variations between exposures of a given pointing become largely irrelevant due to the linearity of the convolution scheme. We note that we perform the convolution at the individual exposure level to maximise the number of individual cubes that end up being combined during the final mosaicking step, hence optimising the robustness of the rejection of spurious pixels.

The PSF homogenisation proceeds as follows: First, for each individual pointing, a three-dimensional model of the best-fit Moffat PSF is created, with constant power index and with a FWHM varying as a function of wavelength as described in Sect. 4.2.6.

Second, a three-dimensional model of the target PSF with constant FWHM at all wavelengths is created, assuming a FWHM strictly larger than the worst value measured within the mosaic (position, wavelength). More specifically, we use the worst FWHM value and add $0\stackrel{″}{.}2$ in quadrature, which is a reasonable compromise to reach a robust Gaussian profile at all positions.

Third, the convolution kernel cube is created via the Python package pypher (see Sect. 4.1). Fourth, the resulting pypher three-dimensional kernel is fed into a fast-Fourier-transform two-dimensional convolution scheme, for each individual MUSE exposure belonging to the given pointing, addressing each wavelength slice independently, thus reducing the computational time needed for the operation. The spectral variance is also propagated slice by slice via the relation ${\sigma }_{\mathrm{c}\mathrm{o}\mathrm{p}\mathrm{t}}^2={\sigma }^2\otimes {ker}^2$, as derived from the usual principles of propagation of uncertainties. ${\sigma }_{\mathrm{c}\mathrm{o}\mathrm{p}\mathrm{t}}^2$ is the variance of the convolved cube, σ² is the variance of the native cube, ker is the kernel and ⊗ denotes a convolution.

Fifth, the convolved data cubes go through a step of binary erosion⁷ to lower the impact of edge effects during the convolution. We do not keep track of the correlation that this procedure generates between individual spaxels in the convolved cube. Finally, PSF-homogenised mosaics are built from the individually convolved data cubes, via the same combination algorithm used to create the native data cubes.

We should emphasise that the convolution process leads to high levels of spaxel-to-spaxel correlations, which in principle should be reflected in covariance terms: we are ignoring those in subsequent analyses. The variances reported from the convolved cubes should therefore be used with caution when binning the data to larger spaxel sizes.

The resulting copt mosaics are further processed to create additional lower-resolution mosaics, where the spatial resolution is set to either a fixed physical scale (in parsec) or to a fixed on-sky resolution (in arcseconds). All native and convolved data cubes are analysed in the same way, using a dedicated PHANGS-MUSE Data Analysis Pipeline (DAP), as described in the next section (Sect. 5).

Fig. 9. PHANGS-MUSE reconstructed images. RGB images combining the three i, r, and g (red, green, and blue channels, respectively) SDSS-band reconstructed maps of the 19 targets that are full mosaics from the PHANGS-MUSE Large Programme (including the pilot galaxy NGC 0628) – produced as part of the DRP+DAP flow. The scale in arcminutes, which roughly corresponds to the extent of one single MUSE pointing, is shown at the top left of the figure as a labelled black arrow, while the orange vertical bar on the left side of each panel represents the distance-dependent 5 kpc scale for each target.

4.2.9. Data reduction products

The data reduction flow described in the previous sections delivers a set of data products in various forms, including data cubes and images. They include measurements of the PSFs, listed in Table A.1, alignment tables in FITS or ASCII formats including the applied RA and Dec offsets for each individual exposures, as well as the flux normalisation (Sect. 4.2.5). The final products that feed the DAP (see next section) include the full mosaic data cubes for each galaxy target (native, and copt; see Sect. 4.2.8). We also generate a set of reconstructed broadband images for an extensive set of given filters (i_SDSS, r_SDSS, g_SDSS etc; see Fig. 9).

5. The data analysis pipeline

5.1. Requirements and design

The aim of the DAP for PHANGS-MUSE is to generate high-level data products (i.e. fluxes, kinematics, etc), from the stellar continuum and absorption lines, as well as from the ionised gas emission lines. Several software tools have been developed to this aim in recent years, especially to process the data associated with large IFU surveys: CALIFA, MaNGA, SAMI, and TIMER. An incomplete list of notable tools includes Pipe3D (Sánchez et al. 2016a,b), the MaNGA DAP (Belfiore et al. 2019; Westfall et al. 2019), LZIFU (Ho et al. 2016), and gist (Bittner et al. 2019).

We impose several requirements to the software framework for the analysis of the PHANGS-MUSE data, namely: (i) integrate an adaptive spatial binning scheme, to ensure a minimum signal-to-noise ratio required for a robust measurement of observable quantities; (ii) perform a well-tested and robust extraction of physical quantities, including gas and stellar kinematics, and stellar population properties; (iii) parallelise the spectral fitting over multiple cores, to allow for efficient processing of ∼10⁶ spectra (possibly several times, when iterating over subsequent data releases); (iv) be sufficiently modular to allow us to implement changes and/or replace individual modules, without affecting the structure of the code; and (v) support the analysis of IFU datasets from other instruments and surveys, to compare with publicly available high-level products, and therefore benchmark the output of our code against best practices in the field.

We judged these requirements to correspond closely to the philosophy behind the gist code⁸ (Bittner et al. 2019), which has a module-based structure and supports parallelisation for the fitting stages. We therefore adopted gist as the starting point for our pipeline environment.

The modular structure of gist allowed us to easily replace several of its constituent modules with algorithms more closely aligned to our goals with PHANGS-MUSE. In particular, with respect to the public implementation of gist, we made changes to virtually every module, and replaced the emission line spectral fitting and the stellar population analysis routines with ones written by members of our team (F. Belfiore and I. Pessa). The version of the code used for the analysis of PHANGS-MUSE is publicly available⁹.

Our pipeline implementation, which we refer to as DAP is described in detail in the next subsections. Several of these pipeline-level software tools share core pieces of software to perform spatial binning and spectral fitting. Two such modules stand out for their wide applicability: vorbin (Cappellari & Copin 2003) and pPXF (Cappellari & Emsellem 2004). These were originally developed for the pioneering IFU work performed as part of the SAURON/ATLAS^3D surveys (de Zeeuw et al. 2002; Cappellari et al. 2011), and subsequently updated and upgraded (Cappellari 2017). They address the key tasks of binning two-dimensional data to reach a specific S/N level, and to fit the stellar continuum (and, optionally, the gas emission lines) with a non-negative linear combination of templates. Below we briefly described these modules, which we utilised in the analysis of the PHANGS-MUSE data.

vorbin. This is a robust and broadly used package to adaptively bin data cubes along the two spatial dimensions. The method uses Voronoi tessellations (Cappellari & Copin 2003), an optimal solution being found via an iterative process constrained by a given parameter representing the targeted signal to noise. For the present datasets, we used the estimate of the signal-to-noise ratio per spaxel as direct input.

pPXF. The extraction of the stellar and gas kinematics, and information pertaining to the stellar population content of spectra, has been popularised via a set of excellent pieces of software that exploit spectral features in various ways. pPXF (Cappellari & Emsellem 2004; Cappellari 2017) is an intensively tested and robust algorithm to perform direct pixel fitting of spectra making use of spectral template libraries. The generic fitting module of pPXF is extremely flexible and supports a wide range of applications, including the simultaneous fitting of absorption and emission lines and the extraction of non-parametric star-formation histories, with the possibility to add multiple kinematic components and generic constraints (e.g., line flux ratios). The DAP consists of several modules wrapping pPXF for specific applications.

5.2. The data analysis flow

In this section we describe the analysis flow developed within the DAP going from the preparatory to the main computational steps delivering specific data products. The DAP workflow consists of a set of modules running in series, some depending on the outcome of previous steps (e.g., derivation of the stellar kinematics is needed as input to the emission lines fitting module). Figure 10 contains a schematic representation of DAP, including the main input and output parameters, libraries and constraints. Each of the individual DAP modules writes to disc a set of intermediate output files, which can be useful to re-run specific modules in case of a failure, and also contains a more extensive set of outputs that may be of interest for specialised analysis. The key set of physical parameters produced by the DAP are then consolidated into a main output file, described in Sect. 5.2.6.

Fig. 10. Analysis flow for the PHANGS-MUSE DAP. After spectral re-binning using a ln(λ) prescription, the resulting mosaicked cubes are used for the extraction of the stellar kinematics, stellar extinction, and stellar population maps, using an appropriate adaptive binning scheme, while the emission line maps (including the corresponding gas kinematics) are derived from the originally sampled cubes. The main information pertaining to the stellar libraries, spectral masks, assumptions and fixed input for each step are mentioned in the upper part of the figure, while the data products are illustrated in its lower part. See Sect. 5.2 for further details.

5.2.1. Preparatory steps

The DAP requires as input a configuration file, which specifies certain pipeline parameters, the location of the input IFU data, the galaxy’s redshift (or systemic velocity), and the Galactic extinction at its position.

The data are read via a bespoke input module, which is instrument-specific (tailored to MUSE in our case), but may be replaced in order to process data from different instruments. By appropriately changing the parameters in the configuration files and writing a new input module, users can easily process data from other surveys. For PHANGS-MUSE, we utilise a common data input model for both the WFM-noAO and WFM-AO observations. The wavelength range corresponding to the AO gap is automatically masked by MUSE DRS in the case of AO observations. This mask is propagated by the DAP. We assume the LSF given by Eq. (3). The systemic velocity of each galaxy is taken from Lang et al. (2020), who derived these from an analysis of the PHANGS-ALMA CO(2–1) kinematic maps. The MUSE data are corrected for foreground Galactic extinction, using the Cardelli et al. (1989) extinction law and the E(B − V) values from Schlafly & Finkbeiner (2011).

Since we chose to produce fully reduced MUSE data cubes on a linear wavelength axis¹⁰, we perform a resampling of the data on a logarithmic (natural log) wavelength axis, as required for input to pPXF using a channel size of 50 km s⁻¹. This channel size is sufficient to Nyquist sample the LSF of the MUSE data with more than two pixels for λ <  7000 Å, but inevitably over-samples it at the blue edge of the MUSE wavelength range. As discussed further below, we fit the wavelength range 4850−7000 Å to avoid strong sky residuals in the redder part of the MUSE wavelength range.

5.2.2. Spatial binning

Fig. 11. Histograms (violin plots) of the distribution of Voronoi bin sizes for the galaxies in our sample. In orange we show bins within 0.2 R25 and in blue those at larger galactocentric radii. Galaxies are ordered from left to right by stellar mass. The median bin size outside the central regions (R >  0.2 R25) is between 10 and 100 spaxels, roughly corresponding to circularised radii of 0.4″ and 1.1″. The 25th, 50th, and 75th percentiles of each distribution are marked with dashed lines.

As extensively discussed in the literature, an accurate and unbiased determination of the stellar kinematics and stellar population properties from full spectral fitting requires a minimum signal-to-noise ratio (S/N; Johansson et al. 2012; Westfall et al. 2019). The DAP computes the continuum S/N in the 5300−5500 Å wavelength range, using the noise vector from the reduction pipeline.

The data are then Voronoi binned to a target S/N of 35 making use of vorbin. This S/N level is used to determine both the stellar kinematics and the stellar population properties. The value of 35 was chosen to ensure that the relative uncertainty in the stellar mass measurement, which we estimated via Monte Carlo realisations of the errors (see Sect. 5.2.4 below), stays below 15%. For comparison, the MaNGA data are re-binned to an S/N = 10 to determine the stellar kinematics by the publicly available run of the dedicated analysis pipeline (Westfall et al. 2019). We opted for a higher S/N threshold for two reasons: a) we aimed to keep the same Voronoi bins for both the stellar kinematics analysis and the determination of stellar population properties via full spectral fitting, which generally require higher S/N, b) a S/N target of 35 still generated bins that are generally of small size, comparable with the scale of the PSF, except for the edges of the maps where bins can be significantly larger.

In Fig. 11, we show histograms of the number of spaxels included in a Voronoi bin across our sample of galaxies. The histograms are smoothed via a kernel density estimator for ease of visualisation. Bins within 0.2 R₂₅ are shown in orange, while those at larger radii are presented in blue. Galaxies are ordered by stellar mass from left to right. The 25th, 50th, and 75th percentiles of each distribution are marked with dashed lines. The figure confirms that the median bin size outside the central regions (R >  0.2 R₂₅) is between 10 and 100 spaxels, corresponding to circularised radii 0.4″ and 1.1″ (assuming bins are roughly circular). The two lowest-mass galaxies in our sample (NGC 5068 and IC 5332), together with NGC 2835, exhibit lower surface brightness levels even in the central regions, and therefore require larger bin sizes. Other galaxies show a significant number of bins consisting of one, or a few spaxels (which appear as isolated peaks in the figure due to the log-scaling of the y-axis and the smoothing of the distribution). We also note that there exists a tail of bins with a very large number of spaxels (> 1000), but these are very uncommon and correspond to the outermost regions.

Fig. 12. Maps of the average noise (left panel), S/N (middle panel), and binning map (right panel) for NGC 4535. The noise and S/N maps are computed by averaging the pipeline noise and flux over the 5300−5500 Å wavelength range. The stripes dividing the surveyed area into six squared subregions correspond to the overlap regions of the six MUSE pointings obtained for this galaxy. The noise map also shows an evident cross-hatch pattern within individual pointings due to the cube-generating resampling step in the MUSE DRP when combining exposures with different rotation angles. This behaviour is also visible in the S/N map but does not significantly affect the results of the binning process. The binning map shows the result of the Voronoi binning procedure with target S/N = 35. The black contour shows the S/N = 12 level on individual spaxels. In this target only the galaxy centre (and a few foreground stars) have S/N >  35 in individual spaxels, which are therefore left unbinned.

The large size of our MUSE mosaics, well beyond the size of images intended for use with the vorbin package by their authors, meant that we had to limit the size of an optimisation loop in the vorbin code in order to ensure convergence. When calculating the errors for the binned spectra, we assume the MUSE spaxels to be independent. This is not strictly correct, because of the inevitable resampling of the raw data in the data-cube-generation process. In short, spaxels nearly congruent with a single detector pixel will have almost no covariance, while spaxels whose flux originates from several pixels at the detector level will have errors that are strongly correlated with their nearest neighbours (see Sect. 4.6 of Weilbacher et al. 2020). In our native data cubes, this effect is visible as weak horizontal and vertical bands (the two orientations are due to the set of 90deg rotations we perform) in the noise maps. An example map of the average noise vector and the resulting S/N distribution in the 5300−5500 Å wavelength range are shown in Fig. 12. The striping pattern we are referring to should not be confused with the much more evident change in the noise properties of the data in the overlap region of the different MUSE pointings, which corresponds to a real reduction in the noise due to longer effective exposure times. Within mosaics where one pointing consists of fewer than the nominal four exposures (e.g., because one exposure was discarded for quality control reasons) the noise is also consequently higher.

In Fig. 12, we also show the resulting Voronoi binning map, demonstrating that the small-scale noise striping pattern does not have a visible effect on the resulting Voronoi bins. In this work, we therefore neglect the issue of small-scale spatial covariance in the MUSE data.

The DAP supports the determination of emission line properties for two different binning schemes: either the same Voronoi bins as the stellar kinematics or for single spaxels. For the first public PHANGS-MUSE data release (DR1.0, hereafter; see Sect. 7), the emission line properties are derived for single spaxels because the Hα emission is detected at the single-spaxel level across most maps.

We also tested the pipeline with different binning schemes. Two such implementations, optimised for the study of H II regions and the DIG, respectively, are discussed in Santoro et al. (2022) and Belfiore et al. (2022).

5.2.3. Stellar kinematics

The stellar kinematics are derived using pPXF, following the same procedure as in gist. Briefly, to fit the stellar continuum we use E-MILES simple stellar population models (Vazdekis et al. 2016), generated with a Chabrier (2003) initial mass function, BaSTI isochrones (Pietrinferni et al. 2004), eight ages (0.15−14 Gyr, logarithmically spaced in steps of 0.22 dex) and four metallicities ([Z/H] = [ − 1.5, −0.35, 0.06, 0.4]), for a total of 32 templates. We again fit the wavelength range 4850−7000 Å in order to avoid strong sky residuals in the redder part of the MUSE wavelength range. The regions around the expected positions of ionised gas emission lines are masked. The mask width is taken to be ±400 km s⁻¹ of the systemic velocity of the galaxy. This mask width is found to be appropriate for the range of velocities and dispersions present in the PHANGS data. In MUSE-noAO observations the region around the Na I D absorption doublet is also masked, because of the potential ISM contribution. Finally, we mask the region around the bright sky lines at the observed wavelengths of 5577, 6300 and 6363 Å.

The spectral resolution of E-MILES¹¹ is higher than that of the MUSE data within the wavelength range we are considering (see Fig. 13), although at ∼7000 Å E-MILES is expected to have virtually the same spectral resolution as MUSE. The templates are therefore convolved to the spectral resolution of the data, using an appropriate wavelength-dependent kernel. No convolution of the data is performed where the E-MILES resolution is worse than or comparable to that of the MUSE data (beyond about 6750 Å). Velocities are computed with respect to the systemic velocity of the galaxy. We fitted four moments of the line-of-sight velocity distribution (i.e. velocity, velocity dispersion, h3 and h4). To derive the stellar kinematics we make use of additive Legendre polynomials (12th order, in the spectral direction), and no multiplicative polynomials. Polynomials are found to be advantageous in the derivation of stellar kinematics with pPXF (Westfall et al. 2019); the choice between additive and multiplicative is purely dictated by computing efficiency in this step.

Fig. 13. Comparison of the spectral resolution of the MUSE data (Eq. (3), taken from Bacon et al. 2017) and the E-MILES stellar templates. The greyed-out part of the wavelength range (λ >  7000 Å) was excluded from the fit in this first DR1.0 PHANGS-MUSE data release.

The code is able to perform a Markov chain Monte Carlo (MCMC) based error estimate for kinematic parameters (a feature inherited from gist), but this step was not performed for the current data release (DR1.0). The errors reported for the kinematic parameters are therefore formal errors, as given by pPXF. Examples of kinematic maps obtained from our data are shown in Sect. 8.3.

5.2.4. Stellar populations

The stellar population module also employs pPXF in a sequence of steps optimised for use with our MUSE data. We do not use the original gist code for this module, because we moved away from ‘regularisation’, as explained below. The stellar population analysis is performed on the Voronoi-binned data (i.e. the same spectra used for the stellar kinematics determination).

As in the derivation of stellar kinematics, we used E-MILES stellar population models using the same setup as for the stellar kinematics (see Sect. 5.2.3), but we increased the number of templates to 78, with ages = [0.03, 0.05, 0.08, 0.15, 0.25, 0.40, 0.60, 1.0, 1.75, 3.0, 5.0, 8.5, 13.5] Gyr and [Z/H] = [−1.49, −0.96, −0.35, +0.06, +0.26, +0.4]. We fit the same wavelength range as used for the stellar kinematics determination. We retain the previously derived stellar kinematics parameters by fixing them in the pPXF fit, and mask the regions around emission lines, as in Sect. 5.2.3. Additionally, we mask further regions of the spectrum particularly affected by sky line residuals (wavelengths within ranges [5861.1−5911.7], [6528.8−6585.1] and [6820.0−6990.9] Å). The fit is performed in two stages. In the first stage we determine the extinction of the stellar continuum, as parametrised by a Calzetti (2001) extinction curve. In the second iteration the stellar extinction is kept fixed, and multiplicative polynomials of 12th order, with a mean of one, are used to correct for residual inaccuracies in the relative flux calibration. This two-tiered procedure was demonstrated to be necessary to overcome inaccuracies in the sky continuum subtraction, which can cause subtle changes in the spectral shape. To the best of our knowledge, such residuals do not have a large effect on the derived stellar population properties, except for the E(B − V) of the stellar component, as discussed in Sect. 6.3.1. Additive polynomials cannot be used as they both modify the absorption line equivalent widths (EWs) and affect our measurements of stellar mass surface density.

Given the degeneracy that exists in spectral fitting (i.e. between attenuation, metallicity, and age), the use of regularisation (i.e. along metallicity and age axes) when trying to recover star-formation histories has been promoted for certain situations (e.g., Cappellari 2017). Given the average spatial resolution of ∼50 pc, our data resolve individual star-forming regions with significant contributions from very young stars, resulting in strong variations in the spatial distribution of stellar ages across the galactic discs. Extensive tests showed that forcing a fixed level of regularisation on our dataset leads to star-formation histories with strong biases, which themselves vary strongly from region to region and are difficult to control. While we could envision a scheme to provide a more controlled bias, this is beyond the scope of the present PHANGS-MUSE DR1.0 public release. We therefore decided to rely on Monte Carlo simulations to estimate the uncertainty in the recovered stellar population parameters, and use un-regularised fitting. For each spectrum, we performed 20 Monte Carlo iterations. In each iteration, we added to the input spectrum Gaussian noise with a mean of zero and a standard deviation corresponding to the error vector at each wavelength bin. The uncertainties of stellar population parameters were calculated as the standard deviation of their distributions produced by the Monte Carlo realisations. This is meant as a first-order estimate of the true uncertainties. We only ran Monte Carlo realisations during the second step of the fitting procedure, once the stellar extinction has been fixed. Hence, no error was computed for the stellar E(B − V). The output of pPXF is a vector with the weights of the templates that, when linearly combined, best reproduce the observed spectrum. These weights represent the fraction of the total stellar mass born with a given age and metallicity, and they can be used to produce the final maps, including stellar mass surface densities, and both light- or mass-weighted ages and metallicities. The stellar mass surface density maps include both, the contributions from live stars and remnants. For each pixel, the average age and metallicity are computed as

$$\begin{array}{c}\hfill ⟨log\mathrm{age}⟩=\frac{{\mathrm{\Sigma }}_{i}log({\mathrm{age}}_i)w_i}{{\mathrm{\Sigma }}_i{w}_i}\end{array}$$(4)

and

$$\begin{array}{c}\hfill ⟨[\mathrm{Z}/\mathrm{H}]⟩=\frac{{\mathrm{\Sigma }}_i{[\mathrm{Z}/\mathrm{H}]}_i{w}_i}{{\mathrm{\Sigma }}_i{w}_i},\end{array}$$(5)

where age_i and [Z/H]_i correspond to the age and metallicity of each template, and w_i is its corresponding weight in the linear combination. To convert mass-weighted quantities to light-weighted quantities, we use the M/L of each template in the V-band. We compute luminosity-fraction weights as

$$\begin{array}{c}\hfill w_i^{\mathrm{LW}}=\frac{w_i}{{({\mathrm{M}/\mathrm{L}}_V)}_i},\end{array}$$(6)

where ${\mathit{w}}_i^{\mathrm{LW}}$ corresponds to the luminosity-fraction weight of a given template, w_i its mass-fraction weight, and (M/L_V)_i correspond to its M/L in the V-band. We can use these luminosity-fraction weights to calculate light-weighted properties, following Eqs. (4) and (5). It is possible to use the stellar population weights we derive to produce maps of, for example, stellar mass in different age ranges (e.g., young, intermediate, and old stars), or to study the age versus metallicity relation within individual regions. Such analyses are beyond the scope of the present paper.

5.2.5. Emission lines

Emission lines are fitted by performing an independent call to pPXF, where emission lines are treated as additional Gaussian templates, and the stellar continuum is fitted simultaneously. We do not use the module provided with gist, based on the gandalf implementation (Sarzi et al. 2006). This choice mirrors the philosophy of the MaNGA DAP, and is motivated by the greater flexibility of the pPXF implementation and the extensive testing and experience documented in Belfiore et al. (2019) and Westfall et al. (2019). Some of the code we use to interface with pPXF in this fitting stage was adapted directly from the MaNGA DAP, and makes use of the analytical Fourier transform implemented in version > 6 of pPXF (Cappellari 2017).

The fits are performed on individual spaxels, fixing the stellar velocity moments to the values obtained during the stellar kinematics fitting step within the associated Voronoi bin (Sect. 5.2.3). We tested the effect of leaving the stellar kinematics free and find largely identical results for spaxels with large S/N in the continuum. We use the same set of 32 stellar templates as in Sect. 5.2.3.

The kinematic parameters of the emission lines (velocity and velocity dispersion) are tied together within three different groups, as follows: (i) hydrogen Balmer lines: Hα, Hβ; (ii) low-ionisation lines: [O I]λλ6300,64, [N I]λλ5197,5200, [N II]λλ6548,84, [N II]λ5754, and [S II]λλ6717,31; and (iii) high-ionisation lines: [He I]λ5875, [O III]λλ4959,5007, and [S III]λ6312 (see Table 2).

We tie the intrinsic (astrophysical) velocity dispersion within each kinematic group, prior to convolution with the instrumental LSF. The measured velocity dispersion of lines belonging to the same kinematic group is therefore different, but the difference is that required to bring the intrinsic velocity dispersion into agreement, given the assumed LSF. Because the MUSE LSF changes with wavelength, this is an important effect to take into account.

Initial testing showed that the kinematics of the [O III] line is sufficiently different from that of the Balmer lines to require an independent kinematic component. On the other hand, we found that leaving the kinematics of the Hβ line free generates a large number of non-physical Balmer decrements (Hα/Hβ) at low S/N. The definition of a third group of low-ionisation metal lines may support specific science cases focusing on the DIG, where such low-ionisation lines are prevalent with respect to the hydrogen Balmer lines. Alternative tying strategies can be trivially implemented by changing simple keywords in a configuration file of the DAP.

During the emission lines fit, pPXF is run with eight order multiplicative Legendre polynomials, but no additive polynomials (polynomials are only applied to the stellar continuum templates). The use of additive polynomials would be inappropriate in this fitting stage as they modify the EW of the Balmer absorption lines, and therefore potentially introduce non-physical corrections to the hydrogen Balmer line fluxes. The use of a different set of polynomials for this fitting stage with respect to the stellar kinematics fitting stage, in addition to the different S/N of the continuum (going from bins to single spaxels) and the absence of a mask in the wavelength regions around emission lines contribute to creating subtle differences in the best-fit continua generated in the two, complementary fitting stages. The effects of the different polynomials treatment and of masking versus simultaneous fitting of the emission lines are discussed in detail in Belfiore et al. (2019, their Sect. 5.2), who find that the effect of polynomials can be large (∼10% of the emission line flux), especially for the high-order Balmer lines (not present, however, in the MUSE spectral range). They find, moreover, that the effects of masking versus simultaneous fitting of the emission lines are only evident in the regions of Balmer absorption, causing small systematic changes in the fluxes of hydrogen Balmer lines (< 2% for Hα). In light of these findings we prefer to refit the continuum in the emission lines fitting stage.

We note that alternative approaches to this problem are possible. In the Pipe3d CALIFA analysis, for example, the best-fit continuum from the bins used for the stellar kinematics extraction is simply re-scaled according to the median flux in the constituent spaxels, so that the spectrum that is subtracted at the spaxel level is left unchanged (Sánchez et al. 2016b). In the SAMI public data release, a new fit is performed at the spaxel level, as done here, but to limit the impact of degeneracies only those stellar templates with non-zero weights in the parent Voronoi bin are retained in the fit of the individual spaxels (Croom et al. 2012). In general, the effect of restricting the template library in this way are small if the stellar populations are reasonably uniform within the Voronoi bin (Westfall et al. 2019, specifically their Sect. 7.4.1). At the spatial resolution of PHANGS this may not be always the case, since the distribution of young stellar populations is stochastic on small scales. We prefer, therefore, not to restrict the range of templates used, also in light of the fact that we use a small set of templates to start with. Aside from the correction for foreground Milky Way extinction, applied by the DAP as part of its preliminary steps, no correction for dust is applied within the target galaxy.

An overview of the line maps for Hα, [S II] and [O III] is presented for the full PHANGS-MUSE sample in Fig. 14. The emission line science enabled by MUSE is unprecedented, due to the remarkable sensitivity and wide wavelength coverage of the instrument. When applied to nearby galaxies, we have the further advantage of achieving the 50−100 pc physical scales necessary to isolate individual ionising sources from each other and from the surrounding diffuse ionised emission. Our wide coverage across the galaxy discs enables us to connect these small scales to the large kiloparsec-scale relations determined in the literature (see also Sect. 8.1).

Fig. 14. RGB images of the 19 targets in the MUSE sample obtained by combining the emission line maps for [S II], Hα, and [O III] (green, red, and blue channels, respectively). A 5 kpc double black arrow indicates the fixed physical scale for all panels, just below the RGB labels respectively associated with each emission line map.

5.2.6. Output files

The DAP output that we make publicly available consists of two-dimensional maps of parameters and physical properties of interest. These maps are consolidated in the so-called MAPS file, a multi-extension FITS file. The full set of maps currently contains 133 extensions. These MAPS employ the exact same pixel grid and world coordinate system as the mosaic data cubes produced by the reduction pipeline. The full list of extensions of MAPS is described in Table 3. The extensions include two-dimensional maps of stellar kinematics, ionised gas fluxes and kinematics, and stellar-population-related (stellar mass and age) measurements. In the first public PHANGS-MUSE DR1.0 release, all maps pertaining to the stellar and ionised gas fluxes and kinematics are available, so far excluding the stellar populations that need further attention (see Sect. 6.3.2). DAP also produces many intermediate products that are not part of the first public data PHANGS-MUSE DR1.0 release, but could be useful to expert users for specific science cases (e.g., the individual weights of the best-fit stellar templates, which can be used to compute the mass fraction within specified age bins). We welcome requests for these products from interested members of the scientific community. Depending on the feedback our team will receive, future public data releases may include a larger set of products.

We also emphasise that there is a difference in the way the velocity dispersions for gas and stars are reported. The stellar velocity dispersion is reported as the astrophysical dispersion, corrected for instrumental broadening, while the velocity dispersion of the lines is reported as the observed velocity dispersion. To obtain the astrophysical velocity dispersion for the emission lines, the relevant instrumental dispersion at the wavelength of the emission line contained in the SIGMA_CORR extension should be subtracted from the observed value in quadrature. Since the kinematics of the emission lines are tied within the groups defined in Sect. 5.2.5 there is a large amount of redundancy in the velocities and the astrophysical velocity dispersions (i.e. these quantities are not independent). We nonetheless decided to maintain the redundancy in the data products files, to ensure a format robust against changes in the line-tying scheme. Users should be warned, however, that in the current implementation there are only three independent velocity and velocity dispersions for the emission lines, corresponding to the kinematic groups defined above.

5.3. Star masks

Depending on the area mapped and Galactic latitude, we can have up to tens of foreground stars in a single PHANGS-MUSE mosaic. Using the Gaia all-sky catalogue we can localise and mask these contaminants. However, the Gaia catalogue also includes bright stellar clusters, and even some bright galaxy centres, that we do not wish to mask. We therefore follow a series of steps to create our final foreground star masks. These masks will be publicly released together with catalogues extracted from miscellaneous PHANGS datasets, including the MUSE data, in the near future.

We use the Gaia DR2 archive (Gaia Collaboration 2018), searching for all sources within an 8′ aperture centred on each galaxy, allowing therefore for stars up to the edges of our MUSE mosaic FoVs. We then define a circular aperture for each source that attempts to capture all spaxels impacted significantly by the source. The radius of the aperture is defined as a linear function of the Gaiag-band magnitude, empirically calibrated to the MUSE white-light image. We set as the minimum aperture size the worst PSF in the mosaic.

For each source, we determine the maximum and mean of the calcium triplet CaT EWs at heliocentric (i.e. zero) velocity of the spaxels contained within the aperture (i.e. the sum of three bands centred at the features as compared to median flux over ∼i-band). We also determine the maximal stellar velocity difference between the spaxels within the aperture and the median stellar velocity in the unaffected surrounding spaxels representing the typical galaxy velocity. These quantities are then used to classify each source as a Galactic (i.e. star) or extragalactic object (e.g., cluster).

We label each Gaia DR2 object as a star if the velocity difference is large (> 50 km s⁻¹), and the CaT EW is greater than a determined threshold (CaT EW mean >  0.03 Å and maximum CaT >  0.05 Å). These empirically defined thresholds were determined by examining several spectra of Gaia sources and classifying these by eye. These thresholds were found to prevent almost all galaxy centres, bright clusters and H II regions from being automatically masked. Some faint stars (g-mag >  19.5), however, are missed by these criteria.

Other contaminants that may affect the stellar population fits (i.e. diffraction spikes from very bright stars, background galaxies) are not masked. We also note that some faint stars have been found not to be in the Gaia catalogue, meaning that they will be missed by our masking procedure. NGC 1566 has its centre masked due to its type 1 AGN (i.e. exhibiting strong and broad emission lines as well as thermal emission).

We note that the spectral cubes or analysis data product MAPS that we provide as part of the public data release (PHANGS-MUSE DR1.0) are not automatically masked at the positions of the foreground stars we have identified. The users are, on the other hand, left to judge the suitability of these masks depending on their science goals. We caution, however, that most two-dimensional maps of physical parameters (e.g., emission line fluxes, stellar masses) will be biased or incorrect at the position of foreground stars, because of the incorrect spectral fitting (see Sect. 6.2.1).

6. Quality assessment

To validate our final data products of the PHANGS-MUSE sample we derive several statistical measures within the million spectra of our Large Programme, and also cross-compare with existing data on our galaxies. In this section we present quality tests performed both on the final mosaicked cubes, and the high-level data products of the data analysis pipelines.

6.1. Quality of the mosaicked data cubes

6.1.1. Validation of the photometric calibration

To validate the overall photometric calibration of the cubes, we compare synthetic broadband images against existing SDSS images for the nine galaxies that lie within the SDSS legacy survey footprint (Fig. 15). To construct the synthetic r-band images we apply the Sloan r-band filter curve to the MUSE spectral cubes. We use the SDSS Science Archive Server Mosaic tool¹² in order to construct large-scale r-band images that cover the MUSE mosaic footprint. We then sample both the SDSS image and our MUSE image with uniformly spaced 5″ × 5″ apertures. This aperture size was chosen to be large enough to mitigate any differences in seeing or astrometric offsets between the two datasets. Across this sample of galaxies, the median photometric offset ranges from −0.06 to 0.01 mag. NGC 4321 in particular shows a skewed distribution, corresponding to a larger offset for the apertures with lower fluxes, which could potentially be attributed to sky subtraction residuals in either dataset. Overall, the median offset of the galaxies in this sample is −0.01 mag, corresponding to slightly brighter magnitudes in the MUSE image. This is roughly consistent with the SDSS photometric calibration uncertainty (Padmanabhan et al. 2008). Typical scatter within any galaxy is ∼0.04 mag.

Fig. 15. Comparison of r-band magnitudes measured over 5″ × 5″ apertures within MUSE synthetic images (rMUSE) and SDSS images (rSDSS) for the nine galaxies with SDSS imaging available. Each panel is separated into a main one-to-one comparison plot and an additional plot below showing the magnitude difference (ΔrMUSE − SDSS) versus the SDSS value. Insets in each panel show the histograms of the magnitude difference, with the median offset indicated with a dashed line. Solid lines indicate the identity line (or the 0 value for the magnitude difference). Across this sample of galaxies, the median offset ranges from −0.06 to 0.01 mag. Typical scatter within any galaxy is ∼0.04 mag.

6.1.2. Validation of the absolute astrometric calibration

As discussed in Sect. 4.2.4, not every MUSE pointing includes enough bright stars to perform an accurate determination of the astrometry, or to robustly check it a posteriori. This was the motivation to rely on the full extended (stellar and gaseous) emission from the galaxy, and on its comparison with independent ground-based imaging covering much larger FoVs. Still, across our mosaics there are a total of 96 stars that can be overall used to partly validate the absolute and relative astrometric calibrations of the PHANGS-MUSE survey. The only two galaxies that do not contain any foreground stars within the FoV covered by MUSE are NGC 1385 and NGC 7496.

To validate the astrometric solution of the MUSE data, we compared the positions of stars in the MUSE mosaics (as defined in Sect. 5.3) both with their WFI and du Pont broadband locations and with their Gaia DR1 locations. We first used the PHOTUTILS routine IRAFSTARFINDER (Bradley et al. 2020) on our MUSE mosaics to accurately determine the centroid of each object classified as a star in our stellar masks. When comparing the MUSE positions with the broadband positions (measured with the same procedure), we obtain ΔRA = 0.026″ ± 0.047″ and ΔDec = −0.013″ ± 0.044″: such values are observed consistently across the full PHANGS-MUSE sample and are well within the accuracy expected from our alignment routine, representing only from about 1/5th to 1/20th of a MUSE spaxel size. This confirms that our alignment process is robust, yielding a good astrometric calibration using an intermediate imaging dataset. We note that the relative astrometry between the WFI ground-based imaging and Gaia DR1 is consistent with a null offset, and a scatter significantly better than 100 mas.

We more directly compared the MUSE and Gaia DR1 positions for stars in common and actually find a non-zero offset of ΔRA = −0.011″ ± 0.067″ and ΔDec = −0.10″ ± 0.09″. Such an offset is basically zero in Right Ascension, but represents about half a MUSE spaxel in declination, which appears significant. When using the same small subset of stars to compare WFI and Gaia DR1 star positions, we consistently find a similar offset. We should emphasise that the stars in the MUSE FoVs represent much fewer than 1% of the total number of stars used for the astrometric solution of WFI versus Gaia DR1. Hence, if we now consider stars in the broadband FoV that are not covered by MUSE, but still represent more than 99% of all available stars, we measure a shift of less than 1 mas, with scatters of 70 mas and 88 mas in RA and Dec, respectively.

All this suggests that we have a robust absolute and relative astrometric solution at the level of a MUSE sub-spaxel. There is, however, a clear residual offset with respect to Gaia DR1, observed only for stars that overlap with the PHANGS-MUSE galaxy targets. Such a residual shift may be due to systematics of the Gaia solutions associated with the presence of an extended emission background, or could be more generally driven by the difficulty to derive centroids of stars superimposed on bright emission. In any case, addressing this issue pertains to the WFI dataset and a more detailed analysis will be included in Razza et al. (in prep.) when presenting the imaging data.

6.1.3. Validation of the wavelength solution

As noted in Weilbacher et al. (2020), temperature variations between the daytime arc lamp-based wavelength solution and nighttime observations can result in a zero-point shift of the wavelength calibration of up to 1/10th of the MUSE spectral resolution. As a validation of our zero-point calibration, we compare the Hα velocity centroids, tracing the ionised gas kinematics, with the molecular gas kinematics, as traced by PHANGS-ALMA observations of the CO emission (Leroy et al. 2021a). While these two different gas phases are not necessarily co-spatial, they are predominantly constrained to the mid-plane of the disc and can reasonably be expected to globally trace the same kinematics. Here we used the first moment maps from the CO observations, and convert them to a heliocentric reference frame using the optical velocity definition. To minimise uncertainties and ensure we exclude most of the DIG, we require an S/N >  20 on the Hα line. We find the Hα velocities are systematically offset by 0.6 km s⁻¹ towards smaller values compared to the CO, corresponding to 1/100th of a MUSE channel, with ∼10 km s⁻¹ scatter. This systematic agreement provides increased confidence in our overall calibration accuracy, while the variations could reflect differences in the physical impact of feedback processes on different phases of the ISM.

Fig. 16. Median offset and percentiles of the difference between spectra of the same regions from overlapping pointings. The statistics have been computed from a set of 249 regions covering the full PHANGS-MUSE sample, and the resulting percentile vectors have been filtered to make it legible (keeping the sky line residuals visible). Top panel: percentiles of the distribution of the bias level, normalised by the typical (overlap-region-averaged) noise level. Ninety percent of the spectra have a normalised bias that is typically between 30 and 50% of the noise level, while a small fraction show up at levels of 60−120% of the noise, specifically in the blue or red part of the spectrum. We note that although beyond 7000 Å residuals are heavily contaminated by sky line residuals, the pipeline still constructs a roughly correct noise vector. Bottom panel: percentiles (50, 68.3, 95.5%, and 99.7%) of the distribution of differences normalised by the individual spectra noise level, after subtraction of a wavelength-independent median bias offset (see text). The dashed (respectively, green, yellow, and red) lines show values of 1.12 (12% above 1), 2.24, and 3.36, showing that the noise level is slightly underestimated (by about 12%). A trend is visible towards the redder and bluer end of the wavelength coverage.

6.1.4. Sky subtraction

As described in Sect. 4.2.5, the sky continuum spectrum subtracted from individual exposures is computed from dedicated sky exposures. Unfortunately, the overall sky brightness as well as the detailed properties of the sky spectrum generally vary between and during exposures on a timescale of a few minutes. We therefore expect systematic differences (shape, absorption features) between the actual sky continuum in a given science exposure and the sky continuum extracted from the closest sky exposure. To control for any overall shift in the normalisation of the sky brightness, we apply to our MUSE data a normalising factor, deduced from a reference broadband image as described in Sect. 4.2.5. This normalisation, however, cannot address any wavelength-dependent changes in the sky spectrum, for example slope or global shape. Using our MUSE data, we can estimate the impact of these changes by comparing either spectra from individual exposures within the same pointing (assuming sky variations are the main source of variations once any astrometric shifts have been corrected) or spectra from overlapping pointings. Here we favoured the latter approach, as we can use pointings made of several merged exposures, thus optimising the S/N of the evaluated spectra. In this section, we briefly report on such a systematic comparison.

We identified all overlapping regions between different MUSE pointings throughout our sample (as a reminder, such regions are generally about 2″ wide; see Sect. 4.2.5). For each such overlap region we selected a set of n_s spaxels and form all n_s pairs of corresponding spectra (Spec${}_{\mathrm{P}1}^i$, Spec${}_{\mathrm{P}2}^i$, i ∈ [1 − n_s]), where P1 and P2 refer to a generic Pointing 1 and Pointing 2. We discarded pairs of spaxels located at the edges of the pointings using a binary erosion process, thus focusing on spectra within the central part of the overlapping region. We further discard overlap regions containing fewer than 500 spaxels. The process finally generated a set of 249 overlap regions with a minimum number of common spaxels.

For each overlap region, we then derived the median of all the (paired) differences ${\mathrm{S}\mathrm{p}\mathrm{e}\mathrm{c}}_{\mathrm{P}1/\mathrm{P}2}^{\mathrm{m}\mathrm{e}\mathrm{d}}=\mathrm{m}\mathrm{e}\mathrm{d}\mathrm{i}\mathrm{a}\mathrm{n}({\mathrm{S}\mathrm{p}\mathrm{e}\mathrm{c}}_{\mathrm{P}1}^i-{\mathrm{S}\mathrm{p}\mathrm{e}\mathrm{c}}_{\mathrm{P}2}^i)$: this represents the median ‘bias’ offset between P1 and P2 as a function of wavelength. We compare the median bias to the typical noise level in the spectra over a given overlap region, taken as the square root of the average variance of each pair ($⟨{\mathrm{V}\mathrm{a}\mathrm{r}}_{\mathrm{P}1/\mathrm{P}2}⟩=\sum _i({\mathrm{V}\mathrm{a}\mathrm{r}}_{\mathrm{P}1}^i+{\mathrm{V}\mathrm{a}\mathrm{r}}_{\mathrm{P}2}^i)/2$). The result is illustrated in the top panel of Fig. 16 where we see the median bias typically represents about 10% (respectively, 15, 30, 50%) of the noise level for 50% (respectively, 75, 90, 100%) of the spectra. There is a very significant trend towards the blue and red part of the wavelength coverage (with biases up to 150% of the typical noise level). As expected, we also note that sky line residuals have a significant impact on this budget.

We then evaluated the residual difference between pairs of spectra, after removal of that median bias, as compared to the noise level (i.e. the variance of those spectra). This exercise is aimed at testing the reliability of the noise variance delivered by our data flow for individual pointings. This is illustrated in the bottom panel of Fig. 16 where we present some percentiles of its distribution, over the 249 overlap regions. Assuming a standard normal distribution of the noise in the spectra, we would expect the 68.3, 95.5 and 99.7% percentiles of the absolute value of the differences to correspond to 1, 2, and 3 times the standard deviation. The figure shows that we generally underestimate the noise level by about 12% (the three coloured dashed lines in the bottom panel corresponding to values of 1.12, 2.24, and 3.36). We also observe a clear trend with higher values towards the blue and red ends of the spectra, where the discrepancy goes from, for example, 12% to ∼30%. This is consistent with previous estimates; for example, Bacon et al. (2017) quote a re-scaling of the variance by a factor of about 1.3, which they attribute to the impact of the interpolation process and its impact on the noise covariance.

Overall, we conclude that we may slightly under-estimate the noise level by 10–30% when using the derived variances (and ignoring the covariance terms), and that the bias due to improper sky continuum subtraction is present, but negligible for most of the spectra, but can be significant for about 10–20% of them, especially towards the blue end of the MUSE wavelength range.

By construction, there are no offsets in the broadband colour reconstructed images of individual pointings. We confirm that we observe no systematic differences between adjacent pointings using such broadband filters, a good a posteriori check of our implementation. This is, however, not necessarily true for colours. The spectral dependence of the median bias suggests that the shape of the sky continuum spectrum used for the sky background subtraction process may sometimes depart from the true one. We interpret jumps in the stellar extinction maps (see Fig. 23) as a direct consequence of that discrepancy. Fixing such an issue would require a spectrally dependent correction of the reference sky spectrum itself. This may be addressed by using photometric reference points (e.g., HST imaging) in several bands (as opposed to the single Rc-band used here), but it is beyond the scope of the present release (DR1.0).

Fig. 23. Stellar and gas E(B − V) for three example galaxies. The stellar E(B − V)stars is determined from the Voronoi binned data cubes used for stellar population analysis, while the gas E(B − V)gas is computed from the Balmer decrement, derived from the spaxel-level analysis of the ionised gas emission lines (only regions with S/N >  4 in the line emission are considered). The colour bars for the gaseous E(B − V) is stretched by a factor of two to account for the average ratio E(B − V)stars ∼ 0.5 E(B − V)gas. E(B − V)stars traces similar structures as E(B − V)gas, although jumps are evident when comparing the average level of some MUSE pointings to that of their neighbours (e.g., the bottom-right pointing in NGC 4535 and the three left-most pointings in NGC 4254).

Fig. 17. Histograms of the distributions of reduced chi-square, ${\chi }_{\mathrm{d}.\mathrm{o}.\mathrm{f}.}^2$, of the resulting spectral fits for all the spaxels in all the galaxies in our sample. The reduced χ2 is calculated over the wavelength range fitted by the DAP (4850−7000 Å). The red horizontal line and shaded area corresponds to the median (16th and 84th) percentiles of the distribution across all galaxies, corresponding to ${\chi }_{\mathrm{d}.\mathrm{o}.\mathrm{f}.}^2=1.{20}_{-0.18}^{+0.23}$.

6.2. Quality assessment of the spectral fitting and analysis products

In the following subsections we provide a brief assessment of the reliability of our spectral fitting procedure carried out by the DAP. We focus on a few specific questions: the reliability of the error vectors and consequent χ² of the spectral fits, the errors and detectability of emission lines, and the comparisons of our derived stellar masses with those obtained from 3.6 μm imaging. Detailed discussion of the SFR (and extinction corrections) derived from our data and comparison with UV+IR estimates of SFR are presented in Belfiore et al. (in prep.), and in a summarised form in Leroy et al. (2021a), and will therefore not be repeated here.

6.2.1. Overall quality of the spectral fits

In order to validate the overall performance of the DAP we investigate the quality of spectral fits in terms of the reduced χ² (${\chi }_{\mathrm{d}.\mathrm{o}.\mathrm{f}.}^2$). If the errors provided by the pipeline are correct, we expect ${\chi }_{\mathrm{d}.\mathrm{o}.\mathrm{f}.}^2\sim 1$ for good spectral fits. Vice versa, assuming that the spectral fits are correct, we may consider to re-scale the average value of $\sqrt{{\chi }_{\mathrm{d}.\mathrm{o}.\mathrm{f}.}^2}$ as a correction factor to apply to the error vectors to bring them in good agreement with the residuals. In Fig. 17 we show histograms of the reduced ${\chi }_{\mathrm{d}.\mathrm{o}.\mathrm{f}.}^2$ distributions over the full fitted wavelength range fitted by the DAP (4850−7000 Å) for the galaxies in our sample. Galaxies are ordered by stellar mass from left to right. The 25th, 50th, and 75th of each distribution are marked with dashed lines. The red horizontal line corresponds to the median value (16th and 84th percentiles shaded) across all galaxies ${\chi }_{\mathrm{d}.\mathrm{o}.\mathrm{f}.}^2=1.{20}_{-0.18}^{+0.23}$. The median reduced ${\chi }_{\mathrm{d}.\mathrm{o}.\mathrm{f}.}^2$ value for individual galaxies ranges from 1.0 to 1.4. This demonstrates that our spectral fits lead to residuals in good agreement with the error vectors provided by the pipeline. In all galaxies, however, a very small number of spaxels show much larger ${\chi }_{\mathrm{d}.\mathrm{o}.\mathrm{f}.}^2$ (only 0.1% of spaxels have ${\chi }_{\mathrm{d}.\mathrm{o}.\mathrm{f}.}^2>3$). We can gain more insight into these extreme regions by locating them spatially within our galaxies.

Fig. 18. Maps of the ${\chi }_{\mathrm{d}.\mathrm{o}.\mathrm{f}.}^2$ for three example galaxies at the individual spaxel level. We show in blue contours the Hα surface brightness and in red the regions masked as foreground stars. The ${\chi }_{\mathrm{d}.\mathrm{o}.\mathrm{f}.}^2$ maps demonstrate the quality of our spectral fits ${\chi }_{\mathrm{d}.\mathrm{o}.\mathrm{f}.}^2\sim 1$, except in a few regions of bright Hα emission (e.g., the central regions in NGC 1672 and NGC 4535). A few regions of higher ${\chi }_{\mathrm{d}.\mathrm{o}.\mathrm{f}.}^2$ can also be seen in the middle of spiral arms, corresponding to the centres of bright H II regions. Foreground stars appear as small spots of very high ${\chi }_{\mathrm{d}.\mathrm{o}.\mathrm{f}.}^2$ and are emphasised with red contours.

Fig. 19. Histogram of the distributions of ${\chi }_{\mathrm{d}.\mathrm{o}.\mathrm{f}.}^2$ (in blue) for the same example galaxies as in Fig. 18. The median ${\chi }_{\mathrm{d}.\mathrm{o}.\mathrm{f}.}^2$ is shown as a blue vertical dashed line. We show as a bar plot in orange and green the stacked histograms of pixels belonging to our foreground star masks (orange) and to the brightest line-emitting regions (> 99th percentile of the Hα emission, in green). Because the bar plots are stacked, their sum must be smaller than the histogram of all spaxels in blue. This allows us to visually assess that the vast majority of high-${\chi }_{\mathrm{d}.\mathrm{o}.\mathrm{f}.}^2$ pixels are due to foreground stars and bright line-emitting regions.

In Fig. 18 we show example reduced ${\chi }_{\mathrm{d}.\mathrm{o}.\mathrm{f}.}^2$ maps for three galaxies in our sample. As expected, the median reduced ${\chi }_{\mathrm{d}.\mathrm{o}.\mathrm{f}.}^2$ is very close to 1 (1.2, 1.1, and 1.2 for NGC 1672, NGC 4535, and NGC 4254, respectively). We also see the expected pointing to pointing variations (due to the slight differences in the noise levels across pointings) and the cross-hatch pattern, which has already been discussed in Sect. 5.2.2 (see Fig. 12). It is also evident that a few localised regions have much higher reduced ${\chi }_{\mathrm{d}.\mathrm{o}.\mathrm{f}.}^2$ than the median. One such class of deviant regions is represented by foreground stars, which are enclosed in red contours in Fig. 18 (as defined by our foreground star masks; see Sect. 5.3). In addition, regions of very bright line emission, both in the middle of spiral arms (traced by the Hα emission, shown with blue contours) and the galaxy centres (in NGC 1672 and NGC 4535 in particular) show enhanced reduced ${\chi }_{\mathrm{d}.\mathrm{o}.\mathrm{f}.}^2$. We interpret this as a consequence of non-Gaussianity in the bright emission lines (see Sect. 6.2.2), causing strong residuals at their specific wavelengths.

Finally, we observe (especially in NGC 4254) a trend of increasing ${\chi }_{\mathrm{d}.\mathrm{o}.\mathrm{f}.}^2$ with galactocentric distance. This effect is likely the consequence of the presence of sky residuals in the spectra. These residuals contribute to the ${\chi }_{\mathrm{d}.\mathrm{o}.\mathrm{f}.}^2$ in regions of low surface brightness, while in regions of higher surface brightness the Poissonian errors from the continuum become dominant. Therefore, at high surface brightness levels pipeline errors are a better representation of the model residuals, resulting in reduced ${\chi }_{\mathrm{d}.\mathrm{o}.\mathrm{f}.}^2$ closer to unity.

We investigate further whether foreground stars and bright emission lines can explain all the high ${\chi }_{\mathrm{d}.\mathrm{o}.\mathrm{f}.}^2$ spaxels by presenting in Fig. 19 a histogram of the ${\chi }_{\mathrm{d}.\mathrm{o}.\mathrm{f}.}^2$ distributions for our three example galaxies (blue histograms), computing such a quantity at the individual spaxel level. We plot the $\log ({\chi }_{\mathrm{d}.\mathrm{o}.\mathrm{f}.}^2)$ on the x-axis and employ a logarithmic y-axis to emphasise the small tail of spaxels at high ${\chi }_{\mathrm{d}.\mathrm{o}.\mathrm{f}.}^2$ values. We also present as stacked barplots the distributions of pixels associated with foreground stars (orange) and bright line emission (parametrised here as the 1% of spaxels with brightest Hα emission, green). Since the orange and green barplots are stacked on each other, their total height would reach the blue histogram if all the high-${\chi }_{\mathrm{d}.\mathrm{o}.\mathrm{f}.}^2$ spaxels fall into these two categories. Overall, we see that for NGC 1672 and NGC 4535 we explain virtually all the spaxels ${\chi }_{\mathrm{d}.\mathrm{o}.\mathrm{f}.}^2>3$ with either foreground stars or bright line emission. The absence of foreground stars in NGC 4254 creates a ${\chi }_{\mathrm{d}.\mathrm{o}.\mathrm{f}.}^2$ with a much more limited tail, which is overall largely explained by bright line emission.

We conclude noting that the residuals from our spectral fits are in excellent agreement with the pipeline noise. If we take the median values of ${\chi }_{\mathrm{d}.\mathrm{o}.\mathrm{f}.}^2=1.20$ and assume that our fits are good, then we would conclude that the noise vectors provided by the pipeline ought to be corrected upwards by 10% on average, in excellent agreement with the 12% estimates obtained from the overlap region statistics in Sect. 6.1.4. Such a level of agreement can be considered extremely satisfactory, and highlights the robustness of the error vectors provided the MUSE DRS. Spaxels with highly deviant ${\chi }_{\mathrm{d}.\mathrm{o}.\mathrm{f}.}^2$ values are explained as either foreground stars, which can be masked efficiently by the masks we provide, or regions of bright line emission, where the emission line residuals dominated the overall spectral ${\chi }_{\mathrm{d}.\mathrm{o}.\mathrm{f}.}^2$.

6.2.2. Reliability of errors for emission line fluxes

We obtain emission line errors directly from the output of pPXF, which obtains them as the output of its non-linear fitting stage, performed via Levenberg-Marquardt minimisation, using a Python version of the well-established mpfit routine, originally written in IDL by Markwardt (2009). Several authors have found that when fitting emission lines with single Gaussian the formal errors obtained from this procedure are accurate to within a few percent (Ho et al. 2016; Belfiore et al. 2019). We performed our own simple recovery simulations, where noise is added to a model spectrum consisting of Gaussian emission lines, to assess the validity of the errors under such idealised conditions, and find, in agreement with previous work, that the flux is recovered with no bias and flux errors are correct to within better than 5% percent at all flux levels.

Users interested in errors for weak lines may wish to consider the underestimation of the noise vector of ∼12%, discussed in Sects. 6.1.4 and 6.2.1. Weak lines that are found in regions of enhanced sky residuals should also be considered carefully. A subtler effect affects our estimates of the errors for very strong emission lines. It is found empirically that for S/N >  30, residuals from the Gaussian fit are generally larger than the errors, leading to an increase in χ² of the fit within the core of the line as a function of increasing S/N (see Fig. 3 in Belfiore et al. 2019).

6.2.3. Emission line detection thresholds

We demonstrate the sensitivity of our data by plotting the fraction of pixels above a given surface brightness threshold, for a representative sample of emission lines mapped within 0.5 R₂₅ and at S/N >  3 (Fig. 20). This provides a more quantitative illustration than multi-line maps (Fig. 14), which we hope may be of use to the reader in proposal and project planning. Our brightest line, Hα, is detected in ∼95% of all pixels for most galaxies. Typical 3σ flux sensitivity in Hα is 4−7 × 10³⁷ erg s⁻¹ kpc⁻² (3−7 × 10⁻¹⁹ erg s⁻¹ cm⁻² per pixel). Hβ is detected in 50−80% of pixels, which has implications for our ability to perform matched-resolution extinction corrections based on the Balmer decrement (see also Sect. 8.2). Low-ionisation lines ([N II]6584, [S II]6717, and [S II]6731) are detected in 60−95% of pixels, while the high-ionisation [O III]5007 line emission is less common (50−70% of pixels). The faint auroral [N II]5754 emission line is detected at a 3σ level in ∼5% of pixels (though fewer than 1% are detected at > 5σ). Even moderate spatial binning can help significantly with the detection of low-surface-brightness emission and fainter auroral lines (e.g., Santoro et al. 2022; Belfiore et al. 2022).

Fig. 21. Two-dimensional histogram of the ratio between the stellar mass surface density derived from full spectral fitting from MUSE (Sect. 5.2.4) and archival stellar mass maps from S4G, based on Spitzer IRAC maps (Querejeta et al. 2015), as a function of the light-weighted age of the underlying stellar population in each pixel, derived from the full spectral fitting. We show three galaxies in our sample, spanning the stellar mass range of the PHANGS-MUSE sample (in ascending order from left to right). The figure shows good agreement between the MUSE and the S4G mass maps in pixels dominated by stellar populations older than ∼4 Gyr. However, in regions hosting younger populations, the masses derived from the near-IR data are systematically larger, possibly due to the contribution of AGB stars to the near-IR flux. The contours mark the iso-probability curves where the probability density drops below 40% and 10% of the maximum. The dashed red line shows the median mass ratio at a given luminosity-weighted age.

6.2.4. Stellar masses

We now compare the values of stellar mass surface density derived from our stellar population analysis (Sect. 5.2.4) with those derived by Querejeta et al. (2015) using the 3.6 μm and 4.5 μm IRAC bands from the Spitzer Survey of Stellar structure in Galaxies (S⁴G; Sheth et al. 2010). The S⁴G work uses the independent component analysis (ICA) presented in Meidt et al. (2014) to separate the contributions of the stellar and the dust emission to the IRAC fluxes and to derive a relation that allows the M/L to be obtained using the [3.6]−[4.5] colour, assuming that the contribution to the stellar light at those wavelengths is dominated by an old population with an almost constant M/L ratio.

We convolved the MUSE maps to an angular resolution of 1.5″, consistent with that from the S⁴G maps, and computed, for each pixel, the ratio M_MUSE/M_S⁴G, where M_MUSE and M_S⁴G correspond to the stellar mass derived from our MUSE data and S⁴G data, respectively. On average, we found that M_MUSE values are about 30% smaller than M_S⁴G, but the differences are strongly dependent on the age of the underlying stellar population estimated from the full spectral fitting.

This can be seen in Fig. 21, where we show galaxies spanning the stellar mass range probed by PHANGS-MUSE, NGC 5068 being the least massive object in our sample, NGC 1300 close to the median mass, and NGC 1365 being the most massive galaxy. We show log(M_MUSE/M_S⁴G) as a function of the light-weighted age of each pixel. As can be seen, both methods agree when the light is dominated by old stars, but start to diverge when the mean light-weighted age is lower than ∼4 Gyr. A similar trend was reported by de Amorim et al. (2017) using CALIFA data. This is not surprising, as the S⁴G calibration assumes that the stellar flux in the near-IR is dominated by old stars. However, in stellar populations with ages of 1−2 Gyr, the contribution of Asymptotic Giant Branch (AGB) stars to the near-IR flux can be dominant (leading to a different M/L than that appropriate for and old population) and all the pixels with a mean light-weighted age below ∼4 Gyr have a contribution of stars in this age range. It is worth mentioning that the age dependence of the mass difference is much stronger when a constant M/L ratio is applied to the original 3.6 μm image (i.e. the ICA has some effect in compensating for the excess of light in regions that host young stellar populations) but does not completely remove the age trend.

Figure 22 shows the median M/L_3.6 μm ratio of a pixel, at a given light-weighted age, calculated as the ratio M_MUSE/L_3.6 μm, where L_3.6 μm corresponds to the luminosity of a pixel in the original IRAC 3.6 μm image (not the ICA version from S⁴G), in solar units. Each line represents a galaxy, coloured by its total stellar mass. The typical dispersion in the M/L ratio measurement at a given age across galaxies is ∼0.03. The figure shows little variations of this trend among galaxies, with no obvious correlation with total stellar mass. The horizontal lines in Fig. 22 correspond to M/L values commonly adopted in the literature. The value of Meidt et al. (2014) was intended to be applied to the dust-corrected images (i.e. after applying the ICA to remove dust emission), and is therefore higher than the other values from the literature and closer to the value expected for a very old stellar population.

Fig. 22. Age dependence of the MMUSE/L3.6 μm ratio (in M⊙/L⊙) for the galaxies in our sample, where L3.6 μm corresponds to the luminosity of a pixel in the original IRAC 3.6 μm image (not the ICA version from S4G). Each line represents a galaxy, coloured by its total stellar mass. The figure shows a positive trend, with pixels hosting older stellar populations having larger M/L ratios. The horizontal lines mark different values adopted in the literature (dot-dashed line M/L = 0.41, Schombert et al. 2019; dotted line, M/L = 0.53, Eskew et al. 2012; dashed line M/L = 0.6, Meidt et al. 2014).

Leroy et al. (2021a) compare the MUSE mass maps presented here with their GALEX+ WISE version from z0MGS (following Leroy et al. 2019). As an alternative to the ICA procedure, they use a M/L that scales with GALEX + WISE colours to compensate for the impact of specific SFR on the 3.6 μm or 3.4 μm flux, in the range of 0.2 <  M/L <  0.5. Our mass maps agree well with this range of values. However, they also find an offset of 0.08 dex between their masses inferred from GALEX + WISE imaging and our MUSE-derived mass maps. The offset is found to be roughly independent of specific SFR, which can be understood as a proxy for stellar age, and is therefore likely associated with other systematic differences in the stellar population modelling (e.g., differences in the Single Stellar Populations – SSP – models).

Fig. 24. Age (top row), metallicity (middle row), and extinction (bottom row) of the inner star-forming ring of NGC 3351 obtained using three different sets of templates: E-MILES (left column), CB07, excluding templates younger than 30 Myr (CB07-A; central column), and CB07 including young templates (CB09-B; right column). Adding younger templates to our template grid slightly alleviates the problem of extremely metal-poor and young regions in the sense that the young clusters are younger and mildly less metal-poor. However, the improvement is marginal, and these regions are still significantly more metal-poor than the surrounding pixels. Panels in the last row show, additionally, an abnormally low extinction in the young and metal-poor regions (marked as white circles in the left column). This is likely a consequence of the same issue, as it is partially improved when younger templates are included. Overall, adding templates younger than 30 Myr to our age-metallicity grid does not provide a solution to the issue of young and extremely metal-poor regions.

6.3. Known systematic errors in the stellar population maps

6.3.1. Imperfect sky subtraction: Effect on stellar extinction

As discussed earlier, the stellar E(B − V) map shows clear systematic jumps between different MUSE pointings. These jumps are caused by spectral differences in the sky continuum subtraction across adjacent pointings in the mosaic (see Sect. 4.2.5). We use the overlap regions among adjacent MUSE pointings to quantify the impact of the different continuum levels on the recovered stellar extinction. We computed the stellar E(B − V) for this set of 2 × 249 spectra (two per overlap region), extracted from the individual overlap regions. For each pair of spectra, we calculated ΔE(B − V), as the absolute difference in the extinction between the two. We therefore measured 249 values of ΔE(B − V) across our full sample obtaining a median value (as well as the 68, 95 and 99.7th percentiles, respectively) of ΔE(B − V) = 0.037 (0.058, 0.133 and 0.261, respectively). This test implies than, although most of the pointings are relatively homogeneous with respect to their neighbours, with ΔE(B − V) <  0.04, a fraction of them show larger spectral differences that lead to systematic offsets in the derived stellar extinction, up to ΔE(B − V) ≈ 0.3, in agreement with the pointing-to-pointing jumps visible in Fig. 23. Fixing such an issue would require the usage of different sky continuum reference spectra for individual exposures, something that is envisioned for future public data releases (see Sect. 7.3).

6.3.2. Systematic errors in the stellar population fits at young ages

While examining the maps associated with the stellar population fitting (see Sect. 5.2.4), we noticed the presence of low metallicity values (LW [Z/H] <  −1.3) in a few regions encompassing very young stellar clusters (LW age <  400 Myr). Such low metallicity values would be inconsistent with an internal and progressive chemical enrichment of the ISM (e.g., Ho et al. 2017). This suggests that the fitting process converges towards a misleading local minimum, the bluest available stellar template, constrained by the youngest age bin (30 Myr) of the implemented template library. These low-metallicity regions usually coincide with strong Hα emission, indicating that these young clusters coincide with active star formation. In addition to the lack of younger templates (due to the low number of young and metal-poor stars in the E-MILES library), contributions from nebular emission as well as unmasked emission lines are expected to further impact the χ² minimisation in such ‘young’ regions. A visual inspection of the spectral fits for some associated spaxels revealed a systematic overestimation of the stellar continuum at wavelengths bluer than ∼5100 Å. We therefore deem these age and metallicity measurements (as well as E(B − V)) unreliable. This issue has already been reported in several studies (e.g., Carrillo et al. 2020; Bittner et al. 2020), who similarly reported unexpected young and metal-poor regions.

We explored whether adding younger templates to our age-metallicity grid is sufficient or not to overcome this issue. To this end, we used SSP models from the Bruzual & Charlot evolutionary population synthesis database (Charlot & Bruzual 2007, priv. comm.; CB07), computed assuming a Padova 1994 isochrone and a Chabrier (2003) IMF. In order to account for differences driven by using a different stellar library, we defined two sets of CB07 templates: CB07-A, with five metallicity bins ([Z/H] = [−1.7, −0.7, −0.4, 0, 0.4]) and 16 age bins, log-spaced from 30 Myr to 20 Gyr (i.e. with the same low-age limit as E-MILES), and CB07-B, with the same metallicity bins, but with a larger age coverage, including 25 age bins ranging from 1 Myr to 20 Gyr . Figure 24 presents the output LW age, LW [Z/H] and stellar E(B − V), obtained with each different set of templates (E-MILES, CB07-A, and CB07-B) for the nuclear star-forming ring of NGC 3351.As expected, adding younger templates has an impact, in the sense that the young clusters are younger and mildly less metal-poor. However, the improvement is marginal, the derived metallicities for the above-mentioned regions being still significantly lower than for their surrounding spaxels. The associated low measured E(B − V) is also likely driven by the same mechanism, biasing the blue end of the best-fit spectrum. We conclude here that including templates younger than 30 Myr unfortunately does not solve the degeneracy. A robust solution is beyond the reach of the present paper, and would require a broader consideration, examining a combined set of factors such as nebular emission (continuum and lines), template mismatch due to the currently poor observational constraints in the young and metal-poor regime, and potential sky-subtraction residuals.

7. Data versions

7.1. The first PHANGS-MUSE public data release (DR1.0)

PHANGS-MUSE is a legacy dataset, with wide applications beyond those pursued by the PHANGS team. To this end, we provide public data releases of science-ready pipeline products, including both the fully reduced mosaics (at native, copt), and the DAP high-level data products (at the same two resolution levels). The science-ready data cubes and data products are provided through two main links, via the ESO Archive tool¹³ and also via the CADC portal, where the PHANGS-ALMA dataset is also available.

The first public PHANGS-MUSE data release (DR1.0) includes the datasets and data products as described in Sects. 4.2.9 and 5.2.6. Maps associated with the stellar population analysis will be made available in subsequent releases, as well as other convolved datasets (e.g., 15″). The products associated with DR1.0 were available as an internal PHANGS data release (with minor modifications), and have been already exploited in the course of several studies (published or submitted at the time of this writing), including Leroy et al. (2021a), Pessa et al. (2021), Turner et al. (2021), Williams et al. (2021), Williams et al. (2022), Barnes et al. (2021), Bešlić et al. (2021), Belfiore et al. (2022). An H II region catalogue based on the PHANGS-MUSE data is presented Santoro et al. (2022).

7.2. Early versions of the PHANGS-MUSE data

Results from the PHANGS-MUSE survey have already been published using preliminary versions of the DRP, coupled with different analysis methods. We briefly summarise here the differences between these early versions of the reduced PHANGS data and that presented in the first public PHANGS-MUSE data release (Sect. 7.1).

An initial sample of eight galaxies (IC 5332, NGC 0628, NGC 1087, NGC 1672, NGC 2835, NGC 3627, NGC 4254, NGC 4535) was reduced manually using esorex pipeline recipes (i.e. without the pymusepipe software framework). Stellar continuum fitting and emission line subtraction was performed using LZIFU (Ho et al. 2016). Further details on the data reduction and data analysis are provided in Kreckel et al. (2019). This data version was used in associated early PHANGS papers (Kreckel et al. 2017, 2018, 2020; Ho et al. 2019; Herrera et al. 2020). Both the astrometric accuracy and the absolute photometric calibration were not yet fully developed for this data version. Absolute flux calibration was determined only for NGC 0628, and we find an overall agreement in the Hα line flux between the older maps and the current DAP products to within 2% for this galaxy. Data from the remaining galaxies was used only to determine line ratios, and pixel-to-pixel comparisons between those maps and our latest version of the DRP and DAP find they underestimate line fluxes by up to 10% at S/N >  5 in the Hα line flux, due to non-photometric observing conditions. Line ratios show significantly better systematic agreement, to within 3% at S/N >  5 in [N II]/Hα. H II region catalogues derived from both versions have astrometric agreement of 0.3″ when cross-matching the unresolved objects.

Finally, the DRP and DAP described in this paper were used, in a preliminary form, for the development of tools internally to the PHANGS collaboration. The associated preliminary dataset for NGC 0628 was specifically used by Andrews et al. (2021). The main limitation of these early data versions is in the astrometric calibration, and it did not significantly impact the results in that paper.

7.3. Future developments

This paper describes the data reduction and data analysis steps employed as a baseline in producing science-ready data products. However, as discussed throughout the paper, we are aware of several limitations of the current analysis, which we plan to address accordingly. Our planned next steps are as follows.

Improved sky subtraction. This both concerns the sky continuum and the sky emission lines. The former will be better constrained via a refined estimate of the sky continuum contribution using overlapping exposures and external constraints from multi-band HST images when available. For the latter, we would consider a principle component analysis (PCA) approach, similar to what has been developed by other groups (e.g., Soto et al. 2016). We hope this would result in significantly reduced stellar E(B − V) residuals between pointings (see Fig. 25) within the mosaic, and further scientific exploitation of the red end of the MUSE spectral coverage.

Fig. 25. Distribution of pixel-based Hα/Hβ line ratios as a function of observed (top) and extinction corrected (bottom) Hα surface brightness (SBHα) for three galaxies of our sample (see text). We employ a S/N cut of 10 for both lines. Most pixels have Hα/Hβ >  2.86 (dashed line, as predicted by case B recombination for Te = 104 K, ne = 100 cm−3), and those that have ratios below this value are consistent within the errors. Because of this, when de-reddening Hα we apply no correction for those pixels with Hα/Hβ < 2.86. At low Hα surface brightness, our Hβ sensitivity sets an upper threshold on the recoverable value of Hα/Hβ. The value of this threshold for the median Hβ error is indicated in all panels by the dotted lines. Histograms of SBHα are included for each set of plots for pixels that meet our S/N cut (in black) and those that do not (in grey). A significant number of pixels are excluded by requiring an S/N >  10 detection of Hβ. A wide range of extinctions is observed, with no obvious correlation with SBHα. This demonstrates the need for matched resolution constraints on the dust extinction. Within PHANGS-MUSE, the use of the Balmer decrement provides a unique opportunity to perform robust, extinction-corrected SFR measurements across the bulk of the galaxy disc.

Improved geometric transformation. The data release described in the present paper relies on fixed geometric and astrometric MUSE calibration files, as delivered with the raw datasets. We would examine whether a time-varying set of calibration files (provided via e.g., the MuseWise system implemented by the MUSE GTO team; Vriend 2015) could improve the registration of individual MUSE exposures of the same FoV.

Improved absolute flux zero-point calibration and astrometric registration. We would make use of an improved photometry and astrometry for the ground-based du Pont and Direct CCD imaging, and further exploit the existing PHANGS-HST imaging when possible.

Use of the auto-calibration functionalities in theMUSE DRS to improve on the residual imprints left by the slicer-to-slicer variations (flat-fielding). While this is clearly a secondary and low-amplitude effect, the exploitation of all blank sky exposures during the night associated with individual exposures seems to give reasonable results and a significant improvement on test exposures. It still requires pre-calculated ‘AUTOCAL_FACTORS’ with carefully defined sky masks and, hence, goes beyond the present data release.

Application of alternative stellar population templates. This would include a systematic study of the choice of alternative stellar population templates and its impact on inferred properties. In particular, we would assess whether the identification of anomalously low-metallicity young stellar population associated with high Hα EWs is physical or reflects limitations in the templates currently employed.

Multi-component Gaussian emission line fits. Secondary components or non-Gaussian profiles are sometimes visible for the emission lines in the MUSE data. Those may be associated with, for example, small-scale outflows, complex dynamical regimes (e.g., shocks) or (spectrally, spatially) unresolved structures. A systematic identification of spectra significantly impacted by our default assumption of single Gaussian profiles for each individual emission line is one of the listed improvements to be added to a future data release (see e.g., Henshaw et al. 2020a).

8. Key science enabled by PHANGS-MUSE

The PHANGS-MUSE survey enables a number of key science goals (see Sect. 2.3). In the following subsections we demonstrate the richness of the PHANGS-MUSE dataset.

8.1. Mapping of star-formation rate

Dust proves an obstacle to inferring the direct physical conditions, as it acts to both extinguish and redden the stellar and nebular light. With PHANGS-MUSE, we can directly parameterise the effect of dust on nebular emission lines through the Balmer decrement (traced by the ratio of Hα/Hβ). Figure 25 demonstrates the range of Hα/Hβ line ratios that we observe at S/N >  10 for three galaxies in our sample, as a function of Hα surface brightness. These three galaxies represent the lowest stellar mass (NGC 5068), an average stellar mass flocculent morphology (NGC 4254) and the highest stellar mass, AGN-dominated, strongly barred morphology (NGC 1365). The majority of pixels are above the canonical Hα/Hβ of 2.86 (applicable for Case B recombination, assuming T_e = 10⁴ K and n_e = 100 cm⁻³), revealing the impact of dust reddening across these spiral galaxies. Those pixels with lower values are predominantly consistent with 2.86 within uncertainties, even accounting for our fairly high S/N >  10 cut.

As also reflected in Fig. 20, a large fraction of pixels have significant detections in Hα, but are not detected in Hβ (Fig. 25, histogram). For pixels where both emission lines are detected, we can apply an extinction correction at matched resolution. Applying an O’Donnell (1994) extinction law, and applying no correction for Hα/Hβ < 2.86, the observed Balmer decrements correspond to extinctions up to A_V ≈ 4 mag in these galaxies, which is the typical maximum detectable extinction value for our sample. Most pixels are consistent with A_V ≈ 1−2 mag, typical for many face-on galaxies (Kennicutt 1998). In NGC 1365, where the central AGN contributes to the high Hα intensity pixels, the assumed intrinsic Hα/Hβ = 2.86 is probably not appropriate, resulting in an overestimation of A_V for these pixels. It is worth emphasising that applying a constant intrinsic Balmer decrement might be restrictive as the Balmer decrement is sensitive to temperature: other approaches and methodologies exist to estimate the extinction correction via the emission lines, that is, allowing for a varying Hα/Hβ ratio via the fitting of photoionisation models (see e.g., Mingozzi et al. 2020; Pellegrini et al. 2020).

Applying an extinction correction imposes a correlation between Hα/Hβ and extinction-corrected Hα surface brightness. Due to our Hβ detection limit, we cannot determine an extinction correction for highly extincted pixels at low Hα surface brightness (represented by the dotted line). The small number of points in excess of this upper limit reflect slight variations in the Hβ error map. Accounting for that, we see no clear correlation between observed Hα brightness and dust extinction. The uniformity in Hα pixel-based surface brightness distribution between these three quite different galaxies reflects the uniformity in the H II region luminosity functions, tracing the underlying photoionisation physics powering the majority of the brightest emission (Santoro et al. 2022). The PHANGS-MUSE dataset thus builds on previous work (e.g., Kreckel et al. 2018) and represents the basis for robust SFR maps using matched resolution extinction correction (Pessa et al. 2021).

Fig. 26. Mapping how the pixel-based BPT classifications change as a function of physical resolution for NGC 1365 and NGC 4254. We convolve line maps of each target from copt resolution to 150 pc and 500 pc scales and apply the Kauffmann et al. (2003) [N II]/Hα demarcation to distinguish photoionisation (blue) from harder ionising sources (shocks and AGN, in red). At high spatial resolution, the DIG also displays line ratios consistent with such shock models (Belfiore et al. 2022). Hα contours at the copt resolution (left panels) are shown at each scale for reference. In NGC 4254, as we convolve to larger spatial scales we see that the luminosity weighting results in a larger area of the map appear consistent with photoionisation. This is more apparent in a zoomed-in view (bottom panel). Individual isolated sources with shock excitation (likely SNRs) and the imprint of the DIG are no longer distinguishable in such disc galaxies as we move to the largest 500 pc scales. In contrast, the strong, extended AGN contamination in NGC 1365 can still be recovered. These figures demonstrate the power of diagnostic line ratios in distinguishing ionisation sources, as well as the importance of high (< 150 pc resolution) imaging to obtain an accurate census of ionising sources in the disc.

8.2. Emission line diagnostics

As seen in Fig. 20, the detection of Hα emission across nearly all pixels in our mosaics is complemented by a high (> 50%) detection fraction for a wide array of high and low-ionisation species of oxygen, nitrogen and sulphur atoms even at the level of individual pixels. With binning, we can achieve near complete detection of these strong emission lines across our mosaics, providing new insights into not just the brightest regions but also into the physical conditions impacting the diffuse gas (Belfiore et al. 2022). In addition, we achieve direct detection of weak, temperature-sensitive auroral emission lines across the ∼5% of pixels that correspond to the centres of the brightest H II regions (Ho et al. 2019). While those detections may be biased towards higher-metallicity regions (as for the metallicity-dependent [N II]λ5754 line), line ratios such as [N II]λ5754/[N II]λ6548 and [S III]λ6312/[S III]λ9058 represent relevant insights on the local conditions. This wealth of emission lines enables a variety of diagnostics that constrain the physical properties of the ionised ISM, including the temperature, density, abundances, pressure, and dust extinction (Kreckel et al. 2019, 2020; Barnes et al. 2021).

Well-established line diagnostics (Baldwin et al. 1981, BPT diagrams, named after Baldwin, Philipps and Terlevich) can also distinguish the contribution of different sources of ionisation, for example SNRs, shocks, and planetary nebulae. Many of those sources act on < 100 pc scales and are less luminous than neighbouring photoionised H II regions. As such, they are distinguishable only when sufficiently high spatial resolution observations are obtained. When binning to larger scales, the intrinsic luminosity weighting means that the line ratios often become dominated by the brighter photoionised H II regions, and all spectral information on other ionising sources is lost.

Figure 26 demonstrates how BPT classifications (here using the Kauffmann et al. 2003 demarcation as an example) appear to change with changing spatial resolution, from our native 50−100 pc scales, to 150 pc and 500 pc scales. In the two galaxies shown, the diagnostic maps indicate widespread photoionisation due to star-formation activity (in blue) in addition to ionisation from harder ionising sources (in red) such as shocks and AGN. Diffuse ionised gas also exhibits line ratios consistent with shock excitation, though the exact explanation for these line ratios is still unclear (Zhang et al. 2017). In NGC 1365 the strong central AGN clearly imprints a signal both in the galaxy centre as well as out to large kiloparsec scales (Venturi et al. 2018). In NGC 4254, a flocculent star-forming disc, the red regions are likely due to a combination of DIG emission and SNRs.

Convolving our maps from the copt resolution (50−100 pc) to 150 pc and 500 pc scales, we see that the number and distribution of ionising sources changes with decreasing spatial resolution. Nearly all individual sources in NGC 4254 are lost at 500 pc. However, the impact of the central AGN on NGC 1365 is recoverable even at this lowest physical resolution. These figures demonstrate the importance of our high (< 100 pc) spatial resolution approach in constructing a complete census of ionising sources in the discs of these galaxies.

These images also demonstrate that while weaker emission lines often cannot be detected in pixel-scale maps certain regions; for example, in the outer part of the disc, a detection can be recovered with rather modest binning. In this sense, spatial binning of the original data cubes holds great promise for providing diagnostics for the full disc, and offers leverage for the study of the pervasive DIG (Belfiore et al. 2022).

The multi-wavelength, multi-phase and multi-tracer approach, as adopted by the PHANGS project and illustrated in Fig. 5, has the potential to provide a key link for further dissecting the baryon cycle, probing those complex gaseous plus stellar galactic ecosystems, and placing quantitative constraints on the feedback physics that regulates the chemical and dynamical evolution of the ISM (Chevance et al. 2022; Kim et al. 2021; Barnes et al. 2021). PHANGS-MUSE more specifically enables us to catalogue individual H II regions (Kreckel et al. 2019; Santoro et al. 2022) along with other ionising sources (Planetary Nebulae – PNe –, SNRs; Scheuermann et al., in prep.), crucially informing us about the impact of the individual young clusters catalogued by HST (Turner et al. 2021), and the evolving parent Giant Molecular Clouds (GMCs Sun et al. 2020a,b; Rosolowsky et al. 2021, Hughes et al., in prep.). This cross-observatory approach will provide new insights into stellar feedback (Herrera et al. 2020) and provide fertile testing grounds for future JWST studies of heavily embedded star formation.

Fig. 27. Stellar mass surface density (ΣM⋆) computed from the MUSE spectroscopy as a function of r-band surface brightness (mr) for three representative galaxies (see caption of Fig. 25 and the associated text). All measurements are performed on Voronoi binned regions that reach at least a stellar continuum S/N = 35. For all histograms, bins are required to contain at least three measurements. Top row: binned distribution revealing a clear correlation between mass and light, but with up to an order of magnitude scatter (the colour scale, in log(N), being set up by the point number density). Second row: colouring bins by the median luminosity-weighted stellar age. We demonstrate that the spread is dominated by variations in stellar age. Third row: colouring bins by the median stellar reddening E(B − V) revealing a secondary dependence, reflecting the long known degeneracy between dust extinction and stellar age. Bottom row: colouring bins by the median stellar metallicity illustrating the global scaling with mass as well as additional local variations with both luminosity and stellar mass surface density throughout the disc. We note that very low metallicity values may correspond to an unsolved degeneracy due to bright and very young stellar clusters (see Sect. 6.3.2). As demonstrated in this set of figures, PHANGS-MUSE enables us to infer stellar masses and ages, breaking the long-standing degeneracy between reddening due to dust and reddening due to an ageing stellar population.

8.3. Star-formation histories and dynamical environments

Turning to the stellar continuum, MUSE data enable us to infer stellar masses and ages, breaking the long-standing degeneracy between reddening due to dust and reddening due to an ageing stellar population, going well beyond what is possible from broadband colours. Figure 27 reflects the power of our stellar spectroscopic fits presenting the multi-dimensionality of the PHANGS-MUSE data products. While it is clear that our inferred stellar masses trace the broadband light (top row), the order-of-magnitude spread reflects predominantly the range in stellar ages (second row). This correlation qualitatively dominates over variations driven by stellar dust reddening (third row) or changes in stellar abundance (bottom row). Ongoing work explores how these high-resolution (∼100 pc) constraints on the stellar population relate to traditional scaling relations, such as the resolved stellar main sequence (Pessa et al. 2021). It also provides essential input when carrying out comprehensive modelling of ionising sources in the disc, which in the literature has been suggested to contain a non-negligible contribution from hot evolved main-sequence stars (Zhang et al. 2017; Belfiore et al. 2022).

One other crucial part of the PHANGS-MUSE science case relies on mapping wide areas in order to sample a variety of dynamical environments (spiral arm, bar, centre Querejeta et al. 2021). These trace the individual impact of the galactic potential on the stars and the ionised gas, as well as perturbations from non-axisymmetry (e.g., bar and spiral structures) in the disc (see Fig. 5).

Fig. 28. Changes in stellar velocity dispersion reflected by their dynamical environment. These three galaxies represent galaxies at the lowest stellar mass (NGC 5068), an average stellar mass flocculent morphology (NGC 4254), and the highest stellar mass, AGN-dominated, strongly barred morphology (NGC 1365). Top: each galaxy broken down into its centre (blue), bar (orange), and disc (grey) environments from Querejeta et al. (2021). Bottom: distribution of values from the individual Voronoi bins plotted (see Sect. 5.1), comparing the stellar velocity dispersion (σstars) to the stellar mass surface density (ΣM⋆). We further colour-code each part of the histogram by the environment (second row), the median luminosity-weighted stellar age (third row), or the number density (fourth row). The disc-dominated galaxies (NGC 5068 and NGC 4254) show moderate σstars with an increase at high ΣM⋆, corresponding to the central environment, and younger stellar ages at lower ΣM⋆ (in the spiral arm and outer disc). NGC 1365 presents a more complicated star-formation history, with older stellar ages dominating across the bar associated with increased stellar dispersions. In this set of figures, we see how our PHANGS-MUSE data allow us to trace the individual impact of the galactic potential on the stars, as well as perturbations from non-axisymmetry (e.g., bar and spiral structures) in the disc.

As an illustration, Fig. 28 shows the pixel-based distribution of stellar velocity dispersion as a function of stellar mass surface density (Σ_M⋆) for three of our galaxies (top row). We apply environmental masks (Querejeta et al. 2021) that classify each target into distinct morphological environments: centre, bar, and disc. Environments and galaxy morphologies are shown directly at the top of Fig. 28, while each bin is colour-coded by the median environment (second row). The highest dispersions are associated with high Σ_M⋆, and map to the locations morphologically classified as centre- or bar-dominated, while disc-dominated regions generally exhibit lower velocity dispersions (as expected). We also find good correspondence between stellar ages, mass surface density and dynamics. Colouring each bin by the median luminosity-weighted age (centre row) reveals an expected correlation between younger ages and lower Σ_M⋆ (e.g., at larger radii) in the most disc-dominated galaxies (NGC 5068 and NGC 4254). NGC 1365 presents a more complicated star-formation history, with older stellar ages dominating across the bar. While we find good consistency between expectations based on morphology and the dynamical structure revealed by MUSE kinematics, it also provides new insight into the nature of various dynamical structures. For example, we detect a turn-over in the stellar velocity dispersion towards the centre of NGC 1365, at the highest densities, suggesting a cold, and potentially older, central dynamical structure.

Figures 27 and 28 are meant as a first glimpse at the richness of the PHANGS-MUSE dataset, emphasising the power and beauty (Fig. 5) of two-dimensional spectroscopic mapping of star-forming galactic discs, and the scientific potential that will be developed in future studies by our team and the community via the public data releases.

9. Summary

In this paper we present a detailed account of the PHANGS-MUSE survey, which delivers MUSE/VLT IFS in the optical for a set of 19 nearby star-forming galactic discs. PHANGS-MUSE both covers and resolves a significant fraction of the optical disc, and in particular maps the region enclosing most of its observed molecular gas content (as detected via the CO(2–1) line). The 19 galaxies of the PHANGS-MUSE sample are part of 38 systems surveyed by PHANGS-HST, all drawn from the parent sample of 90 PHANGS-ALMA galaxies, and will be targeted for eight band 2–21 μm imaging as part of the PHANGS-JWST Cycle 1 Treasure programme. This makes such a campaign a unique tool for addressing the physics of various tracers pertaining to the onset of star formation, stellar evolution, and feedback.

We provide a description of the main science objectives for PHANGS-MUSE. They include a refreshed and resolved view on scaling relations, probing the impact of feedback and its relation to the local environment, better understanding the chemical history of galactic discs, and determining the intricate link between dynamical regimes and star formation. The PHANGS-MUSE data can be further exploited on many scientific fronts, representing a huge legacy programme for the scientific community.

Compared to other spectroscopic surveys of galaxies, PHANGS-MUSE stands out for the sheer number of spectra observed. With millions of independent spatial elements, PHANGS-MUSE covers a rich variety of galactic local environments, all at a typical intrinsic resolution of 50 pc. Compared to other two-dimensional IFU galaxy surveys (e.g., MaNGA, SAMI, and CALIFA), we target many fewer galaxies but in far greater detail, achieving a qualitatively different physical resolution. Compared to planned surveys of the Milky Way and Magellanic Clouds (e.g., SDSS-V/LVM), we more uniformly sample the z = 0 galaxy population. And thanks to the multi-wavelength coverage by HST and ALMA at matched or better resolution, we have an unparalleled multi-wavelength view of the baryon cycle within nearby massive star-forming main-sequence galaxies.

We detail the pipeline flows associated with the data reduction and analysis, demonstrating a robust albeit challenging setup to address our ambitious PHANGS-MUSE dataset. We emphasise the modularity of our setup, strongly motivated by the need to re-process and re-analyse the entire set of raw and reduced data several times, driven by, for example, algorithmic improvements and tests. The reduction pipeline is heavily based on the MUSE DRS routines, wrapped up via the pymusepipe package, while for our DAP we further tuned a PHANGS-specific package inherited from the gist package, focusing on the use of, for example, pPXF for all major fitting stages (stellar kinematics, emission lines, and stellar populations). Challenges on the data reduction side include the sky subtraction, accurate astrometric solution for all individual exposures, and post-processed homogenisation of the PSF for each mosaic: these are all addressed via specific and simple algorithms, sometimes supported by existing software (e.g., pypher or mpdaf). The demanding steps for the data analysis are mostly associated with the individual binning schemes and the extraction of the star-formation-history (and associated extinction) maps.

We then provide quantitative measures of the quality of the MUSE-PHANGS data when it comes to wavelength, astrometry, and photometric calibrations. We also give a detailed account of the quality of the spectral fits that drive the delivered science maps, demonstrating that the fits are overall of excellent quality, while regions heavily contaminated by foreground stars or reflecting the presence of very young sources may need further attention.

The outcome of the data reduction and analysis processes are formatted into a PHANGS-MUSE public data release, namely DR1.0, which includes: full native data cubes of the 19 galaxy mosaics, encompassing both the data spectra and their estimated variances, their reconstructed images for a list of standard broad- and medium-band filters, and a large set of derived maps representing extracted information pertaining to the stellar and gas kinematics, emission line gas, stellar and gas extinction, and stellar populations. Earlier versions of the PHANGS-MUSE dataset were distributed internally and exploited in a series of published works. Going beyond PHANGS-MUSE DR1.0, subsequent data releases will potentially include a number of key improvements concerning, for example, the photometry, astrometry, and parameter fitting.

We close the presentation of the PHANGS-MUSE survey by illustrating the science potential of this new observational window on nearby star-forming discs. We provide brief accounts of the Hα pixel-based surface brightness distribution and how dust correction impacts the associated distribution of line fluxes and ratios (Sect. 8.1), the importance of resolving regions with different emission line regimes (Sect. 8.2), and the variation in gas physical states, dynamical states, and star-formation histories among various dynamical environments within discs (Sect. 8.3). These represent only a glimpse of the richness of the PHANGS-MUSE dataset: we believe that its (first and subsequent) public releases will help in developing a better understanding of how the baryon cycle proceeds and how it connects with nearby galaxy disc hosts and their dynamical structures.

Acknowledgments

This work was carried out as part of the PHANGS collaboration. E. S., F. S., R. M. E., T. G. W. acknowledge funding from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation programme (grant agreement No. 694343). S. C. O. G. and R. S. K. acknowledge financial support from the German Research Foundation (DFG) via the collaborative research centre (SFB 881, Project-ID 138713538) “The Milky Way System” (subprojects A1, B1, B2, and B8). They also acknowledge funding from the Heidelberg Cluster of Excellence “STRUCTURES” in the framework of Germany’s Excellence Strategy (grant EXC-2181/1, Project-ID 390900948) and from the European Research Council via the ERC Synergy Grant “ECOGAL” (grant 855130). F. Bi, A. T. B., I. B. acknowledge funding from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation programme (grant agreement No. 726384/Empire). C. E. acknowledges funding from the Deutsche Forschungsgemeinschaft (DFG) Sachbeihilfe, grant number BI1546/3-1. E. J. W. acknowledges support from the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) – Project-ID 138713538 – SFB 881 (“The Milky Way System”, subproject P2). J. M. D. K. and M. C. gratefully acknowledge funding from the Deutsche Forschungsgemeinschaft (DFG) through an Emmy Noether Research Group (grant number KR4801/1-1) and the DFG Sachbeihilfe (grant number KR4801/2-1) and from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation programme via the ERC Starting Grant MUSTANG (grant agreement number 714907). E. C. acknowledges support from ANID project Basal AFB-170002. K. K., O. E., and F. S. gratefully acknowledge funding from the German Research Foundation (DFG) in the form of an Emmy Noether Research Group (grant number KR4598/2-1, PI: Kreckel). E. W. acknowledges support from the DFG via SFB 881 “The Milky Way System” (Project-ID 138713538; sub-project P2). M. B. gratefully acknowledges support by the ANID BASAL project FB210003. J. Pe acknowledges support by the Programme National “Physique et Chimie du Milieu Interstellaire” (PCMI) of CNRS/INSU with INC/INP, co-funded by CEA and CNES. Based on observations collected at the European Southern Observatory under ESO programmes 1100.B-0651, 095.C-0473, and 094.C-0623 (PHANGS-MUSE; PI: Schinnerer), as well as 094.B-0321 (MAGNUM; PI: Marconi), 099.B-0242, 0100.B-0116, 098.B-0551 (MAD; PI: Carollo) and 097.B-0640 (TIMER; PI: Gadotti). This paper makes use of data gathered with the 2.5 meter du Pont located at Las Campanas Observatory, Chile, and data based on observations carried out at the MPG 2.2m telescope on La Silla, Chile. J. C. L. acknowledges support for PHANGS-HST (program number 15654), provided through a grant from the Space Telescope Science Institute under NASA contract NAS5-26555. PHANGS-HST is based on observations made with the NASA/ESA Hubble Space Telescope, and obtained from the data archive at STScI. STScI is operated by the Association of Universities for Research in Astronomy, Inc. under NASA contract NAS 5-26555.

References

Appendix A: Contributions

The PHANGS–MUSE Survey campaign and the associated data products are the outcome of a large team effort, with major contributions from many and input from the entire team. The VLT/MUSE Large Programme was led by E. Schinnerer. The current paper has benefited from work and insights from many members of the team. In the following, we summarise the key contributions that brought us from the early discussions in 2015, to this point.

Observation Design, Data Processing, and Quality Assurance of the MUSE Data: The specific efforts required for designing the MUSE science strategy and observations, its implementation, monitoring, data gathering, data reduction and analysis, quality checks was channelled through a MUSE-focused Working Group within the broader PHANGS team, led by E. Emsellem since its implementation in 2017. The PHANGS–MUSE working group included F. Belfiore, G. Blanc, E. Congiu, E. Emsellem, I-T. Ho, B. Groves, K. Kreckel, R. McElroy, I. Pessa, P. Sanchez-Blazquez, F. Santoro, E. Schinnerer. The design and submission of the MUSE Observation Blocks to ESO Portal were handled by B. Groves, R. McElroy, and C. Faesi, with help from E. Emsellem. The observing campaign monitoring and book-keeping were conducted by R. McElroy and F. Santoro, who also led the uploading and management of the raw data on the Heidelberg computer nodes, with contributions from E. Emsellem and E. Schinnerer. An initial shell version of a MUSE dataflow was written by G. Blanc and B. Groves. The pymusepipe package was written by E. Emsellem, with team contributions by E. Congiu and F. Santoro. The specific alignment process for the MUSE mosaicking was conducted in parallel by E. Emsellem, I-T. Ho, R. McElroy, K. Kreckel and F. Santoro. The DAPwas developed by F. Belfiore, initially based on the gist package (Bittner et al. 2019), with contributions by E. Emsellem and I. Pessa (an early version of the data analysis framework was written and implemented by I-T. Ho). Extensive tests pertaining to the extraction of the stellar population information and stellar extinction were conducted by I. Pessa with support from F. Belfiore, E. Emsellem, I-T. Ho and P. Sanchez-Blazquez. An early version of the datasets were reduced by R. McElroy, and analysed by I-T Ho and K. Kreckel using LZIFU (Ho et al. 2016). The subsequent releases were reduced mainly by E. Emsellem and F. Santoro, with the efforts on the analysis led by F. Belfiore and I. Pessa. All data organisation on the MPIA computers was led by R. McElroy from 2017 to 2018, and then by F. Santoro.

Management of the PHANGS Collaboration: E. Schinnerer has served as the leader of the PHANGS collaboration since 2015. G. Blanc, E. Emsellem, A. Leroy, and E. Rosolowsky have acted as the PHANGS steering committee. E. Rosolowsky has served as team manager since 2018. The PHANGS ‘core team’ provides key input and oversight to all major collaboration decisions. The core team includes: F. Bigiel, G. Blanc, E. Emsellem, A. Escala, B. Groves, A. Hughes, K. Kreckel, J.M.D. Kruijssen, J. Lee, A. Leroy, S. Meidt, M. Querejeta, J. Pety, E. Rosolowsky, P. Sanchez-Blazquez, K. Sandstrom, E. Schinnerer, A. Schruba, and A. Usero. Scientific exploitation of PHANGS–MUSE has taken place largely in the context of the ‘Ionised ISM and its Relation to Star Formation’ (IonisedISM) and ‘Large Scale Dynamical Processes’ (Dynamics) science working groups. The IonisedISM group was led by B. Groves and K. Sandstrom in 2019, and since 2020 by B. Groves and K. Kreckel. The Dynamics group has been led by S. Meidt and M. Querejeta since 2019.

Preparation of this Paper: E. Emsellem has managed the preparation of the text and figures in this paper in close collaboration with E. Schinnerer. E. Emsellem wrote Sect 1, 3.1, 4.1, 4.2.1-4.2.5, 6.1.3-6.1.5, 6.3.1, 7, and 9, provided Figs. 2, 3, 6, 7, 8, 9, 10, 14, 16, 23, made heavy reviews of all sections until finalisation. Major contributions to the paper were provided by F. Belfiore, K. Kreckel and F. Santoro. F. Belfiore wrote Sects. 5.1, 5.2, 6.2.1, 6.2.2 and generated Figs. 2, 3, 4, 11, 12, 13, 17, 18, 19, 21, Tables 1, 2, 3. I. Pessa provided text and figures for stellar population related subsections (6.2.4, 6.3.2; Fig. 21, 22, 24). E. Congiu provided text for Sect 3.2, 4.2.6, 4.2.7, 4.2.8, 6.1.2 and the basis for Fig 7. B. Groves provided heavy input and edits for the sample and science goals and the quality assessment Sections (2.2), and generated and wrote the section about star masks (5.3). F. Santoro provided text in Sect 4, Fig. 5 and Table A.1. K. Kreckel provided text in 6.1.1, 6.1.3, 6.2.3, 8 and Fig. 15, 20, 25, 27, 28. Beyond these, many in the team provided significant input, text or figures, including: G. Blanc, O. Egorov (Fig. 26), R. Klessen, P. Sanchez-Blazquez.

Observatory and Community Support: PHANGS–MUSE would not exist without the key contributions from many staff at the European Southern Observatory, and the expert and unfailing support from the Observatory. This includes all staff involved in receiving and dealing with OBs, the Observatory Staff in Chile, the User Support Department staff associated with MUSE as well as all the behind-the-scene supporting people whose efforts and dedication made the PHANGS-MUSE Large Programme a reality. We especially thank Elena Valenti of the ESO user-support department for her continuous support. We also wish to acknowledge specific discussions and superb software-related support from Ralf Palsa. Beyond ESO, our programme benefited greatly from insights and support from several members of the community. We acknowledge the specific contributions from individuals in the MUSE GTO team, including Roland Bacon, Simon Conseil, Laure Piqueras, and Peter Weilbacher, who were or are key actors in the development and writing of mpdaf and MUSE DRS, and provided invaluable insights.

All Tables

All Figures

Fig. 1. From kpc to ∼100 pc scale, red-green-blue (RGB) colour images (channels derived from reconstructed MUSE mosaic using the SDSS i, r, and g bands) of NGC 4303 (D = 17 ± 3 Mpc) at various spatial resolutions. The four columns, from left to right: at 64 pc (the homogenised resolution for our MUSE dataset) and then convolved to 250, 500, and 1000 pc. The size of the beam (FWHM) is provided as a hatched circle at the bottom-left corner of each top row panel. The bottom row shows zoomed-in views of the top row images (area shown as a white rectangle).

In the text

Fig. 2. Overview of large spectroscopic surveys of nearby galaxies. Left: large spectroscopic surveys in the plane defined by their spatial resolution (in physical units) and the number of spatial resolution elements surveyed. IFU surveys are shown with diamond symbols (green for those achieving cloud-scale or better resolution and blue if probing at kiloparsec scales). VENGA (Blanc et al. 2013) is not shown because it features fewer than 105 resolution elements. For reference, the single fibre SDSS survey (Abazajian et al. 2009, red triangle) is added. PHANGS-MUSE sits in the top-left of this space, ranking highly on both metrics. Right: large spectroscopic surveys in the plane defined by their spatial resolution (in physical units) and the number of galaxies surveyed. PHANGS-MUSE lies on the overall trend line of other IFU surveys.

In the text

Fig. 3. PHANGS-MUSE sample in the M⋆ − SFR plane. Left: PHANGS sample compared with the population of local galaxies from z0MGS (Leroy et al. 2019, small grey dots). The large red circles represent the PHANGS-MUSE galaxies. We show the overlap with the ALMA (blue dots) and HST (black empty squares) components of the PHANGS project. The dashed line is the best fit to the SFMS from Leroy et al. (2019). Right: PHANGS-MUSE sample compared to two complementary projects, EDGE-CALIFA (Bolatto et al. 2017) and ALMaQUEST (Lin et al. 2020), which also target local galaxies with optical IFS and CO interferometric mapping. The dashed line is the best fit to the SFMS from Leroy et al. (2019) with associated scatter (grey shaded area).

In the text

Fig. 4. Synoptic view of the PHANGS multi-wavelength data using NGC 3351 for illustration. Left: blue and red contours showing the footprints of the PHANGS-ALMA and PHANGS-MUSE data. Within the respective footprints, we show flux maps for Hα from MUSE (light red) and CO(2–1) from ALMA (light blue). The footprint of the HST observations included in the PHANGS-HST data release is also shown as an orange contour. Right: HST imaging (F814W filter) shown in colour. The image covers a larger FoV with respect to the left panel in order to show the entire area imaged by HST. We also indicate several radial metrics: the disc scale length (Rd, as derived in Leroy et al. 2021a) and R25. The ALMA and MUSE footprints are also shown, same as in the left panel. We note that all MUSE, ALMA, and HST footprints of the 19 galaxies can be found at https://archive.stsci.edu/hlsp/phangs-hst.

In the text

Fig. 5. Multi-wavelength, multi-phase view of NGC 4535. The top panels present the stellar and gas distribution and kinematics. Top, second panel from the left: multi-emission line view (Hα in red, [O III] in blue, [S II] in green) tracing the sites of massive star formation along the spiral pattern, with differences in the combined colours highlighting the changes in local physical condition (e.g., abundance, ionisation parameter) and ionising source. Second panel from the left: extracted spectra (marked as white circles) demonstrating the typical characteristics of four (numbered) regions, namely: 1- dominated by stellar continuum (red); 2- AGN (orange); 3- H II regions (light purple), and 4- supernova remnants (SNR; purple). AGN line emission is superimposed on the strong stellar continuum absorption, with distinctive strong [O III] and [N II] emission. In the outer disc, H II regions and SNRs have less contribution from the stellar continuum. Broadened line shapes are apparent in the expanding SNR, in contrast to the narrower H II region line emission, and show strong [S II] and [N II] relative to Hα. Zooming into one section of the spiral arm (white box), spatial offsets between the HST star clusters and ionised gas (left top) and between the ionised gas and ALMA molecular gas (left bottom) demonstrate a time sequence evolution across the spiral pattern. Centre right: the gas and stellar velocity fields, mapped through the MUSE spectroscopy, highlight deviations from regular rotation and indicate dynamically driven radial flows along the spiral and bar structures. Right: the stellar mass surface density (ΣM⋆; see Sect. 5.2.4) highlights the location of these dynamical spiral and bar structures, and provides crucial constraints on the underlying gravitational potential affecting all stellar and gaseous processes in the disc.

In the text

Fig. 6. Footprints for the MUSE observations of PHANGS galaxies. Each panel represents one target of the PHANGS-MUSE sample, with a 5 × 5 arcmin2 FoV from the WFI Rc-band images (r-band du Pont for NGC 7496), and the footprints of the MUSE exposures are overlaid in red. Pointings marked with the ⌀ symbol (in NGC 1365, NGC 1512, NGC 1566, NGC 2835, and NGC 3351) and outlined in blue correspond to observations acquired outside of the main PHANGS campaign but reduced following the same data flow and released as part of PHANGS-MUSE. The vertical white bar on the left side of each panel indicates a scale of 5 kpc.

In the text

Fig. 20. Fraction of 0.2″ spaxels inside of 0.5 R25 that have 3σ detections above a given surface brightness threshold (SB) for a representative sample of emission lines. Hα is detected in upwards of 95% of all pixels in most galaxies, with a lower ∼80% fraction in some strongly barred systems (e.g., NGC 1300, NGC 1433, and NGC 1512). Typical 3σ flux sensitivity in Hα is 4−7 × 1037 erg s−1 kpc−2 (3−7 × 10−19 erg s−1 cm−2 per 0.2″ spaxel). Hβ is typically detected in 50−80% of spaxels. Low-ionisation lines ([N II]6584, [S II]6717, and [S II]6731) are detected in 60−95% of spaxels, while the high-ionisation [O III]5007 line emission is less common (50−70% of pixels). For contrast, the faint auroral [N II]5754 emission line is detected at a 3σ level in ∼5% of spaxels (though fewer than 1% are detected at 5σ).

In the text

Fig. 7. Schematic data flow representing the PHANGS-MUSE data reduction part (DRP) of the pipeline. pymusepipe is used as the global process wrapper. The data flow includes a data organiser that addresses the raw data files, followed by standard MUSE data reduction steps (via MUSE DRS), which leads to the first reduced ‘Pixel Tables’ (PixTables) and cubes. These are aligned using WFI or Direct CCD reference images (via reproject). Aligned individual PixTables are then either resampled on a common grid or directly mosaicked. Optional convolution is performed using the framework from mpdaf and astropy, while pypher provides the required kernel cube.

In the text

Fig. 8. Distribution of the FWHM of the MUSE PSF. Top panel: FWHM per pointing as measured by the algorithm described in Sect. 4.2.6. The histograms respectively identify the pointings observed in AO mode (orange) and noAO mode (purple). The red dotted histogram represents the distribution for all existing pointings, while the dashed black histogram shows the distribution of PSFs of the copt data cubes (i.e. after the convolution with the worst observed PSF for each galaxy). Bottom panel: average FWHM per galaxy and the convolved FWHM (again, using the worst PSF in each target as the baseline), with the same colour scheme as in panel (a). The arrows indicate the median value of each sub-sample (AO, noAO, all exposures, convolved) with the associated colour scheme.

In the text

Fig. 9. PHANGS-MUSE reconstructed images. RGB images combining the three i, r, and g (red, green, and blue channels, respectively) SDSS-band reconstructed maps of the 19 targets that are full mosaics from the PHANGS-MUSE Large Programme (including the pilot galaxy NGC 0628) – produced as part of the DRP+DAP flow. The scale in arcminutes, which roughly corresponds to the extent of one single MUSE pointing, is shown at the top left of the figure as a labelled black arrow, while the orange vertical bar on the left side of each panel represents the distance-dependent 5 kpc scale for each target.

In the text

Fig. 10. Analysis flow for the PHANGS-MUSE DAP. After spectral re-binning using a ln(λ) prescription, the resulting mosaicked cubes are used for the extraction of the stellar kinematics, stellar extinction, and stellar population maps, using an appropriate adaptive binning scheme, while the emission line maps (including the corresponding gas kinematics) are derived from the originally sampled cubes. The main information pertaining to the stellar libraries, spectral masks, assumptions and fixed input for each step are mentioned in the upper part of the figure, while the data products are illustrated in its lower part. See Sect. 5.2 for further details.

In the text

Fig. 11. Histograms (violin plots) of the distribution of Voronoi bin sizes for the galaxies in our sample. In orange we show bins within 0.2 R25 and in blue those at larger galactocentric radii. Galaxies are ordered from left to right by stellar mass. The median bin size outside the central regions (R >  0.2 R25) is between 10 and 100 spaxels, roughly corresponding to circularised radii of 0.4″ and 1.1″. The 25th, 50th, and 75th percentiles of each distribution are marked with dashed lines.

In the text

Fig. 12. Maps of the average noise (left panel), S/N (middle panel), and binning map (right panel) for NGC 4535. The noise and S/N maps are computed by averaging the pipeline noise and flux over the 5300−5500 Å wavelength range. The stripes dividing the surveyed area into six squared subregions correspond to the overlap regions of the six MUSE pointings obtained for this galaxy. The noise map also shows an evident cross-hatch pattern within individual pointings due to the cube-generating resampling step in the MUSE DRP when combining exposures with different rotation angles. This behaviour is also visible in the S/N map but does not significantly affect the results of the binning process. The binning map shows the result of the Voronoi binning procedure with target S/N = 35. The black contour shows the S/N = 12 level on individual spaxels. In this target only the galaxy centre (and a few foreground stars) have S/N >  35 in individual spaxels, which are therefore left unbinned.

In the text

	Fig. 13. Comparison of the spectral resolution of the MUSE data (Eq. (3), taken from Bacon et al. 2017) and the E-MILES stellar templates. The greyed-out part of the wavelength range (λ >  7000 Å) was excluded from the fit in this first DR1.0 PHANGS-MUSE data release.
In the text

	Fig. 14. RGB images of the 19 targets in the MUSE sample obtained by combining the emission line maps for [S II], Hα, and [O III] (green, red, and blue channels, respectively). A 5 kpc double black arrow indicates the fixed physical scale for all panels, just below the RGB labels respectively associated with each emission line map.
In the text

Fig. 15. Comparison of r-band magnitudes measured over 5″ × 5″ apertures within MUSE synthetic images (rMUSE) and SDSS images (rSDSS) for the nine galaxies with SDSS imaging available. Each panel is separated into a main one-to-one comparison plot and an additional plot below showing the magnitude difference (ΔrMUSE − SDSS) versus the SDSS value. Insets in each panel show the histograms of the magnitude difference, with the median offset indicated with a dashed line. Solid lines indicate the identity line (or the 0 value for the magnitude difference). Across this sample of galaxies, the median offset ranges from −0.06 to 0.01 mag. Typical scatter within any galaxy is ∼0.04 mag.

In the text

Fig. 16. Median offset and percentiles of the difference between spectra of the same regions from overlapping pointings. The statistics have been computed from a set of 249 regions covering the full PHANGS-MUSE sample, and the resulting percentile vectors have been filtered to make it legible (keeping the sky line residuals visible). Top panel: percentiles of the distribution of the bias level, normalised by the typical (overlap-region-averaged) noise level. Ninety percent of the spectra have a normalised bias that is typically between 30 and 50% of the noise level, while a small fraction show up at levels of 60−120% of the noise, specifically in the blue or red part of the spectrum. We note that although beyond 7000 Å residuals are heavily contaminated by sky line residuals, the pipeline still constructs a roughly correct noise vector. Bottom panel: percentiles (50, 68.3, 95.5%, and 99.7%) of the distribution of differences normalised by the individual spectra noise level, after subtraction of a wavelength-independent median bias offset (see text). The dashed (respectively, green, yellow, and red) lines show values of 1.12 (12% above 1), 2.24, and 3.36, showing that the noise level is slightly underestimated (by about 12%). A trend is visible towards the redder and bluer end of the wavelength coverage.

In the text

Fig. 23. Stellar and gas E(B − V) for three example galaxies. The stellar E(B − V)stars is determined from the Voronoi binned data cubes used for stellar population analysis, while the gas E(B − V)gas is computed from the Balmer decrement, derived from the spaxel-level analysis of the ionised gas emission lines (only regions with S/N >  4 in the line emission are considered). The colour bars for the gaseous E(B − V) is stretched by a factor of two to account for the average ratio E(B − V)stars ∼ 0.5 E(B − V)gas. E(B − V)stars traces similar structures as E(B − V)gas, although jumps are evident when comparing the average level of some MUSE pointings to that of their neighbours (e.g., the bottom-right pointing in NGC 4535 and the three left-most pointings in NGC 4254).

In the text

Fig. 17. Histograms of the distributions of reduced chi-square, ${\chi }_{\mathrm{d}.\mathrm{o}.\mathrm{f}.}^2$, of the resulting spectral fits for all the spaxels in all the galaxies in our sample. The reduced χ2 is calculated over the wavelength range fitted by the DAP (4850−7000 Å). The red horizontal line and shaded area corresponds to the median (16th and 84th) percentiles of the distribution across all galaxies, corresponding to ${\chi }_{\mathrm{d}.\mathrm{o}.\mathrm{f}.}^2=1.{20}_{-0.18}^{+0.23}$.

In the text

Fig. 18. Maps of the ${\chi }_{\mathrm{d}.\mathrm{o}.\mathrm{f}.}^2$ for three example galaxies at the individual spaxel level. We show in blue contours the Hα surface brightness and in red the regions masked as foreground stars. The ${\chi }_{\mathrm{d}.\mathrm{o}.\mathrm{f}.}^2$ maps demonstrate the quality of our spectral fits ${\chi }_{\mathrm{d}.\mathrm{o}.\mathrm{f}.}^2\sim 1$, except in a few regions of bright Hα emission (e.g., the central regions in NGC 1672 and NGC 4535). A few regions of higher ${\chi }_{\mathrm{d}.\mathrm{o}.\mathrm{f}.}^2$ can also be seen in the middle of spiral arms, corresponding to the centres of bright H II regions. Foreground stars appear as small spots of very high ${\chi }_{\mathrm{d}.\mathrm{o}.\mathrm{f}.}^2$ and are emphasised with red contours.

In the text

Fig. 19. Histogram of the distributions of ${\chi }_{\mathrm{d}.\mathrm{o}.\mathrm{f}.}^2$ (in blue) for the same example galaxies as in Fig. 18. The median ${\chi }_{\mathrm{d}.\mathrm{o}.\mathrm{f}.}^2$ is shown as a blue vertical dashed line. We show as a bar plot in orange and green the stacked histograms of pixels belonging to our foreground star masks (orange) and to the brightest line-emitting regions (> 99th percentile of the Hα emission, in green). Because the bar plots are stacked, their sum must be smaller than the histogram of all spaxels in blue. This allows us to visually assess that the vast majority of high-${\chi }_{\mathrm{d}.\mathrm{o}.\mathrm{f}.}^2$ pixels are due to foreground stars and bright line-emitting regions.

In the text

Fig. 21. Two-dimensional histogram of the ratio between the stellar mass surface density derived from full spectral fitting from MUSE (Sect. 5.2.4) and archival stellar mass maps from S4G, based on Spitzer IRAC maps (Querejeta et al. 2015), as a function of the light-weighted age of the underlying stellar population in each pixel, derived from the full spectral fitting. We show three galaxies in our sample, spanning the stellar mass range of the PHANGS-MUSE sample (in ascending order from left to right). The figure shows good agreement between the MUSE and the S4G mass maps in pixels dominated by stellar populations older than ∼4 Gyr. However, in regions hosting younger populations, the masses derived from the near-IR data are systematically larger, possibly due to the contribution of AGB stars to the near-IR flux. The contours mark the iso-probability curves where the probability density drops below 40% and 10% of the maximum. The dashed red line shows the median mass ratio at a given luminosity-weighted age.

In the text

Fig. 22. Age dependence of the MMUSE/L3.6 μm ratio (in M⊙/L⊙) for the galaxies in our sample, where L3.6 μm corresponds to the luminosity of a pixel in the original IRAC 3.6 μm image (not the ICA version from S4G). Each line represents a galaxy, coloured by its total stellar mass. The figure shows a positive trend, with pixels hosting older stellar populations having larger M/L ratios. The horizontal lines mark different values adopted in the literature (dot-dashed line M/L = 0.41, Schombert et al. 2019; dotted line, M/L = 0.53, Eskew et al. 2012; dashed line M/L = 0.6, Meidt et al. 2014).

In the text

Fig. 24. Age (top row), metallicity (middle row), and extinction (bottom row) of the inner star-forming ring of NGC 3351 obtained using three different sets of templates: E-MILES (left column), CB07, excluding templates younger than 30 Myr (CB07-A; central column), and CB07 including young templates (CB09-B; right column). Adding younger templates to our template grid slightly alleviates the problem of extremely metal-poor and young regions in the sense that the young clusters are younger and mildly less metal-poor. However, the improvement is marginal, and these regions are still significantly more metal-poor than the surrounding pixels. Panels in the last row show, additionally, an abnormally low extinction in the young and metal-poor regions (marked as white circles in the left column). This is likely a consequence of the same issue, as it is partially improved when younger templates are included. Overall, adding templates younger than 30 Myr to our age-metallicity grid does not provide a solution to the issue of young and extremely metal-poor regions.

In the text

Fig. 25. Distribution of pixel-based Hα/Hβ line ratios as a function of observed (top) and extinction corrected (bottom) Hα surface brightness (SBHα) for three galaxies of our sample (see text). We employ a S/N cut of 10 for both lines. Most pixels have Hα/Hβ >  2.86 (dashed line, as predicted by case B recombination for Te = 104 K, ne = 100 cm−3), and those that have ratios below this value are consistent within the errors. Because of this, when de-reddening Hα we apply no correction for those pixels with Hα/Hβ < 2.86. At low Hα surface brightness, our Hβ sensitivity sets an upper threshold on the recoverable value of Hα/Hβ. The value of this threshold for the median Hβ error is indicated in all panels by the dotted lines. Histograms of SBHα are included for each set of plots for pixels that meet our S/N cut (in black) and those that do not (in grey). A significant number of pixels are excluded by requiring an S/N >  10 detection of Hβ. A wide range of extinctions is observed, with no obvious correlation with SBHα. This demonstrates the need for matched resolution constraints on the dust extinction. Within PHANGS-MUSE, the use of the Balmer decrement provides a unique opportunity to perform robust, extinction-corrected SFR measurements across the bulk of the galaxy disc.

In the text

Fig. 26. Mapping how the pixel-based BPT classifications change as a function of physical resolution for NGC 1365 and NGC 4254. We convolve line maps of each target from copt resolution to 150 pc and 500 pc scales and apply the Kauffmann et al. (2003) [N II]/Hα demarcation to distinguish photoionisation (blue) from harder ionising sources (shocks and AGN, in red). At high spatial resolution, the DIG also displays line ratios consistent with such shock models (Belfiore et al. 2022). Hα contours at the copt resolution (left panels) are shown at each scale for reference. In NGC 4254, as we convolve to larger spatial scales we see that the luminosity weighting results in a larger area of the map appear consistent with photoionisation. This is more apparent in a zoomed-in view (bottom panel). Individual isolated sources with shock excitation (likely SNRs) and the imprint of the DIG are no longer distinguishable in such disc galaxies as we move to the largest 500 pc scales. In contrast, the strong, extended AGN contamination in NGC 1365 can still be recovered. These figures demonstrate the power of diagnostic line ratios in distinguishing ionisation sources, as well as the importance of high (< 150 pc resolution) imaging to obtain an accurate census of ionising sources in the disc.

In the text

Fig. 27. Stellar mass surface density (ΣM⋆) computed from the MUSE spectroscopy as a function of r-band surface brightness (mr) for three representative galaxies (see caption of Fig. 25 and the associated text). All measurements are performed on Voronoi binned regions that reach at least a stellar continuum S/N = 35. For all histograms, bins are required to contain at least three measurements. Top row: binned distribution revealing a clear correlation between mass and light, but with up to an order of magnitude scatter (the colour scale, in log(N), being set up by the point number density). Second row: colouring bins by the median luminosity-weighted stellar age. We demonstrate that the spread is dominated by variations in stellar age. Third row: colouring bins by the median stellar reddening E(B − V) revealing a secondary dependence, reflecting the long known degeneracy between dust extinction and stellar age. Bottom row: colouring bins by the median stellar metallicity illustrating the global scaling with mass as well as additional local variations with both luminosity and stellar mass surface density throughout the disc. We note that very low metallicity values may correspond to an unsolved degeneracy due to bright and very young stellar clusters (see Sect. 6.3.2). As demonstrated in this set of figures, PHANGS-MUSE enables us to infer stellar masses and ages, breaking the long-standing degeneracy between reddening due to dust and reddening due to an ageing stellar population.

In the text

Fig. 28. Changes in stellar velocity dispersion reflected by their dynamical environment. These three galaxies represent galaxies at the lowest stellar mass (NGC 5068), an average stellar mass flocculent morphology (NGC 4254), and the highest stellar mass, AGN-dominated, strongly barred morphology (NGC 1365). Top: each galaxy broken down into its centre (blue), bar (orange), and disc (grey) environments from Querejeta et al. (2021). Bottom: distribution of values from the individual Voronoi bins plotted (see Sect. 5.1), comparing the stellar velocity dispersion (σstars) to the stellar mass surface density (ΣM⋆). We further colour-code each part of the histogram by the environment (second row), the median luminosity-weighted stellar age (third row), or the number density (fourth row). The disc-dominated galaxies (NGC 5068 and NGC 4254) show moderate σstars with an increase at high ΣM⋆, corresponding to the central environment, and younger stellar ages at lower ΣM⋆ (in the spiral arm and outer disc). NGC 1365 presents a more complicated star-formation history, with older stellar ages dominating across the bar associated with increased stellar dispersions. In this set of figures, we see how our PHANGS-MUSE data allow us to trace the individual impact of the galactic potential on the stars, as well as perturbations from non-axisymmetry (e.g., bar and spiral structures) in the disc.

In the text