EDP Sciences logo
Subscriber Authentication Point
Search
Menu
Home All issues Volume 660 (April 2022) A&A, 660 (2022) A77 Full HTML
Free Access
Issue
A&A
Volume 660, April 2022
Article Number A77
Number of page(s) 26
Section Extragalactic astronomy
DOI https://doi.org/10.1051/0004-6361/202142315
Published online 13 April 2022

© ESO 2022

1. Introduction

Young massive stars originate from the gravitational collapse of giant molecular clouds (GMCs), and they inject energy and momentum into the surrounding interstellar medium (ISM) via different feedback processes such as thermal feedback from protostars, photoionisation, and mechanical feedback from stellar winds and supernovae (for a review, see Krumholz et al. 2014 and Dale 2015). These combined effects can disrupt the parent GMCs (Dale 2015; Howard et al. 2017), resulting in a self-regulating mechanism which inhibits future star formation. On the other hand, positive feedback can occasionally facilitate further collapse in neighbouring regions, boosting the formation of new stars. Stellar feedback directly affects star-forming regions by carving channels and ‘bubbles’ of ionised gas into the surrounding cool gas and dust, and it also impacts the ISM on wider galactic scales, as it clears the path for winds and outflows to escape the star-forming regions. Overall, stellar feedback can shape the global properties of a galaxy (e.g., Scannapieco et al. 2012; Hopkins et al. 2013), and it plays a key role in the recycling of gas, the regulation of star formation, and the chemical enrichment and mixing of star-forming galaxies (e.g., Leroy et al. 2015a; Maiolino & Mannucci 2019). Modelling feedback has therefore proven to be essential in simulations of GMCs (e.g., Dale et al. 2014) as well as galaxy formation and evolution (Schaye et al. 2015; Hopkins et al. 2018) in order to correctly reproduce key observables and relations between them.

A more detailed understanding of the process of stellar feedback is also key to explain the origin of a warm and diffuse component of the ISM (diffuse ionised gas or DIG, see Haffner et al. 2009 for a review), which has been observed to make up a considerable fraction of the Hα emission (up to 50%) in local spiral galaxies (Ferguson et al. 1996; Zurita et al. 2000; Thilker et al. 2002; Hoopes & Walterbos 2003; Oey et al. 2007). The origin of this ISM component is still under study, and has been linked to radiation leaking from the star-forming regions (Zurita et al. 2002; Weilbacher et al. 2018; Belfiore et al. 2022; Della Bruna et al. 2021), field stars (Hoopes & Walterbos 2000; Zhang et al. 2017; Belfiore et al. 2022), shocks (Collins & Rand 2001), cosmic rays (Vandenbroucke et al. 2018) or scattering by dust (Seon & Witt 2012). Studying stellar feedback in nearby galaxies will allow to probe the role of each of these processes in ionising the DIG.

The advent of integral field spectroscopy (IFS) has been a leap forward in the study of stellar feedback, as it allows to acquire spectral information of multiple lines simultaneously over a relatively wide field of view (FoV) at good angular resolution. This makes it possible for example to obtain at the same time the kinematics of gas and stars, and to disentangle different gas ionisation mechanisms (see Kewley et al. 2019, for a review).

IFS instruments have led the way to multi-scale studies of kinematics and physical properties of ionised gas in nearby galaxies. At the largest (kpc) scales, surveys such as CALIFA (Sánchez et al. 2012), MANGA (Bundy et al. 2015) and SAMI (Croom et al. 2012; Bryant et al. 2015) provided a census of the z ∼ 0 galaxy population. The large sample of H II regions provided by these surveys allowed for a better understanding of the general properties of star forming regions. For example, it was possible to investigate the link between the physical conditions of the ISM and the stellar population in the regions (Sánchez et al. 2015), and between the stellar population and the leakage of ionising photons (Morisset et al. 2016). Large scale surveys also allowed to quantify and characterise the presence of extraplanar DIG (eDIG) in edge-on spiral galaxies (Bizyaev et al. 2017; Jones et al. 2017; Levy et al. 2019).

At intermediate scales, two instruments whose wide field coverage and high sensitivity have been particularly important for the study of feedback are: the MUSE integral field unit (IFU, Bacon et al. 2010), mounted on the ESO Very Large Telescope (VLT) at Cerro Paranal observatory and the SITELLE imaging fourier transform spectrograph at Canada-France-Hawaii Telescope (CFHT, Drissen et al. 2019). The relatively large FoV (1 arcmin2 and 11 × 11 arcmin2, respectively) and high spatial sampling ( 0 . ″ 2 $ 0{{\overset{\prime\prime}{.}}}2 $ and 0 . ″ 32 $ 0{{\overset{\prime\prime}{.}}}32 $ pixel−1) of these instruments are ideal to study global properties of nearby galaxies, while at the same time achieving a high level of detail at small scales. This was exploited by galaxy surveys such as MAD (Erroz-Ferrer et al. 2019; den Brok et al. 2020), SIGNALS (Rousseau-Nepton et al. 2019) and PHANGS-MUSE (Emsellem et al. 2022), that are targeting each ∼20–50 galaxies, at distances where 1″∼50–200 pc, enabling the identification and characterisation of individual H II regions. This makes it possible to study the dependence of region properties on their location in the galaxy (e.g., in arm vs. interarm regions) and on the local ISM conditions (Kreckel et al. 2016, 2019; Erroz-Ferrer et al. 2019).

Finally, a few specific studies of local galaxies with MUSE (LMC, McLeod et al. 2019; NGC 300, McLeod et al. 2020, 2021; NGC 7793, Della Bruna et al. 2020, 2021) are starting to probe scales of tens of parsec, resolving individual H II regions in detail and detecting the most massive individual stars, which gives a way to bridge small and galaxy-scale stellar feedback. For individual regions, this allows one to investigate, for example, the contribution of various types of feedback (McLeod et al. 2019, 2020, 2021), the ionisation structure of the region and to derive an escape fraction by modelling the expected ionising photon flux based on the observed stellar population (McLeod et al. 2019, 2020; Della Bruna et al. 2020, 2021).

In this study, we exploit MUSE IFS data to study the ISM in the nearby galaxy M83 at a scale of ∼20 pc. Our work complements the very high resolution (10 pc) studies and the large scale (100 pc – kpc) studies, by covering a wide area (∼20 arcmin2) at a resolution that still allows the characterisation of individual H II regions while at the same time sampling a large number of regions (∼4700, Della Bruna et al., in prep.) and a wide range of ISM conditions across the galactic disk.

M83 (also known in the literature as NGC 5236) is a nearby spiral galaxy, at a distance of 4.89 Mpc (Jacobs et al. 2009). Figures 1 and 2 show M83 observed with the Hubble Space Telescope (HST) and the Atacama Large Millimeter/submillimeter Array (ALMA), and with MUSE (this work). M83 has a grand design barred spiral morphology, and its full optical disk is ∼20 kpc across; the position of the stellar bar and the two spiral arms1 are sketched in Fig. 1. The galaxy has a stellar mass log10M⋆ = 10.53 M⊙ (Leroy et al. 2021b) and hosts a nuclear starburst ring (Sérsic & Pastoriza 1965; Buta & Crocker 1993; Calzetti et al. 2004; Comerón et al. 2010; Knapen et al. 2010), resulting in a high star formation rate (SFR = 4.2 M⊙ yr−1, Leroy et al. 2021b). Its offset from the SFR-M⋆ main sequence for star forming galaxies (ΔMS = 0.44 dex) is consistent with the typical scatter at this stellar mass (Popesso et al. 2019). M83 has been extensively studied over all wavelengths, and its stellar population has been mapped across the galactic disk. Catalogues of young star clusters (YSCs, Silva-Villa et al. 2014; Adamo et al. 2015), Wolf-Rayet stars (Hadfield et al. 2005) and supernova remnants (Blair et al. 2014; Winkler et al. 2017; Williams et al. 2019; Russell et al. 2020) are publicly available. Overall, this makes M83 an ideal candidate for the detailed study of stellar feedback at scales of 10 s of parsecs.

thumbnail Fig. 1.

HST colour-composite image of M83 (Red: F657N, Blue: F814W, Green: F438W). The footprint of the MUSE mosaic is overlaid in white. The position of the bar and the two spiral arms is sketched in black (see Sect. 1). The white contours trace molecular gas mass surface density from ALMA CO(2-1) (contours of 5 and 25 K km s−1, corresponding to a surface gas density ≃30 and 170 M⊙ pc−2).

thumbnail Fig. 2.

Three-colour composites of the MUSE data. Left: composite of stellar bands (Blue: 4875–4950 Å, Green: 6520–6528 Å, Red: 6750–6810 Å). The green rectangle indicates the central starburst region analysed in Sect. 6. The orange square indicates the position of the Hα cloud highlighted in Fig. 3 (central panels). Right: composite of gas emission (Blue: [O III], Green: [S II] Red: Hα).

The MUSE data presented in this work consist of a large mosaic covering the central 3.8 kpc in radial extent at ∼20 pc resolution. This is the first extensive IFS study of stars and ionised gas across the disk of M83. Previous studies have either targeted the ionised gas using a wide field spectrograph (Poetrodjojo et al. 2019, ∼40 pc resolution) or Fabry–Pérot observations (Fathi et al. 2008), or have focused on spectroscopy of the nuclear region (e.g., Knapen et al. 2010; Piqueras López et al. 2012; Gadotti et al. 2020; Callanan et al. 2021). We complement our observations with high-resolution archival imaging from HST (< 2 pc resolution) and CO(2-1) observations from ALMA at ∼50 pc resolution. An overview of the dataset is shown in Fig. 1. The MUSE data allow us to map for the first time the large-scale stellar kinematics, as well as unprecedentedly detailed Hα kinematics. The wealth of information provided by IFU spectroscopy enables us to investigate spatially resolved physical properties of the gas. At the same time, the HST coverage provides us with a detailed catalogue of massive stars and clusters, and the ALMA data trace the distribution of the molecular gas in the galaxy. The goal of this project is the study of star formation from galactic scales (this paper) to small scale (∼20 pc, Della Bruna et al., in prep.). The data allow also for a census of Wolf-Rayet stars (Smith et al., in prep.), planetary nebulae (Della Bruna et al., in prep.) and supernova remnants (Long et al., in prep.). Ultimately, we aim to provide a comprehensive picture of the star formation cycle in M83 at tens of parsec scales.

This paper is organised as follows: in Sect. 2 we give an overview of the dataset and of the MUSE data reduction. In Sect. 3 we present the kinematics of stars and ionised gas. In Sect. 4 we estimate the extinction from the ionised gas and compare it with the distribution of the molecular gas. In Sect. 5 we identify H II regions complexes, determine the fraction of DIG and inspect the properties of the ionised gas in emission line diagrams. In Sect. 6 we take a closer look at the kinematics and emission of the central starburst region. Finally, in Sects. 7 and 8 we discuss the results and draw our conclusions.

2. Dataset overview and MUSE data reduction

We observed M83 with the MUSE instrument at VLT (observing programmes 096.B-0057(A) and 0101.B-0727(A), PI Adamo). The main properties of the target are summarised in Table 1. The MUSE data cover the inner ∼3.8 kpc in galactocentric radial extent (1.1 × Re in effective radius, see Table 1), and a total physical area of 40.5 kpc2. The dataset consists of 20 pointings (indicated in green in the top panel of Fig. C.1) in Wide Field Mode (WFM) and extended wavelength configuration (4650–9300 Å), for a total of 46 exposures. The pointings were observed with science exposures of 550 s (4 exposures for pointings 1 to 8, a single exposure for remaining pointings). Sky frames of 180 s were obtained such that every science exposure is preceded or followed by a sky frame acquisition. Because of the large extent of the target in the sky, sky frame offsets were larger than 4 arcmin. We complemented our observations with MUSE data available in the ESO archive: observing programmes 097.B-0899(B) (PI Ibar) and 097.B-0640(A) (PI Gadotti). Both programmes were observed in WFM with nominal wavelength range (4800–9300 Å). From these programmes, we excluded some frames that had calibration issues, poor seeing or that were not fully overlapping with the rest of our dataset. From archival dataset 097.B-0899(B) we included five pointings (indicated in blue in Fig. C.1, top panel), each observed during three science exposures of 600 s with two external sky frames with a 180 s exposure. Archival dataset 097.B-0640(A) consists of one pointing (purple square in Fig. C.1, top panel), of which we included four 480 s exposures, and an additional two sky exposures of 300 s.

Table 1.

Physical properties of M83.

The final mosaic consists of 26 pointings for a total of 65 MUSE exposures. The central coordinates and exposure time of all the included pointings can be found in Appendix C (Table C.1 and Fig. C.1). The pointings were imaged over a wide range of observing conditions. We derived a PSF for each pointing by fitting a Moffat profile as a function of wavelength to bright, isolated point sources with PAMPELMUSE (Kamann 2018), and report the values in Table C.1. Across the full mosaic we measure a median PSF of 0 . ″ 7 $ 0{{\overset{\prime\prime}{.}}}7 $ at 7000 Å (17 pc at the distance of our target). The full width at half maximum of the Moffat profile declines with increasing wavelength, with a median difference of 0 . ″ 16 $ 0{{\overset{\prime\prime}{.}}}16 $ between the blue and red end.

M83 was also observed with HST during the WFC3 Early release science programme (GO11360, PI O’Connell), using narrow and broad band imaging ranging from the UV to the NIR. The coverage of the inner 4.5 kpc (1.3 × Re) in galactocentric radius was later completed with the programme GO12513 (PI Blair). In total, the HST mosaic consists of seven contiguous pointings (Blair et al. 2014), with a FWHM ∼ 0 . ″ 08 $ {FWHM} \sim0{{\overset{\prime\prime}{.}}}08 $ (1.9 pc).

Finally, archival ALMA CO (J = 2–1) data2 covering the inner 7.2 kpc (2.1 × Re) of the galaxy in galactocentric radius have been processed as part of the PHANGS-ALMA survey (Leroy et al. 2021b) using the PHANGS-ALMA pipeline (Leroy et al. 2021a). The data have an angular resolution of 2 . ″ 14 $ 2{{\overset{\prime\prime}{.}}}14 $ (50 pc), a spectral resolution of 2.5 km s−1 channel−1 and a rms brightness temperature sensitivity of 0.17 K. The data are described in detail in Leroy et al. (2021b), and the data reduction in Leroy et al. (2021a).

In Fig. 1 we show a 3-colour composite of the HST data. The footprint of the MUSE mosaic is shown in white. The white contours indicate the ALMA CO emission at levels of 5 and 25 K km s−1. Using Eq. (1) in Sun et al. (2020) and assuming a CO(2-1)-to-CO(1-0) line ratio R21 = 0.65 (den Brok et al. 2021; Leroy et al. 2022) and a standard Galactic CO-to-H2 conversion factor αCO = 4.35 M⊙ pc−2 (K km s−1)−1, these values corresponds to a cold gas surface density of ≃30 and 170 M⊙ pc−2. The position of the bar and the spiral arms is shown with black dashed lines (see Sect. 1). The HST and ALMA data cover the inner ∼4.5–7 kpc of the disk, imaging the spiral arms in their entirety, whereas the MUSE mosaic is limited to the central 3.8 kpc.

We reduced the MUSE data with the ESO pipeline (Weilbacher et al. 2014, 2020, v2.8.3), following the standard reduction procedure. In a first phase, the instrumental signature was removed using muse_scibasic. For this purpose, we combined bias frames to account for the readout noise and lamp frames to correct for flat fielding. Wavelength calibration was performed by combining arc-lamp frames, and correction of 3D illumination and a refined flat-fielding was achieved by combining twilight frames. Moreover, we used the line spread function (LSF) and geometry table (describing the spatial location of IFU slicers) distributed with the pipeline. In a second phase, we performed flux calibration using a standard star, and we modelled the sky spectrum from the sky frames with muse_create_sky. We then ran muse_scipost on the object and sky frames. For the object frame, we used the ‘subtract-model’ sky subtraction method, that uses sky lines and continuum estimated from the sky frames. We saved whitelight images as well as individual pixel table, that we later used when aligning and combining the individual exposures. In a third step, we assembled the final mosaic. We used muse_exp_align to compute offsets between multiple exposures of the same pointing. Offsets between different pointings were instead computed by matching each MUSE whitelight image with the B-band HST mosaic obtained from the MAST archive3, and registered to the Gaia DR2 (Gaia Collaboration 2018). We selected in each MUSE image at least ten point-like sources within the FOV that corresponded to isolated bright clusters in the HST frame. The list of offsets in WCS was then used by muse_exp_combine to combine the individual pixel tables into a large mosaic. A colour composite of the stellar continuum in different broad bands (left) and of three line emission maps (right) extracted from the final MUSE mosaic are shown in Fig. 2.

3. Kinematics of stars and gas

3.1. Stellar kinematics

We fitted the stellar continuum with PPXF (Cappellari & Emsellem 2004; Cappellari 2017), using the E-MILES simple stellar populations models (Vazdekis et al. 2016) with an unimodal IMF and Padova 2000 isochrones (Girardi et al. 2000). The data are spatially binned using the weighted adaptation (Diehl & Statler 2006) of the Voronoi tessellation method (Cappellari & Copin 2003). The bins were targeted to a signal to noise ratio (S/N) ≃ 250 in the continuum range 5025–5065 Å, corresponding to a S/N ≃ 40 Å−1. We fitted the range 4600–8740 Å (in order to include the Ca II triplet). All relevant emission lines, as well as sky residuals were masked during the fit; in particular, we excluded entirely the range 7220–8507 Å due to the strong residual sky emission. We used the Gaussian parametrisation of the instrumental LSF by Guérou et al. (2017)4. We note that following this parametrisation, velocity dispersions below ∼50 km s−1 at Hα are undersampled by the MUSE instrument.

The resulting stellar velocity and dispersion maps are shown in Fig. 3 (top row). We observe that the stars exhibit overall a regular rotation, and that the velocity dispersion increases from σ ≃ 30 km s−1 in the outskirts to σ ≃ 80 km s−1 towards the centre. At the very centre, we observe a fast rotating component and a dip in the velocity dispersion (σ ≃ 50 km s−1). We discuss this feature in Sect. 6.1, where we analyse the kinematics of the starburst region.

thumbnail Fig. 3.

Kinematics of the stars and gas. Top row: MUSE stellar kinematics. Centre row: MUSE Hα kinematics. The black box in the bottom right corner indicates the position of what we interpret as an extraplanar Hα cloud. Bottom row: molecular gas kinematics from ALMA CO(2-1) (mom1 and mom2 maps). The footprint of the MUSE data is overlaid in black. All velocity maps have been corrected for systemic velocity; no inclination correction has been applied. The velocity dispersion maps refer to the intrinsic dispersion (instrumental effects have been removed). The white dashed lines sketch the position of the bar and spiral arms as shown in Fig. 1 (see Sect. 1). A zoom-in into the central kinematics is shown in Fig. 11.

3.2. Ionised gas kinematics

We created a stellar continuum subtracted cube (‘gas cube’) by rescaling the best stellar population fit from PPXF in each Voronoi bin by the height of the continuum in each spatial pixel (spaxel). We studied the kinematics of the ionised gas from the Hα emission line, the brightest line in the MUSE spectral range throughout our FoV. We Voronoi binned the gas cube to a S/N ∼ 20 in Hα5 and fitted the line with a single Gaussian profile. We used the python scipy module curve_fit, with an initial guess (vi, σi) = (vsyst, 20) km s−1 for the velocity and velocity dispersion, where vsyst is the systemic velocity estimated from the stellar kinematics. We set the initial flux fi to the integral of the binned spectrum in a window of width 22 Å (±500 km s−1) centred on the line. We then performed the fit in the same wavelength window, in order to avoid being affected by a poorly subtracted stellar continuum. We constrained f ≥ 0 and σ ≥ 1 km s−1, and leave all other parameters free. We removed the instrumental signature using the parametrisation of the MUSE LSF described in Sect. 3.1. The resulting velocity and dispersion maps are shown in Fig. 3 (middle row panels). We remark that while the spectral resolution of MUSE at Hα is coarse (σHα ≃ 120 km s−1), the S/N on the Hα line yields a median centroiding accuracy of ≃2 km s−1. We also note that the artefacts in both maps are arising from minute calibration uncertainties at the scale of a few 10 s of km s−1 and to undersampling of the MUSE LSF (see Sect. 3.1).

The kinematics of the ionised gas are more complex than the stellar kinematics. At galactic scales, we observe an overall disk rotation, while the velocity dispersion increases in the interarm regions (σ ≃ 80 km s−1) with respect to the star-forming regions along the spiral arms (σ ≃ 20 km s−1). The black square in the central panels of Fig. 3 indicates a compact region that features a peculiar Hα velocity which appears to be 100 km s−1 lower than the surrounding disk rotation, and a high velocity dispersion ∼80 km s−1. Given that this region does not stand out in the stellar RGB image (see orange square in Fig. 2) and that it shows a typical spectrum of a low-luminosity Hα region, we interpret it as a cloud of extraplanar gas, possibly falling into or being ejected from the galactic disk.

We also performed a double Gaussian component fit to the Hα line. A double component can better capture the emission in the case of extreme line broadening on top of a single emission line peak (caused e.g., by shocks or unresolved flows) and in the case where the line has a double-peaked profile, tracing gas motions on top of the galactic rotation (e.g., resolved flows). The fit was performed as following: in a first step, we fitted a single component, as described above. This fit was then adopted as initial guess for the first component, whereas for the second component, we set an initial guess (vi, σi) = (v1comp, σ1comp) and f0 = 0.25 × f1comp. All boundary conditions were the same as described above. In order to prevent over-fitting by the broad component (e.g., in the case of a noisy or poorly subtracted continuum), we limited the fit to a window of width 22 Å (±500 km s−1) around the rest wavelength of Hα. We determined the optimal number of parameters based on the Bayesian information criterion (BIC) statistic (Schwarz 1978), following the approach of Koch et al. (2021). The BIC statistic consists of a likelihood term ℒ plus an additional term that penalises models with more free parameters, and helps preventing overfitting:

BIC = ln ( m ) k − 2 ln ( L ) . $$ \begin{aligned} \mathrm{BIC} = \ln (m) k - 2 \ln (\mathcal{L} ). \end{aligned} $$(1)

Here, m is the number of fitted data points and k is the number of free parameters. We selected the fit minimising the BIC statistic; we furthermore ignored double component fits in which one component contributes to less than 5% of the total flux. We separated the resulting components into a first component, tracing the galactic rotation, and a second component having a blue or redshift with respect to the disk. In the case where the two components are closer in velocity than 57 km s−1 (corresponding to the MUSE spectral sampling of 1.25 Å at λHα), we picked as the first component the one with the largest amplitude. We found that a double component Gaussian improves the fit only in the starburst region (black rectangle in the left panel of Fig. 2); we therefore present and discuss the resulting kinematic maps in Sect. 6.1, where we study in detail the central kinematics.

3.3. Molecular gas kinematics

In Fig. 3 (bottom row) we show the kinematics of the molecular gas observed in CO(2-1) emission with ALMA. We selected the mom1 and mom2 data obtained with the high confidence ‘strict’ mask (see Leroy et al. 2021a). Despite the difference in resolution, we observe that the molecular gas velocity compares well with the velocity of the ionised gas (centre left panel in Fig. 3). The velocity dispersion is very low (< 15 km s−1) throughout the disk, as expected for the cold gas component, and is only enhanced (∼25–30 km s−1) around the starburst region. We discuss the kinematics of the central starburst region in Sect. 6.1.

3.4. Kinematic fitting

In order to derive the global kinematic parameters of M83, we fitted the velocity maps of stars and gas with the 2D fitting method described in Epinat et al. (2008) and based on the Levenberg-Marquardt non-linear least-square algorithm. In our work, the method has been upgraded using the MocKinG software6 to take into account uncertainties on the line of sight (LoS) velocities, the flux distribution, the linewidth of the profiles and the spatial resolution.

3.4.1. Analysis of the stellar velocity field

We fitted the raw stellar LoS velocity field with several two-dimensional theoretical velocity distributions, and found that the best fit is obtained using two Plummer components (Plummer 1911), one describing the galactic disk and the other tracing the central structure. The Plummer density profile has a finite density core and falls off as r−5 at large radii; this very steep fall-off was essential to fit at best the central component.

The velocity field is described by nine free parameters (see Table 2). They consist of four geometrical parameters – the disk inclination (incl), the position angle of the major axis (PA), the location of the kinematic centre (xc, yc) – plus the heliocentric systemic velocity of the galaxy (vsyst). Additionally, each Plummer component has two free parameters which describe the velocity amplitude and the turnover radius, accounting for four other free parameters. We used an angular sector of inclusion of 67.5 degrees (in the galaxy plane) around the major axis7, and the LoS velocities were weighted according to their angular distance to the major axis by a cosine function, in order to minimise the contamination due to the predominance of radial motions around the minor axis. Regardless of the initial values for the nine parameters, the best fit model converged towards the same solution (model (A) presented in Table 2). We assessed the robustness of the fit by masking low-S/N regions in the velocity map using Kernel filters, testing various several S/N thresholds (see also Appendix A). We found that even when reducing by a factor two the number of fitted spaxels, the output parameters are remarkably stable. The best-fit inclination incl = 20.3 ± 0.1° is in excellent agreement with the one of 21° listed by Gadotti et al. (2019), the position angle PA = 223.1 ± 0.1° is also very close to the one of 227° found by Sheth et al. (2010) and the systemic velocity vsyst = 514.2 km s−1 is only 6 km s−1 higher than the one from LEDA (Paturel et al. 2003; Makarov et al. 2014). The best-fit centre is consistent – within uncertainties – with the Paβ kinematic centre determined by Knapen et al. (2010) (see Fig. A.3). Given the large number of spaxel used in the fit, the statistical uncertainties derived on the fitted parameters are very low; in order to study the robustness of the fit in Appendix A we test three additional models against model (A), adopting different centres from literature values. Using the parameters found with model (A), we computed the rotation curve shown in the top panel of Fig. 4, on top of which we overplot the contribution of the two Plummer components. The agreement between both sides of the rotation curve is very good for r ≳ 50 arcsec (1.1 kpc). Within the first kpc, large discrepancies are observed, with the approaching side rotating faster than the receding one. Regardless of the chosen parameters, it is not possible to make the slopes on both sides of the galaxy coincide within the first 100 pc, indicating that the central stellar component is non-axisymmetric. The discrepancy at the very end of the rotation curve, on the other hand, is due to the decreasing S/N at the edge of the disk.

thumbnail Fig. 4.

Rotation curves for the stellar (top), Hα (centre) and CO (bottom) velocity fields (models (A), (d) and (γ) in Table 2). We remark on the larger radial extent of the CO data on the bottom plot. The contributions of the two best-fit Plummer components – modelling the galactic and circumnuclear disk – are shown in cyan and light green, and the sum of their contributions in purple. The red and blue curves correspond to the receding and approaching side of the velocity map. In Fig. A.1, we additionally show the Hα rotation curve for model (e), where we enforce a more symmetrical curve. In the middle panel, the radial location of features A and B discussed in Sect. 6.1 is indicated.

Table 2.

Parameters derived from the kinematic fitting described in Sect. 3.4.

3.4.2. Analysis of the ionised gas velocity field

We fitted the observed Hα velocity field with the same theoretical velocity distribution as for the stars. We tested several models, whose parameters are also given in Table 2. The best-fit model is model (a). However, the Hα velocity field is so perturbed that the hypothesis of an axisymmetric disk used by the model provides a χ r 2 $ \chi^2_r $ much higher than the one computed for the stellar disk, and does not allow to correctly determine the inclination of the galaxy, notoriously the most difficult parameter to fit (e.g., Epinat et al. 2008). The problem persisted even when masking the central region. To overcome this, we fitted a second model (model b) in which we fixed the inclination to the one of the stellar disk. This change only affected the χ r 2 $ \chi^2_r $ parameter (which increases by < 1%), clearly indicating that the inclination is relatively decoupled from the other parameters. The bottom row of Table 2 shows the average difference between the two sides of the rotation curve. In the literature, the systemic velocity is often determined by fitting a one-dimensional rotation curve. In this case, the amplitude of each side of the curve is matched by adjusting the centre location. On the other hand, when fitting a two-dimensional velocity field as in this work, if the receding and approaching sides display an asymmetric LoS velocity distribution, a difference in rotational velocity (ΔVrot) may appear between the two sides of the rotation curve. We therefore fitted a third model (model c), where we additionally fixed ΔVrot = 0, in order to make the rotation curve more symmetric. This resulted in a moderate increase in χ r 2 $ \chi^2_r $ (5%). Finally, we produced two more models (models d and e) in which we fixed the centre of rotation to the one determined from the stellar disk. This was motivated by the fact that the potential well of the galaxy is dominated by the mass of the stars and of the dark matter halo. In addition, this will facilitate the comparison between the stellar and the gaseous disks. With respect to model (d), in model (e), we also fixed ΔVrot = 0. The only parameter marginally affected was the systemic velocity (vsyst). We show the rotation curve corresponding to model (d) in Fig. 4 (middle panel). In the figure, we also indicate the radial location of the central kinematic features that will be discussed in Sect. 6.1 (A and B, black arrows). In Fig. A.1, we include the Hα velocity rotation curve produced with model (e) as a comparison.

3.4.3. Analysis of the CO gas velocity field

Finally, we fitted the CO data (‘mom1 with prior’ map8, see Leroy et al. 2021b) with the same two Plummer profiles, obtaining the best fit model (α) in Table 2. However, the best fit model exhibits a large inclination incl = 35°, which is interestingly very close to the one derived from the best-fit Hα model (37.9°); this indicates that the warm and cold gas kinematics are similar. For the same reasons mentioned above, in model (β), we fixed the inclination of the CO component to the one of the stellar disk, leaving the other parameters free. None of the fitted parameters were affected by this; only the χ r 2 $ \chi^2_r $ increases by < 2%. For model (γ), we also fixed the centre of rotation of the CO disk to that of the stars. This increased the χ r 2 $ \chi^2_r $ by a further ∼7% with respect to model (β). To facilitate the comparison with the stellar and ionised gas component, in Fig. 4 (bottom panel) we show the rotation curve obtained from model (γ). We remark that since the area spanned by the ALMA data is much more extended than that of MUSE, the CO rotation curve extends ∼50% farther. The first kpc is as difficult to fit as for the Hα component, but we observe a fairly good agreement between the warm and the cold gaseous component. The receding side of both components shows a larger velocity up to the end of the MUSE rotation curve, Beyond this radius, the CO rotation curve indicates that the receding side velocities are lower than the approaching side ones.

4. Extinction traced by the molecular and ionised gas

We estimated the extinction from the Hβ/Hα ratio with PYNEB (Luridiana et al. 2015). We assumed an intrinsic ratio Hα/Hβ = 2.863 (corresponding to case B recombination with Te = 104 K and ne = 100 cm−3, Osterbrock & Ferland 2006) and a Cardelli et al. (1989) extinction law. The Hβ/Hα ratio map was obtained by spatially binning the gas cube to a S/N ≃ 20 in Hβ with the Voronoi tessellation technique. The lines were fitted with a single component Gaussian profile, as described in Sect. 3.2.

The resulting extinction map is shown in Fig. 5 (right panels). Because dust and gas are usually well mixed (e.g., Bohlin et al. 1978), we expect extinction and emission from gas to trace one another, modulo complications due to geometry and the presence of foreground and background sources, as observed for example with CALIFA by Barrera-Ballesteros et al. (2020). By comparing the extinction map with the CO intensity emission traced by ALMA (left panels in Fig. 5), we observe that indeed regions with high extinction correspond to dense molecular gas, with a surface density larger than 100 M⊙ pc−2 (black contours in the top right panel Fig. 5). This is particularly evident along the spiral arms, where the extinction is clearly higher (E(B − V)≳0.5) with respect to interarm region, due to the presence of dense gas. In the bottom right panel of Fig. 5, we show a zoom-in of the central region with CO emission overlaid in black on the extinction map. We observe two high density peaks in the molecular gas density distribution. These regions spatially coincide with the circumnuclear ring and dust inner bar reported in Elmegreen et al. (1998), and also studied by Callanan et al. (2021) in dense gas tracers. The northern peak in CO coincides with the approaching part of the circumnuclear disk observed in Sect. 3.1, and with a peak in the gas extinction map. The southern CO peak (receding part of the circumnuclear disk) shows very low extinction values. Given the inclination of the galaxy, we interpret this as a perspective effect on the vertical separation of the gas, due to the denser molecular gas residing deeper in the disk: the low extinction traced by the ionised gas would then be estimated based on line ratios from gas that is on top of the CO layer.

thumbnail Fig. 5.

Extinction traced by the ALMA and MUSE data. Left panels: CO(2-1) molecular gas emission from ALMA. Right panels: extinction map derived from the MUSE Hβ/Hα ratio. Bottom panels: zoom-in into the central region, where the white lines mark the position of the features discussed in Sect. 6. The contours in the top right panel indicate CO(2-1) emission at 5 (white) and 15 (black) K km s−1 (corresponding to Σmol ∼ 30 and 100 M⊙ pc−2); the black contours in the bottom right panel range from 90 to 800 K km s−1 in steps of ∼90 (corresponding to Σmol ∼ 600−5350 M⊙ pc−2).

5. Properties of the ionised gas

5.1. H II region identification and fraction of DIG

In Fig. 6 (left panel), we show a map of the intensity of the Hα line, obtained by fitting the line in the gas cube spatially binned to a S/N ≃ 20 in Hα. We separated the Hα emission into H II regions and DIG using the Python package ASTRODENDRO9. We remark that within the scope of this work we only require the outer boundaries of star-forming complexes. The detection of individual H II regions will be presented in an upcoming work (Della Bruna et al., in prep.). The data were organised into a hierarchical tree structure of a given depth and minimum leaf size. We set the minimum size of the leaves based on the typical PSF measured at Hα (FWHM = 4.3 pixels ≡ 0 . ″ 86 $ \equiv 0{{\overset{\prime\prime}{.}}}86 $) and we fixed the optional min_delta parameter (minimum difference in flux between two separate structures) to zero. We set the depth of the tree to a surface brightness (SB) threshold corresponding to an H II region ionised by a single low-luminosity O star. We determined this low luminosity threshold using the models of Martins et al. (2005) (Table 1 in their work), which predict an ionising photon flux log Q(H0) = 47.56 photons s−1 for a O9.5 class V star. In the case B approximation10, the Hα luminosity is related to the ionising photon flux as

L ( H α ) [ erg s − 1 ] = α H α eff α B h ν H α · Q ( H 0 ) [ s − 1 ] , $$ \begin{aligned} L(\mathrm{H}\alpha )\ [\mathrm{erg\,s}^{-1}] = \frac{\alpha ^\mathrm{eff}_{\rm H\alpha }}{\alpha _B} h \nu _{\rm H\alpha } \cdot Q(H^0)\ [\mathrm{s}^{-1}], \end{aligned} $$

thumbnail Fig. 6.

Maps of Hα (left panel) and [S II] λ6716,31/Hα (right panel) emission. Both maps have been corrected for extinction and deprojected. The contours indicate the outer limit of star-forming regions and correspond, respectively, to a cut in SBHα = 1.23 × 10−15 erg s−1 cm−2 arcsec−2 and [S II]/Hα = 0.29. In green we indicate the annular sectors for which we compute the radial trends in Fig. 7.

where α H α eff $ \alpha^{\mathrm{eff}}_{\mathrm{H\alpha}} $ is the effective recombination coefficient at Hα and αB is the case B recombination coefficient. Assuming an electron temperature and density Te ∼ 10 000 K and ne = 103 cm−3, this gives (Draine 2011):

L ( H α ) = 1.37 × 10 − 12 Q ( H 0 ) . $$ \begin{aligned} L(\mathrm{H}\alpha ) = 1.37 \times 10^{-12}\ Q(H^0). \end{aligned} $$

We obtained thus L(Hα) = 4.97 × 1035 erg s−1, corresponding to a SB threshold of 1.23 × 10−15 erg s−1 cm−2 arcsec−2, where we assumed an H II region diameter of 10 pc (lower limit estimate). We note that our SB cut is slightly deeper than the one applied in the recent work by Poetrodjojo et al. (2019), which used a cut-off of 1.86 × 10−15 erg s−1 cm−2 arcsec−2 for their H II regions sample. The outer contours of the resulting tree are shown in white in Fig. 6 (left panel). By inverting the H II mask obtained from the dendrogram, we found an (extinction corrected) DIG fraction fDIG = F(Hα)DIG/F(Hα)TOT ∼ 13%. In order to get a first order correction for DIG contamination, we obtained the distribution of the Hα emission over the entire FoV, and estimated its median value via sigma clipping. We then subtracted this value from the H II regions emission. Correcting for the diffuse emission coincident with H II regions results in a negligible increase of ≲0.1%.

We also estimated the fraction of DIG based on a cut in [S II]/Hα, as recently done by Kreckel et al. (2016) and Poetrodjojo et al. (2019). The advantage of this ratio over a simple cut in Hα SB is that it is sensitive to the ionisation state of the gas, as shown in Fig. 6 (right panel): a high ratio (in black) traces regions where the ionising photons have IS+ = 10.4 eV < hν < 13.6 eV = IH+, where I denotes the ionisation potential. An intermediate ratio (orange to purple) indicates gas with IH+ = 13.6 eV < hν < 23.3 eV = IS++ and a low ratio (yellow to orange) marks regions where sulfur is largely doubly ionised (hν > 23.3 eV). As remarked by Kreckel et al. (2016), a cut in [S II]/Hα is better to detect fainter regions, but the generally lower S/N of the [S II] lines can result in irregular boundaries. Figure 6 (right panel) shows a map of the [S II] λ6716,31/Hα ratio, obtained by fitting the emission lines in the gas cube tessellated to S/N ≃ 20 in the [S II] λ6731 line. We observe that the ratio is enhanced in the interarm regions, consistent with DIG observations in the Milky Way and in nearby galaxies (Madsen et al. 2006; Haffner et al. 2009). If the DIG is solely ionised by radiation leaking from H II regions, the enhanced ratios can be explained with the fact that photons escaping from density bounded regions have a harder spectrum between the H I and He I ionisation energies, due to partial absorption, and a softer spectrum at shorter wavelengths, as has been shown for instance by the simulations of a stratified ISM from Wood & Mathis (2004). In their work, Poetrodjojo et al. (2019) applied a cut [S II] λ6716,31/Hα = 0.29, based on typical ratios observed in H II regions and DIG in the MW (Madsen et al. 2006). We adopted the same limit, which results in the contours shown in Fig. 6 (right panel). We observe that the regions contours are similar but somewhat more conservative than the one selected from the Hα map. We recovered a (reddening corrected) fDIG ∼ 20%; also in this case correcting for the diffuse emission coincident with H II regions has a negligible impact (≲0.01%). The recovered value is somewhat lower than the value of 30% estimated by Poetrodjojo et al. (2019) in M83 in the radial range R ≤ 2 Re (using the distance and Re adopted in this work). The discrepancy could be due to the more accurate fitting of the stellar continuum in the MUSE data or to the difference in spatial coverage, depth and spatial and spectral resolution of the two datasets. For the remainder of this work, we adopt the H II region contours selected on the Hα map as a reference.

In Fig. 7, we show the radial flux profile for the H II regions and the DIG. The profile was obtained by radially binning the reddening corrected and de-projected Hα map in annular sectors of width dr = 0.15 Re (indicated in green in Fig. 6, left panel), and matching it with the H II regions contours. The DIG flux was corrected for the diffuse background emission as described above. The first point refers to the emission in the starburst region, 0 ≤ R ≤ 0.15 Re. We see that the luminosity of the H II regions is highest in the centre, due to the starburst activity, and then reflects the configuration of the spiral arms. The SB of the DIG follows a similar trend, although the minimum is offset by 0.15 Re. The resulting fDIG ratio is shown in black in Fig. 7; we see that the DIG contribution to the total luminosity varies between 0.8% and 46%, peaking in the interarm region. The radial trends are in good agreement with what observed by Poetrodjojo et al. (2019, Fig. 6 in their work)11.

thumbnail Fig. 7.

Hα normalised flux as function of radius for the H II regions (red) and DIG (blue), and resulting DIG fraction (black). The data are radially binned in annuli of width dr = 0.15 Re (shown in Fig. 6), and the flux has been normalised by the number of spaxels in each annulus. The H II regions flux has been corrected for diffuse emission.

5.2. BPT diagram analysis

We studied the physical conditions of the ionised gas in ‘BPT’ emission line diagrams (Baldwin et al. 1981; Veilleux & Osterbrock 1987). This set of diagrams offers a powerful tool to interpret ionised gas emission, as the location in the diagram is sensitive to parameters such as the electron density, metallicity, strength and hardness of the radiation field (see e.g., the review by Kewley et al. 2019). We inspected the full set of diagrams, showcasing [N II]/Hα, [S II]/Hα and [O I]/Hα as function of [O III]/Hβ. In each diagram the fluxes were obtained by single component Gaussian fitting of the emission lines in the gas cube tessellated to a S/N = 20 in the weakest line of interest ([O III] λ5007 for the N2- and S2-BPT diagrams, and [O I] λ6300 for the O1 diagram). Typical Voronoi bin sizes in the H II regions vs. the DIG are: 0 . ″ 3 $ 0{{\overset{\prime\prime}{.}}}3 $ vs. 0 . ″ 9 $ 0{{\overset{\prime\prime}{.}}}9 $ in the O III binning and 0 . ″ 7 $ 0{{\overset{\prime\prime}{.}}}7 $ vs. 2″ in the O I binning. We caution that bins located in H II regions are typically below the seeing limit, especially in the blue end of the spectrum (where PSF sizes range between 0.7–0.9 arcsec). However, this should not have a relevant impact within the scope of this work.

The resulting diagrams are shown in Fig. 8: each point corresponds to a Voronoi bin, and is colour-coded as located in an H II region (purple) or in the DIG (green). The black dashed line indicates the location of the ‘extreme starburst’ line by Kewley et al. (2001), corresponding to the upper limit for gas excited purely by SF. Emission above this limit likely originates from shocks or active galactic nuclei (AGN) activity. In Fig. 9 (top panel) we show again the N2-BPT diagram, with each point colour-coded according to its orthogonal distance from the extreme starburst line (Δ N II). Δ N II increases from the bottom left (dark blue, purely SF gas) towards the top right corner (orange-red, purely ionised by shocks or AGN) of the diagram. In the figure, we additionally show the empirical line of Kauffmann et al. (2003, black dotted line), denoting a more stringent limit for gas excited by pure photoionisation. Points located between this line and the extreme starburst line (yellow-green) are likely excited by a mix of SF and shocks or AGN. In the bottom panel of Fig. 9 we show the corresponding ‘2D-BPT’ diagram, where each spatial bin is colour-coded according to Δ N II. We observe that the spiral arms regions stand out as purely SF, the diffuse gas immediately surrounding the regions shows a composite emission, and some of the interarm regions – especially at 0.2–0.4 Re – show a clear signature of shocks.

thumbnail Fig. 8.

BPT emission line diagrams. Each point corresponds to a Voronoi bin, and is colour-coded as located in an H II region (pink) or in the DIG (green). The contours overplotted on the data correspond to iso-proportions of the density, with a probability mass ≤0.25, 0.5, 0.75 and 1. The data are compared to the extreme starburst line from Kewley et al. (2001, black dashed line), the fast shock models from Allen et al. (2008, black grid), and slow shock models from Rich et al. (2011, blue grid).

Finally, we assessed the overall fraction of Hα luminosity originating from SF (regions with Δ N II ≤ 0) and shocks (Δ N II > 0). We observe that, as expected, in H II regions SF accounts for most of the Hα flux (99.8% of the flux originates from regions where photoionisation is the dominant mechanism), whereas bins classified as DIG have a mixed contribution from both photoionisation-dominated regions (accounting for 94.9% of the flux) and shock-dominated regions (accounting for the remaining 5.1%).

We compared our observations with models of fast and slow shocks (black and blue grid in Figs. 8 and 9). The black grid corresponds to the fast shock models of Allen et al. (2008). We use the models that include the photoionising shock precursor, with metallicity Z = 2 Z⊙ and electron density ne = 1 cm−3. We show the model grid spanning b = (0.001–100) μG in magnetic field strength and vs = 290–1000 km s−1 in shock velocity. The blue grid displays the slow shock models described in Farage et al. (2010) and Rich et al. (2011). These models describe shocks driven into the galactic disk by ram pressure originating as a cloud of cool gas (possibly a filament from a merger remnant) falls into the hot ISM halo of a galaxy. The full model grid covers the range 12 + log(O/H) = 7.39–9.39 (∼0.05–5 Z⊙) in metallicity and vs = 100–200 km s−1 in shock velocity; we show however only the super-solar range, Z ≥ 8.69. We observe that the points located beyond the SF limit overlap with shock model grids. In particular in the [N II] diagram, we observe a cloud of points with log[N II]/Hα > 0 that only overlaps with the slow shock models: we investigate more closely these regions in Sect. 6, where we study the kinematics and ionisation state of the starburst region.

thumbnail Fig. 9.

1D (top panel) and 2D (bottom panel) N II-BPT diagram, with each point colour coded according to its distance from the extreme starburst line of Kewley et al. (2001, black dashed line in the top panel). The distance (Δ N II), ranges from dark blue for pure SF to red for pure shock or AGN emission. In the top panel, we also show the empirical demarcation line from Kauffmann et al. (2003, black dotted line), indicating a more stringent limit for photoionised gas. In the bottom panel, the black contours correspond to Δ N II = 0 and white dashed lines indicate the position of the bar and spiral arms.

6. Analysis of the starburst region

We now analyse in more detail the kinematics and ionisation state of the starburst region, indicated by a black rectangle in Fig. 2 (left panel) and shown in more detail in Fig. 10. In the figure, we outline the approximate location of the outer dust ring (of radius ≃9″) and inner bar as determined by Elmegreen et al. (1998) (purple lines) and mark the kinematic features that will be discussed in the rest of this section (black lines).

thumbnail Fig. 10.

Morphology of the M83 starburst region. The background image shows a zoom-in into the three colour stellar composite from MUSE (Fig. 2). The purple dashed lines indicate the schematic location of the outer circumnuclear ring and dust inner bar from Elmegreen et al. (1998). The brightest region of star formation (‘starburst arc’) is visible inside the ring. Black lines indicate the location of the kinematic features discussed in Sect. 6. The scalebar of 20″ corresponds to ≃500 pc at the distance of our target.

6.1. Kinematics of the central starburst region

In Fig. 11 we show a close up view of the kinematics of the central region. In Fig. 12 we additionally inspect the ionised gas kinematics obtained from the 2-component Gaussian analysis. Grey shaded areas in the maps of the second component indicate regions in which the line was best fit by a single component (see description of the fitting method in Sect. 3.2).

thumbnail Fig. 11.

Close up view of Fig. 3 showcasing the kinematics of the starburst region. All velocity maps (left panels) have been corrected for the systemic velocity; no inclination correction has been applied. Black lines mark the kinematic features discussed in Sect. 6.

In the stellar kinematics (top panels in Fig. 11) we observe a fast rotating nuclear component (∼30″ ≃ 700 pc in diameter), already reported by Gadotti et al. (2020). We stress that the alignment of the rotation axis of this inner component (black dashed line) and the stellar bar is purely coincidental, and is expected to significantly vary over secular timescales. We also observe a dip in velocity dispersion (σ < 60 km s−1) at the location of the starburst arc and along the dust lane west of the arc. Similar ‘central dispersion drops’ were reported by Emsellem et al. (2001) and Emsellem (2004) as being possibly due to a dynamically cold stellar component that has formed as a consequence of a bar-driven gas accretion episode.

The ionised gas kinematics are – as already observed in Sect. 3 – more complex. In the velocity map of the single component fit (centre left in in Fig. 11) and in the first component of the double Gaussian fit (top left in Fig. 12), we observe a similar signature of a rotating circumnuclear disk, as already traced by the stellar kinematics. On the east of the nucleus, we observe a stream of gas (labelled as feature C), extending for 50″ (≃1250 pc) and having a velocity difference Δv ≃ +30 km s−1 with respect to the surrounding disk rotation. The stream is surrounded by an extended region (∼1000 × 1600 pc) having an enhanced velocity dispersion ≃60–80 km s−1. This feature was already reported in molecular gas by Lundgren et al. (2004) and in ionised gas by Fathi et al. (2008) as a potential inflowing stream of gas into the central starburst. Also Piqueras López et al. (2012) observe a global velocity gradient in the central region that is pointing to the possible presence of an inflow.

thumbnail Fig. 12.

Results from the 2-component Gaussian fit to the Hα line. First and second row: Hα velocity (left panels) and velocity dispersion (right panels) maps, of the two components. The velocity maps have been corrected for the systemic velocity; no inclination correction has been applied. The velocity dispersion has been corrected for instrumental effects. We draw attention on the different velocity and dispersion scales spanned by the first and second component. Third row: velocity difference between the two components (blue and red indicate, respectively, a second component that is more strongly blueshifted and redshifted with respect to the first one). Black lines mark the kinematic features discussed in Sect. 6. Bottom row: typical spectra of features A and B and C (see labels in the figure). In panels A and B, both the Hα line and the [N II] doublet are clearly doubly peaked, with Δv ∼ 300 km s−1. In panel C, both components are at similar velocity, but the second component has an extremely high σ ∼ 140 km s−1.

Surrounding the nucleus in the map of the second velocity component (second row on the left in Fig. 12), we also observe two conical features on each side of the stellar bar (labelled as features A and B). The two cones are 20″ ≃ 500 pc in size, have a Δv ≃ ±100 km s−1 on top of the disk rotation and a high velocity dispersion σcomp2 ≳ 80 km s−1 (second row on the right in Fig. 12). Cone A appears blueshifted along our LoS, while cone B appears redshifted. In Fig. 5 (bottom panels), we also observed that at the location of cone A there is a peak in the molecular gas emission, whereas the ionised gas traces low extinction. We interpreted this as the fact that the CO emitting gas is located ‘behind’ the ionised gas along our LoS; together with the fact that cone A appears blueshifted, this could indicate that the gas is moving towards us. On the other hand, we do not see a mismatch between CO and E(B − V) along cone B, which might indicate that the ionised gas is moving away from us. Features A and B also clearly stand out in the Hα rotation curve12, as remarked in Fig. 4 (black arrows). On the third and fourth row of Fig. 12, we show the velocity difference between the two Gaussian components, as well as three example spectra, corresponding to features A, B and C. We observe that at the centre of the two cones, the velocity difference between first and second component reaches 200 km s−1.

In the molecular gas (third row in Fig. 11) we observe again the signature of a circumnuclear disk. We also observe some clouds of redshifted molecular gas coincident with feature C. Along these lines of sight, we observe two or three distinct bright peaks of emission, as illustrated in Fig. 13. All but one of these components correspond to gas with a velocity compatible with the galactic rotation, whereas the remaining component is redshifted by ∼100 km s−1 with respect to the disk CO emission. For these complex spectra, we remark that the corresponding mom1 velocity value corresponds to an average of the features. In fitting Gaussian models to some of these clouds, we find that the velocity dispersion of the CO emission at the velocity expected from the circular rotation field is 7–20 km s−1, whereas the redshifted component has a larger velocity dispersion ≃25 km s−1.

thumbnail Fig. 13.

ALMA CO(2-1) spectrum of one of the CO ‘blobs’ along feature C in Fig. 11 (one ALMA spaxel centred on RA = 204.266, Dec = −29.8665). The spectrum has been corrected for the systemic velocity. We observe three separate Gaussian velocity components. The component at vCO ∼ −30 km s−1 is associated with the spiral arm in the galaxy, and the component at −20 km s−1 traces broadening around the first one, likely associated with a bar orbit. We interpret the third component at vLSR ∼ 80 km s−1 as being associated with infalling material.

The peak brightness of the redshifted CO gas blobs is 0.3 to 0.5 K, their velocity dispersions are 15–25 km s−1, and their the clouds are marginally resolved by the ALMA beam (diameters ≳50 pc). Using our adopted conversion factors, the redshifted blobs would have equivalent surface densities of 70–200 M⊙ pc−2 and a total mass of 4 × 106 M⊙. At these surface densities and velocity dispersions, such gas would be typical of what is found in the central regions of barred galaxies13 (e.g., Sun et al. 2020). We remark that features A and B, on the other hand, do not stand out in the CO map. In the remainder of this section, we study features A, B and C in more detail, and in Sect. 7 we discuss their possible origin.

6.2. BPT analysis of the starburst region

In Fig. 14 we show a close up view of the 2D N II-BPT diagram presented in Fig. 9. Blue shaded regions in the map correspond to areas whose emission is overlapping purely with slow shock models (area spanned exclusively by the blue grid in the top panel of Fig. 8). We observe that most of the region surrounding stream C, as well as the far end of cone B (at d ≥ 20″ from the galactic centre) are consistent with pure slow shocks.

thumbnail Fig. 14.

Zoom-in into the 2D N II-BPT diagram from Fig. 9. The grey contours correspond to ΔNII = 0. Blue shaded areas indicate regions whose emission is compatible with slow shock models only (blue grid in Fig. 8).

We also performed a BPT analysis analogous to Sect. 5.2 for the two Gaussian components separately. We fitted the relevant emission lines in the cube binned to a S/N ≃ 20 in the [O III] λ5007 line. Given that for weaker emission lines double components might be harder to disentangle, we first performed a fit to the Hα line and then fixed the kinematic parameters obtained from this fit for all emission lines.

In the top panels of Fig. 15 we show the resulting N2-BPT diagrams. Hereby, points with an uncertainty on either [N II]/Hα or [O III]/Hβ greater than 50% of the ratio are masked, in order to remove bad fits. We observe that both components are tracing both SF and shocks, although the second component extends to more extreme values of [N II]/Hα. This is more clearly visible in the 2D-BPT diagram shown in the bottom panels of Fig. 15. We see that cones A and B are – in both components – consistent with star formation near the stellar bar, and only become shocked at a projected distance d ≥ 20″ from the galactic centre, perhaps tracing an outflow originating from the starburst region that is shocking into the surrounding gas. We discuss this further in Sect. 7.

thumbnail Fig. 15.

N II-BPT diagrams obtained from a double component Gaussian fit to the emission lines. The line and model grids in the top panels are the same as in Fig. 8. Each Voronoi bin is colour coded according to its orthogonal distance from the extreme starburst line (Δ N II), ranging from blue (pure SF) to red (pure shock or AGN emission). Bins with an uncertainty > 50% on either line ratio are excluded from the top plots. In the bottom plots, black lines mark the kinematic features discussed in Sect. 6.

6.3. Shock-sensitive emission line ratios

The top two panels of Fig. 16 show a map of the [O III] λ4959,5007/Hα ratio for the central region (maps of the full FoV can be found in Appendix B). The maps were obtained by fitting the emission lines with a double Gaussian component in the gas cube spatially binned to a S/N ≃ 20 in the [O III] λ5007 line. A high [O III]/Hα ratio is indicative of gas with a high ionisation state, where emission from doubly ionised oxygen (tracing photons with hν ≥ 35.1 eV) is non-negligible with respect to ionised hydrogen (e.g., Veilleux & Osterbrock 1987). We observe that both the region surrounding stream C and the far end of cone B (at d ≥ 20″ from the galactic centre) have a high ratio both in the first ([O III]/Hα≃ 0.4–0.8) and in the second Gaussian component ([O III]/Hα≃ 0.5–1.1). Cone A, on the other hand, only shows a very locally enhanced ratio ([O III]/Hα ≃ 1) in the second component. This could however in part be due to the extremely high extinction at this location, as traced both by the ionised and molecular gas (bottom panel of Fig. 5).

thumbnail Fig. 16.

Ratios of strong emission lines observed with MUSE. Top and centre: map of [O III] λ4959,5007/Hα obtained from a double component Gaussian fit to the emission lines. Bottom panel: map of [O I] λ6300/Hα obtained from a single Gaussian component fit. All maps have been corrected for extinction. The labelled features are discussed in Sect. 6.3. Maps of the full FoV are shown in Fig. B.2.

In the bottom panel of Fig. 16 we show a map of the [O I] λ6300/Hα ratio, obtained from a single Gaussian component; we do not perform a 2-component fit to the [O I] line due to its weak nature. A high ratio of [O I] emission with respect to Hα is indicative of the presence of shocks (Veilleux & Osterbrock 1987). We observe that also in this map both the region surrounding stream C and the far end of cone B have an enhanced ratio ([O I]/Hα≃ 0.5–0.7).

7. Discussion

7.1. Kinematic features in the starburst region

Both models and observations have shown how bars in massive disk galaxies are responsible for physical processes that result in new stellar structures such as nuclear discs or rings and inner bars (Gadotti et al. 2020, and references therein). These processes are driven by bar-induced resonances resulting from the non-axisimmetric potential (e.g., Binney & Tremaine 1987), such as Lindblad resonances14. Dynamical models for the evolution of gas in barred spiral galaxies in 1D (Krumholz & Kruijssen 2015; Krumholz et al. 2017) and 2D (Simkin et al. 1980; Regan & Teuben 2003) are able to model the creation of a nuclear ring within the inner Lindblad resonance (ILR). The main processes involved are described in detail in Krumholz & Kruijssen (2015); see also Renaud et al. (2015); we briefly summarise them here.

In a first phase, the bar exerts torques on the orbiting material. As the gas looses angular momentum, it moves inwards and by energy conservation, the gravitational potential energy is transformed into turbulent energy, resulting in an increase in velocity dispersion. In a second phase, due to the increased velocity dispersion and the mostly flat rotation curve (low shear) within the ILR, acoustic instabilities develop in the gas. This allows for a more efficient transport of angular momentum, and leads to an inflow of gas with high turbulent pressure that is extremely gravitationally stable and has a low SFR. In a third phase, gas starts to build up on the stable ILR orbit, eventually leading in gravitational instabilities that cause fragmentation and collapse, resulting in a circumnuclear ring. The ring can form stars episodically (e.g., in the 1D dynamical models of Krumholz & Kruijssen 2015), or having an initially steady fuelling rate before fragmenting after ∼10 Myr (e.g., in the hydrodynamical simulations from Emsellem et al. 2015).

In Sect. 6.1, we confirmed the presence of a circumnuclear disk, both in the stellar, ionised gas and molecular gas kinematics (Fig. 11). This feature had already been recently observed with MUSE (Gadotti et al. 2020), and has been postulated to coincide with an ILR. Recently, Callanan et al. (2021) analysed high resolution ALMA data mapping the central 500 pc at scales of ∼10 pc, and put forward a model in which the gas revolves around the centre in eccentric orbits. The study proved that the gas in the ring features strong azimuthal variations in velocity dispersion and intensity which are consistent with the expectations from 1D dynamical models (Krumholz & Kruijssen 2015; Krumholz et al. 2017). Furthermore, in their simple model scenario the starburst phase resulting from the disk instability is constrained to be highly localised, both in space and in time, resulting in very efficient stellar feedback.

We furthermore observed that both the ionised and molecular gas are tracing a flow of gas east of the nucleus (feature C in Fig. 11). Lundgren et al. (2004) were the first to report this feature in CO(2-1) and (1-0). Their observations traced – on top of a regular rotating disk – streaming motions along the spiral arms, with the strongest deviation on the NW side of the nucleus. This feature was later confirmed by Fathi et al. (2008) using Fabry-Perot data of the Hα line across the disk. More recently, the Piqueras López et al. (2012) observed evidence of what they intepreted as a gas inflow also in high-resolution NIR IFS data mapping the central ∼200 × 200 pc. In our dataset, we observe that the stream has ΔvHα ≃ +30 km s−1 with respect to the surrounding disk rotation (Fig. 11, centre left) and multiple peaks of emission at different velocities in the molecular gas (Fig. 13), one of which is redsfhited by ∼100 km s−1 with respect to the disk rotation and has an enhanced velocity dispersion (≃25 km s−1). The MUSE data furthermore trace an extended region surrounding the stream (∼1000 × 1600 pc) featuring: (1) a high velocity dispersion (≃80 km s−1, Fig. 11, centre right); (2) [N II]/Hα and [O III]/Hβ ratios that situate the gas clearly above the line separating SF from shocks in a BPT diagram (Figs. 9 and 15), in a region consistent for most part with slow shock models only (Fig. 14); (3) high ratios of [O III] and [O I] with respect to Hα (Fig. 16), indicative of gas with a high ionisation state and of shocks; (4) the presence of bright molecular gas but weak Hα emission (Fig. B.1).

We interpret feature C as the result of two main scenarios. A first possible physical picture is the superposition along the LoS of the disk and an extraplanar layer of DIG. Boettcher et al. (2017) obtained a high-resolution (σHα = 23 km s−1 vs. the MUSE resolution of 120 km s−1) single slit spectrum cutting through feature C. Their data show the presence of two distinct Gaussian components: a narrow component tracing Galactic rotation and a broad component (σ ≃ 95 km s−1) having a velocity lag of ∼25 km s−1 (similar to what we observe in the MUSE data) and high ratios of [N II]/Hα ∼ 1.0. The authors interpret this as the presence of an extraplanar layer of DIG, as broadly observed in other star-forming disks (see e.g., Rossa & Dettmar 2003; Lacerda et al. 2018; Levy et al. 2019; Rautio et al. 2022). This scenario is supported by the multiple peaks observed in the CO spectrum in Fig. 13.

In a second scenario, feature C could be a bar-driven inflow of gas located in the same plane as the disk; in this case the increased velocity dispersion would be tracing excess turbulence as the flow is shocked within the bar. An excess turbulence in the molecular gas could explain the lack of Hα emission despite the bright CO emission along the stream, as the gas would be unable to collapse and form stars.

A final possibility could be a past interaction, possibly with the neighbouring galaxy NGC 5253 (1.8° in projected distance). This option has been taken into account by many authors in order explain the peculiar morphology and kinematics of the central region, where the optical nucleus is offset from the kinematic centre (Thatte et al. 2000; Díaz et al. 2006; Mast et al. 2006; Houghton & Thatte 2008; Rodrigues et al. 2009; Knapen et al. 2010; Piqueras López et al. 2012) and the structure of the H I disk (Miller et al. 2009; Heald et al. 2016). However, given the extremely regular stellar rotation field and the general lack of global scale perturbances, we discard this hypothesis.

In the ionised gas kinematics resulting from the double component Gaussian fit (Fig. 12), we also observed two kinematic features (labelled as cones A and B) where the Hα line is composed of two peaks, separated by a velocity Δv ≤ 200 km s−1 (third and fourth row in Fig. 12) and a high velocity dispersion (up to 200 km s−1, top and centre right in Fig. 12). The two cones appear, respectively, blue- and redshifted along our line of sight (v ≃ ±100 km s−1). Cone B features at its far end (d ≥ 20″ from the galactic centre): (1) [N II]/Hα and [O III]/Hβ ratios above the SF limit in a N2-BPT analysis (Figs. 9 and 15), in a region of the diagram that overlaps mostly with slow shock models (Fig. 14); (2) high ratios of [O III] and [O I] with respect to Hα (Fig. 16), indicative of gas with a high ionisation state and/or shocks; (3) a relatively bright CO emission in a region with little SF (bottom left in Fig. 5 and left panel in Fig. 6). Cone A stands out less clearly in these tracers, perhaps owing to the high extinction at this location (see bottom right panel of Fig. 5). Nonetheless, we observe BPT line ratios indicative of shocks in the second Gaussian component (Fig. 15, for d ≥ 20″) and a locally enhanced [O I]/Hα ratio (central panel in Fig. 16).

So far, no AGN activity has been confirmed in M83, and hard X-ray observations (Yukita et al. 2016) constrain any AGN to be either highly obscured or to have an extremely low luminosity. Given the fact that the origin of cones A and B is offset from the galaxy’s centre, the lack of shock-related emission immediately surrounding the central region and the presence of slow shocks, we favour over the presence of an AGN the hypothesis of a starburst-driven outflow cone shocking into the surrounding ISM. Similar starburst-driven bi-conical outflows have been observed in M82 (Shopbell & Bland-Hawthorn 1998; Leroy et al. 2015b), in which the outflow cones also feature a double component Hα line with Δv = 300 km s−1, low [N II]/Hα ratios compatible with photoionisation and an [O III]/Hα ratio that increases along the outflow, as well as in NGC 253 (Bolatto et al. 2013) and ESO 338-IG04 (Bik et al. 2018). The scenario of a starburst- (rather than AGN-) driven outflow is also in better agreement with the low [S II] λ6716/6731 line ratio observed throughout cones A and B (see e.g., Bik et al. 2018). The ratio is very close to the sensitivity limit ([S II] λ6716/6731 = 1.45, corresponding to ne ≤ 1 cm−3, Draine 2011), and within the measurement errors is largely consistent with ne ≤ 30 cm−3. A starburst-driven outflow would be compatible with the scenario of a short-lived and localised starburst activity by Callanan et al. (2021), resulting in very efficient stellar feedback and, potentially, in nuclear outflows.

An alternative explanation could be the combination of past AGN activity followed by starburst emission. In this scenario, the AGN could have carved a path for the starburst feedback, and could explain the ‘jet-like’ morphology of the outflow and its extent.

Finally, what appear as redshifted and blueshifted cones along our line of sight could simply simply be the result of the superposition of bar-driven orbits, resulting from the non-axisymmetric and time varying potential. However, this effect alone would be difficult to reconcile with the large separation between the two peaks of Hα emission (up to 200 km s−1) at the centre of the cones.

In the near future, the increasing availability of IFS datasets comparable to the present study, covering large portions of local spiral galaxies at intermediate to high resolution, will allow for a better understanding of the signatures imprinted on the ionised gas by in- and outflows, as well as the relative importance of extraplanar gas in low inclination spirals.

7.2. Origin of the DIG

In Sect. 5 we separated the Hα emission in our FoV into compact H II regions and diffuse ionised gas, finding a DIG fraction fDIG ∼ 13% (20%) based on a cut in Hα surface brightness ([S II]/Hα line ratio). This fraction is on the low end of the range with respect to observations in nearby spiral galaxies, where it was constrained to be between 30 and 50% (Ferguson et al. 1996; Hoopes et al. 1996; Zurita et al. 2000; Thilker et al. 2002; Belfiore et al. 2022), and in some cases even up to 60% (Oey et al. 2007). This might be due to the fact that the MUSE data only cover the inner portion of the optical disk (R ≤ 1.1 Re). The wide field spectrograph data by Poetrodjojo et al. (2019) imaging most of the optical disk (R ≤ 2 Re at our assumed distance and Re) revealed indeed that the DIG contribution increases (up to 40%) in the outskirts of the disk. Overall, we observe that despite the enhanced ratios of [N II], [S II] and [O I] with respect to Hα, 94.9% of the Hα luminosity in the DIG is consistent with originating from star formation (BPT diagrams in Figs. 8 and 9). This could support the hypothesis that the diffuse gas is for most part ionised by radiation escaping the SF regions. DIG featuring stronger ratios of [N II]/Hα, [S II]/Hα can instead be partially explained with being extraplanar, as observed by Boettcher et al. (2017). The [N II] and [S II]/Hα ratios are indeed known to increase with distance from the galactic midplane (Madsen et al. 2006; Jones et al. 2017; Levy et al. 2019), due for example to shocks (Rand 1998), turbulent mixing layers (Rand 1998; Collins & Rand 2001) or ionisation by hot, old, low-mass evolved stars (HOLMES, Lacerda et al. 2018; Levy et al. 2019; Weber et al. 2019). We will investigate the origin of the diffuse ionised gas in more detail in a second publication (Della Bruna et al., in prep.).

8. Conclusions

We have presented a large MUSE mosaic covering the central 3.8 kpc × 3.8 kpc of the nearby barred spiral galaxy M83 with a spatial resolution ∼20 pc. We obtained the kinematics of the stars and the ionised gas, and compared them with molecular gas kinematics from ALMA CO(2-1). We observed that the stellar kinematics trace regular rotation along the main SF disk, plus a fast rotating circumnuclear disk of 30″ (≃700 pc) in diameter, likely originating from secular processes driven by the galactic bar.

The gas kinematics are rich in substructures, and the ionised and molecular gas substantially match one another. At large scales, the gas maps rotation along the disk, as well as the fast rotating circumnuclear disk component. On top of the disk rotation, the gas traces a flow east of the nucleus (50″ ≃ 1250 pc in size), that appears redshifted on top of the disk rotation. In the ionised gas, we observe ΔvHα ≃ +30 km s−1 with respect to the surroundings, and in the molecular gas we observe multiple peaks of emission, one of which is redshifted by ΔvCO ≃ +100 km s−1 with respect to the emission tracing galactic rotation. We also observe an enhanced velocity dispersion both in the ionised gas (σHα ≃ 80 km s−1) and molecular gas (σCO ≃ 25 km s−1 for the redshifted peak). The ionised gas showcases an extended region (∼1000 × 1600 pc) surrounding the flow, whose emission lies above the upper limit for star formation in a N2-BPT emission line diagram, and is in large part consistent with models of slow shocks. The extended region also features high ratios of [O III]/Hα (≳0.5) and [O I]/Hα (≳0.5), indicative of a high ionisation state and of the presence shocks. We interpreted this feature as either the superposition of non-collisional flows originating from multiple vertical layers of gas or a bar-driven inflow of shocked gas.

A double component Gaussian fit to the Hα line moreover revealed the presence of two distinct cone-shaped velocity components (20″ ≃ 500 pc in size) on either side of the stellar bar, where the line features two distinct peaks that are up to 200 km s−1 apart. The two cones appear blue- and redshifted along our line of sight, with v = ±100 km s−1 and have a velocity dispersion > 80 km s−1 and up to 200 km s−1. At the far end of both cones, the gas emission lies above the star formation limit in a N2-BPT diagram, in an area of the diagram consistent with slow shock models, and features an enhanced [O III]/Hα ratio (≳0.4 in one or – in the case of cone B – in both of the Gaussian components). One of the two cones also features enhanced [O I]/Hα ratios (≳0.4). We postulate that these two components are tracing a starburst-driven outflow perpendicular to the stellar bar, shocking into the surrounding ISM.

We estimated the gas extinction from the MUSE Hβ/Hα ratio, and found that the regions more strongly affected by extinction (E(B − V) > 0.5) in general have a high density of molecular gas (> 100 M⊙ pc−2). Finally, we separated the ionised gas into H II regions and DIG, using a cut in Hα surface brightness (S II/Hα line ratio). We obtained a DIG fraction fDIG ∼ 13% (20%) in our FoV, and observed that the DIG contribution varies radially between 0.8 and 46%, peaking in the interarm region (at R ∼ 0.3 Re). We inspected the emission of the H II regions and DIG in BPT diagrams, finding that in H II regions 99.8% of L(Hα) originates from photoionisation-dominated regions, whereas the DIG has a mixed contribution from both photoionisation- (94.9%) and shock-dominated regions (5.1%).

1

We determine the outline of spiral arms as the inner radius of the star-forming regions identified in Sect. 5, and the bar position angle based on archival Spitzer S4G data (Sheth et al. 2010; Salo et al. 2015, https://irsa.ipac.caltech.edu/data/SPITZER/S4G/galaxies/NGC5236.html).

2

Project codes 2013.1.01161.S, 2015.1.00121.S and 2016.1.00386.S, PI Sakamoto.

4

FWHMLSF[Å] = 6.266 × 10−8λ2 − 9.824 × 10−4λ + 6.286.

5

Throughout this work, we estimate the S/N of Voronoi maps using a window of 10 Å centred on λobs.

7

The reason for this choice is the fact that the LoS velocities along and close to the minor axis are dominated by radial motions (expansion and contraction along the radius) rather than by galactic rotation.

8

This map uses a high completeness, in order to cover the larger possible number of sightlines.

10

Case B recombination assumes that photons emitted during recombination are immediately reabsorbed, so that recombinations directly to n = 1 are ignored.

11

We note that using the distance and Re adopted in Poetrodjojo et al. (2019), the MUSE data presented in this work only span the range R ≤ 0.8 Re in their Fig. 6.

12

Both on the receding side, due to their position with respect to the angular sector used to compute the curve.

13

We estimate however that these clouds are likely not self-gravitating – if our assumed αCO value is correct – by evaluating the ratio of the kinetic energy (K) to the gravitational binding energy (Ugrav) under a simple spherical model (e.g., Sun et al. 2020; Rosolowsky et al. 2021) along a single line of sight with a radius of half the beam size (R = 25 pc), σCO = 15 km s−1 and a surface density of 100 M⊙ pc−2. The total mass in a synthesised beam is then 2 × 105 M⊙ and K ∼ 30 Ugrav. Alternatively, their opacity and corresponding conversion factors may be lower than for disk clouds, as observed in some galaxy centres (Sandstrom et al. 2013).

14

Lindblad resonances occur for 2(Ω − Ωp) = ± κ, where κ is the frequency of the radial oscillation and Ωp and Ω are, respectively, the angular velocity of the bar and of the stars/gas (neglecting the effect of spiral arms).

Acknowledgments

We thank the anonymous referee for their very detailed and helpful feedback on this manuscript. We thank Jeff Rich and Lisa Kewley for sharing with us the slow shock models presented in Rich et al. (2011). This work is based on observations collected at the European Southern Observatory under ESO programmes 096.B-0057(A), 0101.B-0727(A), 097.B-0899(B), 097.B-0640(A). A.A. acknowledges the support of the Swedish Research Council, Vetenskapsrådet, and the Swedish National Space Agency (SNSA). F.R. acknowledges support from the Knut and Alice Wallenberg Foundation. This research made use of Astropy (http://www.astropy.org), a community-developed core Python package for Astronomy (Astropy Collaboration 2013, 2018).

References

  1. Adamo, A., Kruijssen, J. M. D., Bastian, N., Silva-Villa, E., & Ryon, J. 2015, MNRAS, 452, 246 [NASA ADS] [CrossRef] [Google Scholar]
  2. Allen, M. G., Groves, B. A., Dopita, M. A., Sutherland, R. S., & Kewley, L. J. 2008, ApJS, 178, 20 [Google Scholar]
  3. Astropy Collaboration (Robitaille, T. P., et al.) 2013, A&A, 558, A33 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  4. Astropy Collaboration (Price-Whelan, A. M., et al.) 2018, AJ, 156, 123 [Google Scholar]
  5. Bacon, R., Accardo, M., Adjali, L., et al. 2010, in The MUSE second-generation VLT instrument, SPIE Conf. Ser., 7735, 773508 [NASA ADS] [Google Scholar]
  6. Baldwin, J. A., Phillips, M. M., & Terlevich, R. 1981, PASP, 93, 5 [Google Scholar]
  7. Barrera-Ballesteros, J. K., Utomo, D., Bolatto, A. D., et al. 2020, MNRAS, 492, 2651 [NASA ADS] [CrossRef] [Google Scholar]
  8. Belfiore, F., Santoro, F., Groves, B., et al. 2022, A&A, 659, A26 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  9. Bik, A., Östlin, G., Menacho, V., et al. 2018, A&A, 619, A131 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  10. Binney, J., & Tremaine, S. 1987, Galactic Dynamics (Princeton, N.J: Princeton University Press) [Google Scholar]
  11. Bizyaev, D., Walterbos, R. A. M., Yoachim, P., et al. 2017, ApJ, 839, 87 [NASA ADS] [CrossRef] [Google Scholar]
  12. Blair, W. P., Chandar, R., Dopita, M. A., et al. 2014, ApJ, 788, 55 [NASA ADS] [CrossRef] [Google Scholar]
  13. Boettcher, E., Gallagher, J. S. I., & Zweibel, E. G. 2017, ApJ, 845, 155 [NASA ADS] [CrossRef] [Google Scholar]
  14. Bohlin, R. C., Savage, B. D., & Drake, J. F. 1978, ApJ, 224, 132 [Google Scholar]
  15. Bolatto, A. D., Wolfire, M., & Leroy, A. K. 2013, ARA&A, 51, 207 [CrossRef] [Google Scholar]
  16. Bresolin, F., Kudritzki, R.-P., Urbaneja, M. A., et al. 2016, ApJ, 830, 64 [NASA ADS] [CrossRef] [Google Scholar]
  17. Bryant, J. J., Owers, M. S., Robotham, A. S. G., et al. 2015, MNRAS, 447, 2857 [Google Scholar]
  18. Bundy, K., Bershady, M. A., Law, D. R., et al. 2015, ApJ, 798, 7 [Google Scholar]
  19. Buta, R., & Crocker, D. A. 1993, AJ, 105, 1344 [NASA ADS] [CrossRef] [Google Scholar]
  20. Buta, R. J., Sheth, K., Athanassoula, E., et al. 2015, ApJS, 217, 32 [Google Scholar]
  21. Callanan, D., Longmore, S. N., Kruijssen, J. M. D., et al. 2021, MNRAS, 505, 4310 [NASA ADS] [CrossRef] [Google Scholar]
  22. Calzetti, D., Harris, J., Gallagher, J. S. I, et al. 2004, AJ, 127, 1405 [NASA ADS] [CrossRef] [Google Scholar]
  23. Cappellari, M. 2017, MNRAS, 466, 798 [Google Scholar]
  24. Cappellari, M., & Copin, Y. 2003, MNRAS, 342, 345 [Google Scholar]
  25. Cappellari, M., & Emsellem, E. 2004, PASP, 116, 138 [Google Scholar]
  26. Cardelli, J. A., Clayton, G. C., & Mathis, J. S. 1989, ApJ, 345, 245 [Google Scholar]
  27. Collins, J. A., & Rand, R. J. 2001, ApJ, 551, 57 [NASA ADS] [CrossRef] [Google Scholar]
  28. Comerón, S., Knapen, J. H., Beckman, J. E., et al. 2010, MNRAS, 402, 2462 [Google Scholar]
  29. Croom, S. M., Lawrence, J. S., Bland-Hawthorn, J., et al. 2012, MNRAS, 421, 872 [NASA ADS] [Google Scholar]
  30. Dale, J. E. 2015, New Astron. Rev., 68, 1 [CrossRef] [Google Scholar]
  31. Dale, J. E., Ngoumou, J., Ercolano, B., & Bonnell, I. A. 2014, MNRAS, 442, 694 [Google Scholar]
  32. Della Bruna, L., Adamo, A., Bik, A., et al. 2020, A&A, 635, A134 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  33. Della Bruna, L., Adamo, A., Lee, J. C., et al. 2021, A&A, 650, A103 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  34. den Brok, M., Carollo, C. M., Erroz-Ferrer, S., et al. 2020, MNRAS, 491, 4089 [NASA ADS] [Google Scholar]
  35. den Brok, J. S., Chatzigiannakis, D., Bigiel, F., et al. 2021, MNRAS, 504, 3221 [NASA ADS] [CrossRef] [Google Scholar]
  36. de Vaucouleurs, G., de Vaucouleurs, A., Corwin, H. G. J., et al. 1991, Third Reference Catalogue of Bright Galaxies (New York: Springer) [Google Scholar]
  37. Díaz, R. J., Dottori, H., Aguero, M. P., et al. 2006, ApJ, 652, 1122 [CrossRef] [Google Scholar]
  38. Diehl, S., & Statler, T. S. 2006, MNRAS, 368, 497 [NASA ADS] [CrossRef] [Google Scholar]
  39. Draine, B. T. 2011, Physics of the Interstellar and Intergalactic Medium (Princeton University Press) [Google Scholar]
  40. Drissen, L., Martin, T., Rousseau-Nepton, L., et al. 2019, MNRAS, 485, 3930 [NASA ADS] [CrossRef] [Google Scholar]
  41. Elmegreen, D. M., Chromey, F. R., & Warren, A. R. 1998, AJ, 116, 2834 [NASA ADS] [CrossRef] [Google Scholar]
  42. Emsellem, E. 2004, in Coevolution of Black Holes and Galaxies, ed. L. C. Ho, 11 [Google Scholar]
  43. Emsellem, E., Greusard, D., Combes, F., et al. 2001, A&A, 368, 52 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  44. Emsellem, E., Renaud, F., Bournaud, F., et al. 2015, MNRAS, 446, 2468 [Google Scholar]
  45. Emsellem, E., Schinnerer, E., Santoro, F., et al. 2022, A&A, 659, A191 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  46. Epinat, B., Amram, P., Marcelin, M., et al. 2008, MNRAS, 388, 500 [Google Scholar]
  47. Erroz-Ferrer, S., Carollo, C. M., den Brok, M., et al. 2019, MNRAS, 484, 5009 [Google Scholar]
  48. Farage, C. L., McGregor, P. J., Dopita, M. A., & Bicknell, G. V. 2010, ApJ, 724, 267 [NASA ADS] [CrossRef] [Google Scholar]
  49. Fathi, K., Beckman, J. E., Lundgren, A. A., et al. 2008, ApJ, 675, L17 [NASA ADS] [CrossRef] [Google Scholar]
  50. Ferguson, A. M. N., Wyse, R. F. G., Gallagher, J. S. I., & Hunter, D. A. 1996, AJ, 111, 2265 [NASA ADS] [CrossRef] [Google Scholar]
  51. Gadotti, D. A., Sánchez-Blázquez, P., Falcón-Barroso, J., et al. 2019, MNRAS, 482, 506 [Google Scholar]
  52. Gadotti, D. A., Bittner, A., Falcón-Barroso, J., et al. 2020, A&A, 643, A14 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  53. Gaia Collaboration (Brown, A. G. A., et al.) 2018, A&A, 616, A1 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  54. Girardi, L., Bressan, A., Bertelli, G., & Chiosi, C. 2000, A&AS, 141, 371 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  55. Guérou, A., Krajnović, D., Epinat, B., et al. 2017, A&A, 608, A5 [Google Scholar]
  56. Hadfield, L. J., Crowther, P. A., Schild, H., & Schmutz, W. 2005, A&A, 439, 265 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  57. Haffner, L. M., Dettmar, R. J., Beckman, J. E., et al. 2009, Rev. Mod. Phys., 81, 969 [CrossRef] [Google Scholar]
  58. Heald, G., de Blok, W. J. G., Lucero, D., et al. 2016, MNRAS, 462, 1238 [NASA ADS] [CrossRef] [Google Scholar]
  59. Hoopes, C. G., & Walterbos, R. A. M. 2000, ApJ, 541, 597 [CrossRef] [Google Scholar]
  60. Hoopes, C. G., & Walterbos, R. A. M. 2003, ApJ, 586, 902 [NASA ADS] [CrossRef] [Google Scholar]
  61. Hoopes, C. G., Walterbos, R. A. M., & Greenwalt, B. E. 1996, AJ, 112, 1429 [NASA ADS] [CrossRef] [Google Scholar]
  62. Hopkins, P. F., Narayanan, D., & Murray, N. 2013, MNRAS, 432, 2647 [CrossRef] [Google Scholar]
  63. Hopkins, P. F., Wetzel, A., Kereš, D., et al. 2018, MNRAS, 480, 800 [NASA ADS] [CrossRef] [Google Scholar]
  64. Houghton, R. C. W., & Thatte, N. 2008, MNRAS, 385, 1110 [NASA ADS] [CrossRef] [Google Scholar]
  65. Howard, C. S., Pudritz, R. E., & Harris, W. E. 2017, MNRAS, 470, 3346 [NASA ADS] [CrossRef] [Google Scholar]
  66. Jacobs, B. A., Rizzi, L., Tully, R. B., et al. 2009, AJ, 138, 332 [NASA ADS] [CrossRef] [Google Scholar]
  67. Jones, A., Kauffmann, G., D’Souza, R., et al. 2017, A&A, 599, A141 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  68. Kamann, S. 2018, PampelMuse: Crowded-field 3D Spectroscopy (Astrophysics Source Code Library) [Google Scholar]
  69. Kauffmann, G., Heckman, T. M., Tremonti, C., et al. 2003, MNRAS, 346, 1055 [Google Scholar]
  70. Kewley, L. J., Dopita, M. A., Sutherland, R. S., Heisler, C. A., & Trevena, J. 2001, ApJ, 556, 121 [Google Scholar]
  71. Kewley, L. J., Nicholls, D. C., & Sutherland, R. S. 2019, ARA&A, 57, 511 [Google Scholar]
  72. Knapen, J. H., Sharp, R. G., Ryder, S. D., et al. 2010, MNRAS, 408, 797 [NASA ADS] [CrossRef] [Google Scholar]
  73. Koch, E. W., Rosolowsky, E. W., Leroy, A. K., et al. 2021, MNRAS, 504, 1801 [NASA ADS] [CrossRef] [Google Scholar]
  74. Kreckel, K., Blanc, G. A., Schinnerer, E., et al. 2016, ApJ, 827, 103 [NASA ADS] [CrossRef] [Google Scholar]
  75. Kreckel, K., Ho, I. T., Blanc, G. A., et al. 2019, ApJ, 887, 80 [NASA ADS] [CrossRef] [Google Scholar]
  76. Krumholz, M. R., & Kruijssen, J. M. D. 2015, MNRAS, 453, 739 [Google Scholar]
  77. Krumholz, M. R., Bate, M. R., Arce, H. G., et al. 2014, in Protostars and Planets VI, eds. H. Beuther, R. S. Klessen, C. P. Dullemond, & T. Henning, 243 [Google Scholar]
  78. Krumholz, M. R., Kruijssen, J. M. D., & Crocker, R. M. 2017, MNRAS, 466, 1213 [Google Scholar]
  79. Lacerda, E. A. D., Cid Fernandes, R., Couto, G. S., et al. 2018, MNRAS, 474, 3727 [Google Scholar]
  80. Leroy, A. K., Bolatto, A. D., Ostriker, E. C., et al. 2015a, ApJ, 801, 25 [NASA ADS] [CrossRef] [Google Scholar]
  81. Leroy, A. K., Walter, F., Martini, P., et al. 2015b, ApJ, 814, 83 [Google Scholar]
  82. Leroy, A. K., Schinnerer, E., Hughes, A., et al. 2021a, ApJS, 257, 43 [NASA ADS] [CrossRef] [Google Scholar]
  83. Leroy, A. K., Hughes, A., Liu, D., et al. 2021b, ApJS, 255, 19 [NASA ADS] [CrossRef] [Google Scholar]
  84. Leroy, A. K., Rosolowsky, E., Usero, A., et al. 2022, ApJ, 927, 149 [NASA ADS] [CrossRef] [Google Scholar]
  85. Levy, R. C., Bolatto, A. D., Sánchez, S. F., et al. 2019, ApJ, 882, 84 [NASA ADS] [CrossRef] [Google Scholar]
  86. Lundgren, A. A., Olofsson, H., Wiklind, T., & Rydbeck, G. 2004, A&A, 422, 865 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  87. Luridiana, V., Morisset, C., & Shaw, R. A. 2015, A&A, 573, A42 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  88. Madsen, G. J., Reynolds, R. J., & Haffner, L. M. 2006, ApJ, 652, 401 [CrossRef] [Google Scholar]
  89. Maiolino, R., & Mannucci, F. 2019, A&ARv, 27, 3 [Google Scholar]
  90. Makarov, D., Prugniel, P., Terekhova, N., Courtois, H., & Vauglin, I. 2014, A&A, 570, A13 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  91. Martins, F., Schaerer, D., & Hillier, D. J. 2005, A&A, 436, 1049 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  92. Mast, D., Díaz, R. J., & Agüero, M. P. 2006, AJ, 131, 1394 [NASA ADS] [CrossRef] [Google Scholar]
  93. McLeod, A. F., Dale, J. E., Evans, C. J., et al. 2019, MNRAS, 486, 5263 [Google Scholar]
  94. McLeod, A. F., Kruijssen, J. M. D., Weisz, D. R., et al. 2020, ApJ, 891, 25 [NASA ADS] [CrossRef] [Google Scholar]
  95. McLeod, A. F., Ali, A. A., Chevance, M., et al. 2021, MNRAS, 508, 5425 [NASA ADS] [CrossRef] [Google Scholar]
  96. Miller, E. D., Bregman, J. N., & Wakker, B. P. 2009, ApJ, 692, 470 [NASA ADS] [CrossRef] [Google Scholar]
  97. Morisset, C., Delgado-Inglada, G., Sánchez, S. F., et al. 2016, A&A, 594, A37 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  98. Oey, M. S., Meurer, G. R., Yelda, S., et al. 2007, ApJ, 661, 801 [NASA ADS] [CrossRef] [Google Scholar]
  99. Osterbrock, D. E., & Ferland, G. J. 2006, Astrophysics of Gaseous Nebulae and Active Galactic Nuclei (Sausalito, CA: University Science Books) [Google Scholar]
  100. Paturel, G., Petit, C., Prugniel, P., et al. 2003, A&A, 412, 45 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  101. Piqueras López, J., Davies, R., Colina, L., & Orban de Xivry, G. 2012, ApJ, 752, 47 [CrossRef] [Google Scholar]
  102. Plummer, H. C. 1911, MNRAS, 71, 460 [Google Scholar]
  103. Poetrodjojo, H., D’Agostino, J. J., Groves, B., et al. 2019, MNRAS, 487, 79 [Google Scholar]
  104. Popesso, P., Concas, A., Morselli, L., et al. 2019, MNRAS, 483, 3213 [Google Scholar]
  105. Rand, R. J. 1998, ApJ, 501, 137 [NASA ADS] [CrossRef] [Google Scholar]
  106. Rautio, R. P. V., Watkins, A. E., Comerón, S., et al. 2022, A&A, 659, A153 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  107. Regan, M. W., & Teuben, P. 2003, ApJ, 582, 723 [Google Scholar]
  108. Renaud, F., Bournaud, F., Emsellem, E., et al. 2015, MNRAS, 454, 3299 [NASA ADS] [CrossRef] [Google Scholar]
  109. Rich, J. A., Kewley, L. J., & Dopita, M. A. 2011, ApJ, 734, 87 [NASA ADS] [CrossRef] [Google Scholar]
  110. Rodrigues, I., Dottori, H., Díaz, R. J., Agüero, M. P., & Mast, D. 2009, AJ, 137, 4083 [NASA ADS] [CrossRef] [Google Scholar]
  111. Rosolowsky, E., Hughes, A., Leroy, A. K., et al. 2021, MNRAS, 502, 1218 [NASA ADS] [CrossRef] [Google Scholar]
  112. Rossa, J., & Dettmar, R. J. 2003, A&A, 406, 493 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  113. Rousseau-Nepton, L., Martin, R. P., Robert, C., et al. 2019, MNRAS, 489, 5530 [NASA ADS] [CrossRef] [Google Scholar]
  114. Russell, T. D., White, R. L., Long, K. S., et al. 2020, MNRAS, 495, 479 [NASA ADS] [CrossRef] [Google Scholar]
  115. Salo, H., Laurikainen, E., Laine, J., et al. 2015, ApJS, 219, 4 [Google Scholar]
  116. Sánchez, S. F., Kennicutt, R. C., Gil de Paz, A., et al. 2012, A&A, 538, A8 [Google Scholar]
  117. Sánchez, S. F., Pérez, E., Rosales-Ortega, F. F., et al. 2015, A&A, 574, A47 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  118. Sandstrom, K. M., Leroy, A. K., Walter, F., et al. 2013, ApJ, 777, 5 [Google Scholar]
  119. Scannapieco, C., Wadepuhl, M., Parry, O. H., et al. 2012, MNRAS, 423, 1726 [Google Scholar]
  120. Schaye, J., Crain, R. A., Bower, R. G., et al. 2015, MNRAS, 446, 521 [Google Scholar]
  121. Schwarz, G. 1978, Ann. Stat., 6, 461 [Google Scholar]
  122. Seon, K.-I., & Witt, A. N. 2012, ApJ, 758, 109 [NASA ADS] [CrossRef] [Google Scholar]
  123. Sérsic, J. L., & Pastoriza, M. 1965, PASP, 77, 287 [CrossRef] [Google Scholar]
  124. Sheth, K., Regan, M., Hinz, J. L., et al. 2010, PASP, 122, 1397 [Google Scholar]
  125. Shopbell, P. L., & Bland-Hawthorn, J. 1998, ApJ, 493, 129 [NASA ADS] [CrossRef] [Google Scholar]
  126. Silva-Villa, E., Adamo, A., Bastian, N., Fouesneau, M., & Zackrisson, E. 2014, MNRAS, 440, L116 [NASA ADS] [CrossRef] [Google Scholar]
  127. Simkin, S. M., Su, H. J., & Schwarz, M. P. 1980, ApJ, 237, 404 [NASA ADS] [CrossRef] [Google Scholar]
  128. Sun, J., Leroy, A. K., Schinnerer, E., et al. 2020, ApJ, 901, L8 [NASA ADS] [CrossRef] [Google Scholar]
  129. Thatte, N., Tecza, M., & Genzel, R. 2000, A&A, 364, L47 [NASA ADS] [Google Scholar]
  130. Thilker, D. A., Walterbos, R. A. M., Braun, R., & Hoopes, C. G. 2002, AJ, 124, 3118 [NASA ADS] [CrossRef] [Google Scholar]
  131. Vandenbroucke, B., Wood, K., Girichidis, P., Hill, A. S., & Peters, T. 2018, MNRAS, 476, 4032 [NASA ADS] [CrossRef] [Google Scholar]
  132. Vazdekis, A., Koleva, M., Ricciardelli, E., Röck, B., & Falcón-Barroso, J. 2016, MNRAS, 463, 3409 [Google Scholar]
  133. Veilleux, S., & Osterbrock, D. E. 1987, ApJS, 63, 295 [Google Scholar]
  134. Warner, P. J., Wright, M. C. H., & Baldwin, J. E. 1973, MNRAS, 163, 163 [NASA ADS] [Google Scholar]
  135. Weber, J. A., Pauldrach, A. W. A., & Hoffmann, T. L. 2019, A&A, 622, A115 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  136. Weilbacher, P. M., Streicher, O., Urrutia, T., et al. 2014, in Astronomical Data Analysis Software and Systems XXIII, eds. N. Manset, & P. Forshay, ASP Conf. Ser., 485, 451 [Google Scholar]
  137. Weilbacher, P. M., Monreal-Ibero, A., Verhamme, A., et al. 2018, A&A, 611, A95 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  138. Weilbacher, P. M., Palsa, R., Streicher, O., et al. 2020, A&A, 641, A28 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  139. Williams, B. F., Hillis, T. J., Blair, W. P., et al. 2019, ApJ, 881, 54 [NASA ADS] [CrossRef] [Google Scholar]
  140. Winkler, P. F., Blair, W. P., & Long, K. S. 2017, ApJ, 839, 83 [Google Scholar]
  141. Wood, K., & Mathis, J. S. 2004, MNRAS, 353, 1126 [NASA ADS] [CrossRef] [Google Scholar]
  142. Yukita, M., Hornschemeier, A. E., Lehmer, B. D., et al. 2016, ApJ, 824, 107 [NASA ADS] [CrossRef] [Google Scholar]
  143. Zhang, K., Yan, R., Bundy, K., et al. 2017, MNRAS, 466, 3217 [NASA ADS] [CrossRef] [Google Scholar]
  144. Zurita, A., Rozas, M., & Beckman, J. E. 2000, A&A, 363, 9 [Google Scholar]
  145. Zurita, A., Beckman, J. E., Rozas, M., & Ryder, S. 2002, A&A, 386, 801 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]

Appendix A: Additional material on the kinematics fitting analysis

A.1. Rotation curves and 2D velocity fields

In Fig. A.1 we present the rotation curve for model (e) (same as in Fig. 4 but enforcing ΔVrot = 0). In Fig. A.2 we show two-dimensional (2D) residual velocity fields for the best-fit models, (A) and (a), respectively. 2D–residual fields are used to optimise the free parameter determination of the rotation curves that are computed from 2D–velocity fields. The residual velocity fields (Vres) shown in the figure are obtained after subtracting the best fit models (models (A) and (a) reported in Table 2) from the observed velocity field. Because rotation curves are based on axisymmetric models, specific patterns observed on the residual velocity field would be due to an incorrect parameter determination (Warner et al. 1973). As a consequence, when the free parameters of the rotation curve are optimised, the residual velocity field is typical-pattern-free and the amplitude in the 2D residual velocity is minimum. No artefacts are observed in the residual velocity fields, indicating that the models are optimal. The main residual patterns that remain are due to non-axisymmetirc features (e.g. the large-scale spiral structure and the central component which features an asymmetric X-shape pattern).

thumbnail Fig. A.1.

Same as the middle panel of Fig. 4, but for model (e) in Table 2, which enforces a symmetrical curve.

thumbnail Fig. A.2.

2D residual velocity for model (A) (left) and (d) (right), as described in Table 2.

A.2. Uncertainty estimation

The raw stellar LoS velocity field consists of 1,678,471 velocity measurements. To remove spaxels possibly containing spurious information, the raw stellar field was masked accounting for spaxels with low S/N and/or large errors in the velocity or velocity dispersion determination, which resulted in masking of 0.16% of the spaxels. Despite the nine free parameters, due to the large number of spaxels used in the fit, the associated statistical uncertainties are very low. Therefore, in order to study the robustness of the fit we inspect three additional models for which we set centre locations estimated in the literature. Fig. A.3 shows these locations, all found within a box of ∼3×3 arcsec2, close to the brightest central clumps of M83. For models (B), (C) and (D) we fix the centre to the coordinates computed in two different studies for which uncertainties have been computed. In model (B), we utilise the kinematic centre computed by Fathi et al. (2008) that used the Hα line and 2D-spectroimagery (Fabry-Perot); for model (C) we use the one determined by Knapen et al. (2010) from kinematic Pa β data. Finally, to estimate our own uncertainty on the centre, for model (D) we consider the corner of the uncertainty box from Knapen et al. (2010) (grey box in Fig. A.3) that is the furthest away from the centre determined in model (A). The fits lead to very similar parameters (see table 2). Even if model (A) provides the smallest χ r 2 $ \chi^2_r $, model (B) cannot be discarded just on the relevance of its χ r 2 $ \chi^2_r $ which is only less than 2% higher than those of model (A). Model (C) and (D) lead to χ r 2 $ \chi^2_r $ respectively 5 and 12% larger than model (A). In conclusion we can estimate that our centre uncertainty within a χ r 2 $ \chi^2_r $ variation of 2% is Δα±4 arcsec and Δδ±1 arcsec. We observe that the disagreement in the rotation curves between the approaching and receding sides within the first kpc - and especially within the first 100 pc - becomes increasingly large with models (B), (C) and (D), confirming that model (A) is the best fit.

thumbnail Fig. A.3.

Different centre positions from this work and the literature. Boxes indicate uncertainty determinations from centres as calculated by other authors. The background images show: continuum emission around the Hα line, rescaled in intensity (left panel, arbitrary units) and the observed stellar velocity field (right panel and horizontal colour bar, not corrected for the systemic velocity). The coordinates of the symbols and the box centres are given in Table A.1. Yellow box (Fabry-Perot kinematic centre, Fathi et al. 2008, Knapen et al. 2010, size of the box: 4.0×4.2 arcsec2); grey box (Pa β kinematic centre, Knapen et al. 2010, size of the box: 9.6×8.0 arcsec2); grey ” (corner of the Pa β uncertainty box that is the furthest away from the centre determined using the MUSE stellar velocity map); purple ” (optical nucleus, Díaz et al. 2006); white ” (MUSE stellar kinematic centre from model (A), this work); black ” (CO kinematic centre from model (β), this work).

Table A.1.

Centre positions from this work and the literature shown in Fig. A.3.

Appendix B: Supplementary linemaps

thumbnail Fig. B.1.

Central starburst region observed in ALMA CO(2-1) (left panel) and in MUSE Hα (reddening corrected, right panel). The white lines mark the position of the features discussed in Sect. 6.

thumbnail Fig. B.2.

Same as Fig. 16, but for the full FoV.

Appendix C: MUSE dataset coverage and properties

thumbnail Fig. C.1.

Layout of the MUSE mosaic. Top panel: location of the MUSE mosaic tiles. In green we indicate the data acquired as part of this project (PI Adamo); in blue and purple archival data that we have included in the final mosaic (PI Ibar and Gadotti). Bottom panel: total exposure time of each mosaic tile.

Table C.1.

Central coordinates, exposure time, and seeing of the MUSE mosaic tiles. The seeing has been measured at 7000 Å, as described in Sect. 2.

All Tables

Table 1.

Physical properties of M83.

In the text
Table 2.

Parameters derived from the kinematic fitting described in Sect. 3.4.

In the text
Table A.1.

Centre positions from this work and the literature shown in Fig. A.3.

In the text
Table C.1.

Central coordinates, exposure time, and seeing of the MUSE mosaic tiles. The seeing has been measured at 7000 Å, as described in Sect. 2.

In the text

All Figures

thumbnail Fig. 1.

HST colour-composite image of M83 (Red: F657N, Blue: F814W, Green: F438W). The footprint of the MUSE mosaic is overlaid in white. The position of the bar and the two spiral arms is sketched in black (see Sect. 1). The white contours trace molecular gas mass surface density from ALMA CO(2-1) (contours of 5 and 25 K km s−1, corresponding to a surface gas density ≃30 and 170 M⊙ pc−2).
In the text
thumbnail Fig. 2.

Three-colour composites of the MUSE data. Left: composite of stellar bands (Blue: 4875–4950 Å, Green: 6520–6528 Å, Red: 6750–6810 Å). The green rectangle indicates the central starburst region analysed in Sect. 6. The orange square indicates the position of the Hα cloud highlighted in Fig. 3 (central panels). Right: composite of gas emission (Blue: [O III], Green: [S II] Red: Hα).
In the text
thumbnail Fig. 3.

Kinematics of the stars and gas. Top row: MUSE stellar kinematics. Centre row: MUSE Hα kinematics. The black box in the bottom right corner indicates the position of what we interpret as an extraplanar Hα cloud. Bottom row: molecular gas kinematics from ALMA CO(2-1) (mom1 and mom2 maps). The footprint of the MUSE data is overlaid in black. All velocity maps have been corrected for systemic velocity; no inclination correction has been applied. The velocity dispersion maps refer to the intrinsic dispersion (instrumental effects have been removed). The white dashed lines sketch the position of the bar and spiral arms as shown in Fig. 1 (see Sect. 1). A zoom-in into the central kinematics is shown in Fig. 11.
In the text
thumbnail Fig. 4.

Rotation curves for the stellar (top), Hα (centre) and CO (bottom) velocity fields (models (A), (d) and (γ) in Table 2). We remark on the larger radial extent of the CO data on the bottom plot. The contributions of the two best-fit Plummer components – modelling the galactic and circumnuclear disk – are shown in cyan and light green, and the sum of their contributions in purple. The red and blue curves correspond to the receding and approaching side of the velocity map. In Fig. A.1, we additionally show the Hα rotation curve for model (e), where we enforce a more symmetrical curve. In the middle panel, the radial location of features A and B discussed in Sect. 6.1 is indicated.
In the text
thumbnail Fig. 5.

Extinction traced by the ALMA and MUSE data. Left panels: CO(2-1) molecular gas emission from ALMA. Right panels: extinction map derived from the MUSE Hβ/Hα ratio. Bottom panels: zoom-in into the central region, where the white lines mark the position of the features discussed in Sect. 6. The contours in the top right panel indicate CO(2-1) emission at 5 (white) and 15 (black) K km s−1 (corresponding to Σmol ∼ 30 and 100 M⊙ pc−2); the black contours in the bottom right panel range from 90 to 800 K km s−1 in steps of ∼90 (corresponding to Σmol ∼ 600−5350 M⊙ pc−2).
In the text
thumbnail Fig. 6.

Maps of Hα (left panel) and [S II] λ6716,31/Hα (right panel) emission. Both maps have been corrected for extinction and deprojected. The contours indicate the outer limit of star-forming regions and correspond, respectively, to a cut in SBHα = 1.23 × 10−15 erg s−1 cm−2 arcsec−2 and [S II]/Hα = 0.29. In green we indicate the annular sectors for which we compute the radial trends in Fig. 7.
In the text
thumbnail Fig. 7.

Hα normalised flux as function of radius for the H II regions (red) and DIG (blue), and resulting DIG fraction (black). The data are radially binned in annuli of width dr = 0.15 Re (shown in Fig. 6), and the flux has been normalised by the number of spaxels in each annulus. The H II regions flux has been corrected for diffuse emission.
In the text
thumbnail Fig. 8.

BPT emission line diagrams. Each point corresponds to a Voronoi bin, and is colour-coded as located in an H II region (pink) or in the DIG (green). The contours overplotted on the data correspond to iso-proportions of the density, with a probability mass ≤0.25, 0.5, 0.75 and 1. The data are compared to the extreme starburst line from Kewley et al. (2001, black dashed line), the fast shock models from Allen et al. (2008, black grid), and slow shock models from Rich et al. (2011, blue grid).
In the text
thumbnail Fig. 9.

1D (top panel) and 2D (bottom panel) N II-BPT diagram, with each point colour coded according to its distance from the extreme starburst line of Kewley et al. (2001, black dashed line in the top panel). The distance (Δ N II), ranges from dark blue for pure SF to red for pure shock or AGN emission. In the top panel, we also show the empirical demarcation line from Kauffmann et al. (2003, black dotted line), indicating a more stringent limit for photoionised gas. In the bottom panel, the black contours correspond to Δ N II = 0 and white dashed lines indicate the position of the bar and spiral arms.
In the text
thumbnail Fig. 10.

Morphology of the M83 starburst region. The background image shows a zoom-in into the three colour stellar composite from MUSE (Fig. 2). The purple dashed lines indicate the schematic location of the outer circumnuclear ring and dust inner bar from Elmegreen et al. (1998). The brightest region of star formation (‘starburst arc’) is visible inside the ring. Black lines indicate the location of the kinematic features discussed in Sect. 6. The scalebar of 20″ corresponds to ≃500 pc at the distance of our target.
In the text
thumbnail Fig. 11.

Close up view of Fig. 3 showcasing the kinematics of the starburst region. All velocity maps (left panels) have been corrected for the systemic velocity; no inclination correction has been applied. Black lines mark the kinematic features discussed in Sect. 6.
In the text
thumbnail Fig. 12.

Results from the 2-component Gaussian fit to the Hα line. First and second row: Hα velocity (left panels) and velocity dispersion (right panels) maps, of the two components. The velocity maps have been corrected for the systemic velocity; no inclination correction has been applied. The velocity dispersion has been corrected for instrumental effects. We draw attention on the different velocity and dispersion scales spanned by the first and second component. Third row: velocity difference between the two components (blue and red indicate, respectively, a second component that is more strongly blueshifted and redshifted with respect to the first one). Black lines mark the kinematic features discussed in Sect. 6. Bottom row: typical spectra of features A and B and C (see labels in the figure). In panels A and B, both the Hα line and the [N II] doublet are clearly doubly peaked, with Δv ∼ 300 km s−1. In panel C, both components are at similar velocity, but the second component has an extremely high σ ∼ 140 km s−1.
In the text
thumbnail Fig. 13.

ALMA CO(2-1) spectrum of one of the CO ‘blobs’ along feature C in Fig. 11 (one ALMA spaxel centred on RA = 204.266, Dec = −29.8665). The spectrum has been corrected for the systemic velocity. We observe three separate Gaussian velocity components. The component at vCO ∼ −30 km s−1 is associated with the spiral arm in the galaxy, and the component at −20 km s−1 traces broadening around the first one, likely associated with a bar orbit. We interpret the third component at vLSR ∼ 80 km s−1 as being associated with infalling material.
In the text
thumbnail Fig. 14.

Zoom-in into the 2D N II-BPT diagram from Fig. 9. The grey contours correspond to ΔNII = 0. Blue shaded areas indicate regions whose emission is compatible with slow shock models only (blue grid in Fig. 8).
In the text
thumbnail Fig. 15.

N II-BPT diagrams obtained from a double component Gaussian fit to the emission lines. The line and model grids in the top panels are the same as in Fig. 8. Each Voronoi bin is colour coded according to its orthogonal distance from the extreme starburst line (Δ N II), ranging from blue (pure SF) to red (pure shock or AGN emission). Bins with an uncertainty > 50% on either line ratio are excluded from the top plots. In the bottom plots, black lines mark the kinematic features discussed in Sect. 6.
In the text
thumbnail Fig. 16.

Ratios of strong emission lines observed with MUSE. Top and centre: map of [O III] λ4959,5007/Hα obtained from a double component Gaussian fit to the emission lines. Bottom panel: map of [O I] λ6300/Hα obtained from a single Gaussian component fit. All maps have been corrected for extinction. The labelled features are discussed in Sect. 6.3. Maps of the full FoV are shown in Fig. B.2.
In the text
thumbnail Fig. A.1.

Same as the middle panel of Fig. 4, but for model (e) in Table 2, which enforces a symmetrical curve.
In the text
thumbnail Fig. A.2.

2D residual velocity for model (A) (left) and (d) (right), as described in Table 2.
In the text
thumbnail Fig. A.3.

Different centre positions from this work and the literature. Boxes indicate uncertainty determinations from centres as calculated by other authors. The background images show: continuum emission around the Hα line, rescaled in intensity (left panel, arbitrary units) and the observed stellar velocity field (right panel and horizontal colour bar, not corrected for the systemic velocity). The coordinates of the symbols and the box centres are given in Table A.1. Yellow box (Fabry-Perot kinematic centre, Fathi et al. 2008, Knapen et al. 2010, size of the box: 4.0×4.2 arcsec2); grey box (Pa β kinematic centre, Knapen et al. 2010, size of the box: 9.6×8.0 arcsec2); grey ” (corner of the Pa β uncertainty box that is the furthest away from the centre determined using the MUSE stellar velocity map); purple ” (optical nucleus, Díaz et al. 2006); white ” (MUSE stellar kinematic centre from model (A), this work); black ” (CO kinematic centre from model (β), this work).
In the text
thumbnail Fig. B.1.

Central starburst region observed in ALMA CO(2-1) (left panel) and in MUSE Hα (reddening corrected, right panel). The white lines mark the position of the features discussed in Sect. 6.
In the text
thumbnail Fig. B.2.

Same as Fig. 16, but for the full FoV.
In the text
thumbnail Fig. C.1.

Layout of the MUSE mosaic. Top panel: location of the MUSE mosaic tiles. In green we indicate the data acquired as part of this project (PI Adamo); in blue and purple archival data that we have included in the final mosaic (PI Ibar and Gadotti). Bottom panel: total exposure time of each mosaic tile.
In the text

Current usage metrics show cumulative count of Article Views (full-text article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.

Data correspond to usage on the plateform after 2015. The current usage metrics is available 48-96 hours after online publication and is updated daily on week days.

Initial download of the metrics may take a while.