NOEMA High-fidelity Imaging of the Molecular Gas in and around M82

Nico Krieger; Fabian Walter; Alberto D. Bolatto; Pierre Guillard; Matthew Lehnert; Adam K. Leroy; Jérôme Pety; Kimberly L. Emig; Rebecca C. Levy; Melanie Krips; Hans-Walter Rix; Dragan Salak; Axel Weiss; Sylvain Veilleux

doi:10.3847/2041-8213/ac01e9

1. Introduction

Galaxy-wide outflows driven by star formation are thought to be crucial drivers in galaxy evolution. Stellar feedback caused by intense central star formation activity can launch such outflows, leading to significant fractions of baryons (ionized, atomic, and molecular gas) that escape the main body of the galaxy (e.g., Veilleux et al. 2020). Outflows in starburst galaxies are a multiphase phenomenon and have been observed across the electromagnetic spectrum from X-ray (e.g., Strickland & Heckman 2007), UV (e.g., Hoopes et al. 2005), optical Hα (e.g., Shopbell & Bland-Hawthorn 1998; Westmoquette et al. 2009), IR (e.g., Veilleux et al. 2009), cold dust (Roussel et al. 2010), polycyclic aromatic hydrocarbon emission (e.g., Engelbracht et al. 2006), warm H₂ (e.g., Beirão et al. 2015), and (sub)millimeter and radio emission (e.g., Walter et al. 2002; Bolatto et al. 2013b; Leroy et al. 2015; Martini et al. 2018). While evidence for galactic outflows is manifold, a detailed characterization is restricted to only a few local systems, where the relevant processes can be spatially resolved at high sensitivity. In particular, the physical characterization of the outflowing gas mass is important, as it influences a galaxy's ability to form stars in the future. In this context, the molecular gas phase is particularly relevant because it often carries the dominant mass fraction of all baryons (e.g., Krieger et al. 2019). The fate of the molecular gas in the outflow is itself not clear. Outflowing molecular clumps may be shocked and evaporated by the fast and hot outflowing gas becoming part of the hot phase (Scannapieco & Brüggen 2015; Schneider & Robertson 2017), or they may act as condensation seeds that gain mass, momentum, and velocity by strongly cooling hot gas (Gronke & Oh 2018; Fielding et al. 2020; Abruzzo et al. 2021). Verifying which of these possibilities actually takes place in molecular outflows would strongly impact the interpretation of the observational data, as well as our understanding of the physical processes driving the cool phases of galaxy outflows.

The galaxy M82 is one of the few that show an extended dusty outflow (e.g., Engelbracht et al. 2006; Veilleux et al. 2009; Salak et al. 2013; Beirão et al. 2015; Chisholm & Matsushita 2016). Because of its close proximity (D = 3.5 Mpc) and almost edge-on orientation (inclination i ∼ 80°; e.g., McKeith et al. 1993), M82 offers a unique laboratory to study galactic winds. The CO emission associated with the outflow was first indicated by observations at the Nobeyama 45 m telescope by Nakai et al. (1987) and confirmed in Taylor et al. (2001) based on observations with the Five College Radio Astronomical Observatory single-dish telescope. A more detailed view of the central disk of M82 was achieved by the first interferometric map of M82's CO emission (using the millimeter interferometer of the Owens Valley Radio Observatory; Walter et al. 2002). More recently, a wide-field single-dish CO map (obtained at the IRAM 30 m telescope) demonstrated that the molecular outflow is indeed as extended as observed in other tracers (Hα, X-rays) out to distances of 3 kpc (Leroy et al. 2015). These previous single-dish studies could constrain the overall dynamics and the amount of outflowing molecular gas, but given their effective resolution of ≳300 pc, they could not resolve the actual structure within the outflow.

To characterize the spatial structure of the molecular outflow in M82, we obtained high-resolution observations over a large field of view in M82 with the Northern Extended Millimeter Array (NOEMA). These interferometric observations were complemented with zero-spacing information from the IRAM 30 m telescope. Compared to previous interferometric observations (Walter et al. 2002), we achieve a sensitivity that is three times deeper and a synthesized beam area that is seven times smaller and cover an area on the sky that is three times larger. At a distance to M82 of 3.5 Mpc, 1ʺ corresponds to 17.0 pc (i.e., 1' to 1.02 kpc).

In this paper, we describe the observations and data reduction in Section 2 and present the imaging of the CO(1–0) data in Section 3. We then characterize the small-scale structure of the molecular gas in the outflow, M82's disk, and the surrounding tidal features in Section 4. Finally, our results are summarized in Section 5.

2. Observations

2.1. NOEMA

The NOEMA mosaic covers an area of ∼25 arcmin² over 7 farcm 7 (7.9 kpc) along the major axis and out to ±2 farcm 8 (±2.9 kpc) along the minor axis of M82 . This area covers the regions in the disk and the outflows/streamers that host significant CO emission (I_CO(2−1) > 1.5 K km s⁻¹) as mapped by single-dish observations (Leroy et al. 2015). The mosaic consists of 154 pointings with a hexagonal arrangement with half-width primary beam spacing (21 farcs 5) to achieve an approximately uniform sensitivity across the field of view. While the main focus of the observations is the rotational ground-state transition of carbon monoxide (CO(1–0) at ν_rest = 115.271 GHz), the wide bandwidth correlator at NOEMA also covers other molecular lines, such as ¹³CO(1–0), CS(2–1), or CN(1–0) that are not discussed in this paper.

The NOEMA observations were carried out under project ID W18BY between 2019 March and 2020 January. Out of the 38 observed runs, four had to be dropped entirely due to poor atmospheric conditions, yielding 34 runs with a combined 44 hr of on-source time using an equivalent 10-antenna array. Of the total time, 70% was observed in NOEMA's C configuration and 30% in its D configuration. The flux calibrators were LKHA 101 or MWC 349, and the complex gain calibration was performed on 0836+710 and 0954+658. We tuned the correlator PolyFix to cover the frequency ranges 92.6–100.3 and 107.8–115.6 GHz at 2.03 MHz spectral resolution.

2.2. IRAM 30 m

Short-spacing data were observed with the IRAM 30 m telescope in two runs in April (28/29/30 April) and June/July (30 June, 1/2/3 July) of 2020 for a total of 22 hr of on-source time. The focus and pointing calibrators were 0954+658, 0923+392, and 0836+710. The area of the NOEMA mosaic was covered by a 560'' × 560'' on-the-fly map using two sets of scans along R.A. and decl. The Eight MIxer Receiver was tuned to cover the frequency ranges 94.4–102.4 and 110.5–118.1 GHz, which includes the CO(1–0) line.

2.3. Data Reduction and Imaging

Calibration of the NOEMA interferometric data was done in gildas (version jul20a) using the clic tool. The data were calibrated using the standard pipeline. From the calibrated visibilities, a subset of 800 km s⁻¹ width at 5.0 km s⁻¹ resolution around the CO(1–0) line was extracted for further analysis, as presented in this paper. Continuum emission was fitted in line-free channels away (>∣ ± 250 km s⁻¹∣) from the line center at 210 km s⁻¹ systemic velocity. The obtained continuum fit was then subtracted from each visibility to provide a line-only data set. Self-calibration slightly improved the image fidelity in three of the central pointings that contain bright emission, i.e., most of M82's disk. In these cases, we applied the self-calibration solutions after four iterations.

The 30 m single-dish observations were also reduced in gildas (version jul20a) using the mrtcal and class tools. The mrtcal tool automatically removed the atmospheric contribution and calibrated the intensity scale in units of antenna temperature. In class, we first extracted a frequency window of 500 MHz centered around the rest frequency of 115271.203 MHz from each spectrum, and we converted the intensity scale to the main-beam temperature using using the task beam_efficiency/ruze. We then subtracted a first-order spectral baseline fitted on line-free channels from the line center (same velocity ranges as for the interferometer data) and filtered out all spectra whose baseline noises are larger than three times the standard deviation of the noise distribution. We then resampled the spectra to the same spectral grid as the NOEMA data. Finally, the spectra were gridded through convolution with a Gaussian kernel whose FWHM is one-third of the natural resolution of the IRAM 30 m telescope using the xy_map task. We choose as the projection center the phase center of the NOEMA mosaic and a pixel size of 4'' × 4''. Visual inspection of the obtained position–position–velocity cube reveals well-behaved intensity and noise distribution.

We imaged the NOEMA data first, and then we combined the single-dish and interferometer maps with the casa feather tool. We used the mapping tool in the parallelized version of gildas (version feb20a) to speed up imaging and deconvolution by a factor proportional to the number of cores. In an effort to retain as much of the faint emission away from the central disk as possible, we produced interferometer-only images with natural weighting. As the synthesized beam slightly varies over the field of view, we regularized the deconvolved flux with a Gaussian clean beam whose FWHM was chosen as the largest measured synthesized beam. This gives a uniform spatial resolution of 2 farcs 08 × 1 farcs 65 at a position angle of 51 fdg 2 (0 farcs 4 × 0 farcs 4 pixels). This corresponds to 36.3 × 28.9 pc at the distance of M82. We cleaned the emission (using the Steer clean algorithm, sdi) to an absolute flux level (ares) of 3.0 mJy beam⁻¹ (∼0.6σ) with as many iterations as required (niter=0) inside a clean mask. Corrections for the primary beam response pattern are applied. For feathering, we converted the single-dish map to flux density units and ran feather in casa (version 5.6.1-8). The conversion from flux density to brightness temperature in the final combined map was done using a conversion factor of 26.8 K (Jy beam⁻¹)⁻¹.

The median rms noise per pixel in a 5.0 km s⁻¹ wide channel of the final cube is 138 mK (5.15 mJy beam⁻¹). We achieve a typical noise level of 100–150 mK over the entire map, except for a single pointing toward the south that lacks observing time causing locally increased noise values of ∼250 mK. For all further analysis, we mask out the edge of the mosaic at an offset of ∼16'' from the map edge.

3. Data Presentation

3.1. Data Cube

In Figure 1, we present the CO(1–0) peak brightness map of the NOEMA M82 mosaic. Figure 2 shows the integrated intensity (moment zero), intensity-weighted velocity (moment 1), intensity-weighted velocity dispersion (moment 2), and velocity at peak intensity. At our spatial resolution of ∼30 pc, the observations reveal a high degree of filamentary and clumpy structure of the molecular gas emission in and around the central starburst disk (see discussion in Section 4). The maps also show a complex velocity structure beyond the central rotating disk.

farcs — **Figure 2.** The CO(1–0) moment maps of M82: integrated intensity map (moment zero; top left), intensity-weighted velocity map (moment one; top right), and intensity-weighted velocity dispersion map (moment two; bottom left) of M82. We also show the velocity at peak CO(1–0) intensity in the bottom right panel. The color map of the moment zero map is chosen to increase the visibility of the fainter emission in the molecular streamers and outflow, leading to saturation of the pixels in the central disk (peak ∼3000 K km s⁻¹). All kinematic maps were created from the CO(1–0) data cube blanked at a signal-to-noise ratio (S/N) = 5 threshold. The systemic velocity is ∼210 km s⁻¹. The synthesized beam (∼30 pc, ∼19) is shown in the bottom left corner of each map.
Download figure:
Standard image High-resolution image

The CO(1–0) channel maps are shown in Figure 3. As is already evident from Figure 1, CO(1–0) is detected across the full area mapped by our observations. This includes the central starburst disk, the northern and southern outflow cones, and two tidal arms ("streamers") toward the east and west (see Section 3.2 and Figure 5). For reference, some CO spectra are shown in Figure 4.

**Figure 4.** Single-pixel CO spectra of the selected position. The background CO peak intensity map (Figure 1) is shown to indicate the locations of the spectra. Note that the scale on the y-axis of the spectra changes as a function of position.
Download figure:
Standard image High-resolution image

Both the peak moment map (Figure 1) and the channel maps (Figure 3) confirm the velocity asymmetry of the molecular outflow, but now on smaller spatial scales than seen before (Walter et al. 2002; Leroy et al. 2015). The northern and southern outflows are not symmetric with respect to the galaxy's center (see also, e.g., Shopbell & Bland-Hawthorn 1998); while the northern outflow shows signs of rotation following that of the central disk, the southern outflow appears to be blueshifted and does not feature a clear velocity gradient parallel to the major axis of the disk. This may be caused by contamination of the northern outflow by the foreground disk. The southern outflow breaks up into individual clouds at an approximately constant velocity v ∼ 150 km s⁻¹, interspersed with some clouds at significantly higher (v ∼ 250 km s⁻¹; red) and lower (v ∼ 75 km s⁻¹; dark blue) velocities. This large range of velocities is significantly higher than the dispersion within an individual cloud (<10 km s⁻¹) and suggests that these molecular clouds are associated with the front and back of the ionized/neutral outflow cone (suggested by observations of the ionized outflow; e.g., Lopez et al. 2020). We see a similar range of velocity offsets, albeit less pronounced, in the molecular clouds associated with the northern outflow.

The tidal streamers toward the east and west connect smoothly to the outer disk but show complex internal structure. The bulk of the gas does not display a strong velocity gradient along the streamers out to distances of several kiloparsecs from the center. However, as in the outflow regions, there are some regions where the velocities of individual neighboring clouds (in projection) differ by up to ∼150 km s⁻¹ in the eastern streamer. We attribute these projected velocity differences to distinct velocity components created by the tidal interaction. Toward the southeast, the velocity structure is complex as the eastern streamer (v ∼ 200 km s⁻¹) and southern outflow (v ∼ 150 km s⁻¹) start to overlap in projection.

3.2. Region Definition

For further discussion, we here define five regions that trace the different environments in and around M82: the disk, the northern and southern outflows, as well as the eastern and western tidal streamers (broadly following Walter et al. 2002). This is shown in Figure 5, where we overplot these regions on top of the CO(1–0) peak intensity (molecular gas) and IRAC 4.5 μm (tracing old stars as well as hot dust and near-IR line emission in the outflow) maps.

**Figure 5.** Region masks for the analysis in Section 4 compared to the NOEMA CO(1–0) map (top) and IRAC band 2 image 4.5 μm (bottom). We separate M82 into five regions, focusing on the disk, outflows (north and south), and tidal streamers (east and west), as indicated in the figures. North is up, and east is to the left (cf. Figure 1.)
Download figure:
Standard image High-resolution image

Disk. We define M82's disk region based on its stellar disk, adopting a cut that is brighter than 20 MJy sr^–1 in the IRAC band 2 observations (The SIRTF Nearby Galaxy Survey (SINGS); Kennicutt et al. 2003).

Outflow. We define the outflow regions as two flat symmetric cones (biconical frustum) with a 40° opening angle and a width of 600 pc at the base (excluding the area assigned to the disk), following previous studies of ionized and molecular tracers that have suggested opening angles in the range ∼30° to ∼60° for the outflow (Bland & Tully 1988; Heckman et al. 1990; McKeith et al. 1995; Shopbell & Bland-Hawthorn 1998; Seaquist & Clark 2001; Walter et al. 2002; Engelbracht et al. 2006; Leeuw & Robson 2009; Kaneda et al. 2010; Leroy et al. 2015).

Streamers. The remaining areas toward the east and west of the outflow and disk consist of the disturbed outer disk and two tidal streamers. Following Walter et al. (2002), we denote these as the eastern and western tidal streamers, respectively. These structures are believed to be tidal in nature; the western streamer connects to the neighboring galaxy M81, whereas the eastern streamer points in the opposite direction. These streamers are also seen in H i imaging of significantly larger fields (e.g., de Blok et al. 2018).

3.3. Second Moment Map

We show the second moment map of the CO(1–0) emission in Figure 2 (bottom left). This map was created after blanking individual channels at an S/N ≤ 5 threshold.¹³ We caution the reader that our intrinsic velocity resolution is 5.0 km s⁻¹. Therefore, line widths smaller than this value will have significant uncertainties, as discussed in the deconvolution analysis below.

The central starburst disk displays high velocity dispersion values of typically 20–39 km s⁻¹ (16th–84th percentiles) with a median of σ = 29 km s⁻¹ (mean 30 km s⁻¹) but has peaks exceeding 60 km s⁻¹. The molecular gas that is associated with the northern and southern outflows shows significantly lower values, with median velocity dispersions of 11 (mean 15) and 12 (mean 16) km s⁻¹, respectively. The eastern and western streamers are found at even lower dispersion, with medians of 5 (mean 8) and 8 (mean 11) km s⁻¹, respectively. The transition regions between the disk and the southern outflow, as well as the disk and the eastern streamer, show extended areas of enhanced velocity dispersion. An inspection of the channel maps shows that this is caused by overlapping structures with distinct kinematic components. In the following, we discuss the properties of the individual molecular clouds in more detail.

4. Properties of Molecular Clouds

With 30 pc resolution and 5 km s⁻¹ channels, we identify and characterize structures with sizes similar to individual Galactic giant molecular clouds (see, e.g., Bolatto et al. 2008). In this section, we decompose our cube into individual compact structures, which we call "clouds" for convenience. We measure the size, line width, and luminosity of each cloud and use these to estimate a series of closely related physical properties. Then we measure how the properties of these clouds vary among the different environments in and around M82 and as a function of projected distance from the central starburst. Because we focus on comparative analysis, many of the uncertainties related to calculating physical quantities effectively "divide out" of the analysis. However, we offer a general caution that the absolute values of mass, surface density, and size–line width coefficient implied by our measurements have significant associated systematic uncertainty.

4.1. Cloud Decomposition

We decomposed the CO(1–0) data cube using Fellwalker (Berry 2015), a watershed algorithm, to segment data into a collection of clouds. For each pixel above a background level, Fellwalker follows the steepest gradient to a local maximum. All pixels that end up at the same maximum define one cloud. It has been demonstrated that, compared to other cloud identification algorithms, Fellwalker shows high completeness and accuracy (Li et al. 2020), and we find that the results agree well with what we expect by eye.

We ran Fellwalker using the pyCupid¹⁴ implementation on the NOEMA M82 CO(1–0) data cube. We start the search at the rms noise level (noise), whereas valid peaks must reach twice the noise (minheight). Each cloud's volume must be greater than twice the beam area times the channel width to be considered a distinct cloud (minpix). The dip between two peaks (mindip) must be larger than the noise level for them to be considered disjoint, and the distance must be more than 4 pixels away from each other (maxjump). These seemingly low thresholds of a few times the noise level allow us to quantify clouds in the fainter regions of the outflow and streamers. In these regions, the S/N per individual pixel per channel is typically low (S/N ≲ 3), but the ensemble of dozens of such pixels is significant.

The output of Fellwalker is a cloud assignment array that labels each pixel with the ID of the cloud that it belongs to. In the left panel of Figure 6, we show the 1891 clouds obtained by Fellwalker and project them onto the CO(1–0) peak intensity map, color-coded by their respective systemic velocity. The projection results in many clouds that are situated on top of each other. The superposition of these overlapping clouds then leads to the high dispersion measurements discussed in Section 3.3, but the individual clouds often have narrower line widths. We indicate the velocity dispersions (Section 4.2.2) of individual clouds as colors in the right panel of Figure 6, again overplotted on top of the CO(1–0) peak intensity map.

The cloud catalog is available in machine-readable format (Appendix, Table 2).

4.2. Cloud Property Measurements

4.2.1. Cloud Radii

We start by using the Fellwalker output to calculate the radius for each cloud. We measure the circularized radius, R, from $R=\sqrt{A}/\pi$ , where A is the projected two-dimensional area of the pixels in the Fellwalker-generated cloud assignment. As we discussed in detail in Appendix A of Krieger et al. (2020), this definition of radius differs by a constant factor from alternative definitions when considering cloud ensembles. Note that for individual clouds, the picture is more complex and needs to take the internal cloud structure into account. These effects, however, are below the resolution limit of this study. Therefore, the choice of radius definition does not affect the comparative study of statistical properties presented here. To account for our instrumental resolution, we deconvolve the radius by subtracting the half-width at half-maximum beam size in quadrature. Note that this represents an approximation, because the area of a cloud in the data cube is a function of both the cloud size and S/N; e.g., a cloud with a peak S/N of 50 has a larger footprint than an otherwise identical cloud with a peak S/N of 3. The FWHM beam size used in convolution would be appropriate for a peak S/N = 2 and so should be approximately correct for the faint clouds in the extended streamers and outflows.

We show the resulting distributions of radius in the top panel of Figure 7. Most of our clouds have radii in the range of 40–60 pc, well above the observational limit. Larger clouds are found in the disk (median ∼60 pc, with values up to ∼150 pc) compared to the outflows and streamers (medians 40–45 pc). There appears to be no statistical difference in the radii of an average cloud found in the outflow compared to those in the streamers. We note that by-chance superpositions of clouds are more likely in edge-on disks and may partially contribute to the difference between the disk and outflows or streamers.

4.2.2. Cloud Dispersions

We also measure the line width, or rms velocity dispersion, of each cloud. To do this, we calculate the intensity-weighted second moment for each pixel in the cloud. Then we take the characteristic rms velocity to be the median of this pixel-by-pixel intensity-weighted second moment map over all pixels in the cloud. Similar to the radius, other definitions of line width differ for ensembles of clouds from the one used here only by constant factors (Appendix A in Krieger et al. 2020). We correct the line widths for the effect of instrumental resolution by deconvolving the measured width with the Gaussian equivalent of the instrumental top-hat profile. Formally, the rms velocity dispersion of a top-hat channel profile is $\delta v/\sqrt{2\pi }$ , so our correction is slightly conservative.

In Figure 7 (second from top), we show the distribution of the rms velocity dispersion σ of individual clouds. The median dispersion differs between the disk (∼9.8 km s⁻¹) and the other regions (5.2–6.5 km s⁻¹) but not so much between the outflows (north: median 6.0 km s⁻¹; south: median 5.2 km s⁻¹) and the streamers (east: median 5.2 km s⁻¹; west: median 6.5 km s⁻¹). In all regions, a value of 5–8 km s⁻¹ is the most common velocity dispersion; however, the disk has a wider distribution with prevalent values up to ∼20 km s⁻¹. As for radii, the cloud velocity dispersions in the disk are more likely to be enhanced by chance superpositions than in the outflows or streamers.

Compared to the numbers derived from the second moment map (Section 3.3), the dispersion values calculated from the individual clouds identified by Fellwalker are significantly smaller. This reflects the fact that when calculating second moment measurements based on the full data cube, overlapping clouds will lead to higher line-of-sight velocity dispersions.

We note that our choices for the derivations of cloud radius and velocity dispersion are conceptually different (based on total area versus intensity-weighted dispersion). This is due to the need for robust methods that work across a wide range of sizes, S/N, and environment. However, a statistical comparison of cloud properties is largely insensitive to the details of how cloud properties are defined.

4.2.3. Cloud Luminosities and Masses

We calculate the integrated flux of each cloud by summing the CO flux inside the cloud assignment created by Fellwalker. Combining this flux with the adopted distance to M82, we calculate a CO(1–0) luminosity for each cloud.

We also make a rough estimate of the mass of each cloud by scaling its CO luminosity. This requires assuming a CO-to-H₂ conversion factor. This factor is uncertain, so we adopt two approaches that likely span the range of possible values. First, we use a common conversion factor of α_CO = 1.0 M_⊙ (K km s⁻¹ pc²)⁻¹ for all regions. This is the value for the starburst disk and the bright regions of the outflow derived by Leroy et al. (2015) comparing dust and CO emission. This estimate is consistent with the typical starburst value (α_CO = 0.8 M_⊙ (K km s⁻¹ pc²)⁻¹; Bolatto et al. 2013a). The masses derived by applying this starburst conversion factor to the whole map appear as hatched violins in Figure 7. We might expect the CO-to-H₂ conversion factor to vary across the galaxy. As the conversion factor in the streamers and outflow is unknown, we also calculate a version of the masses using a higher conversion factor of α_CO = 2.5 M_⊙ (K km s⁻¹ pc²)⁻¹ as proposed by Leroy et al. (2015) for these regions. The resulting numbers are shown as dashed, lightly colored violins in Figure 7.

The distribution of luminosity-based cloud masses for each region appears in Figure 7. In the disk, we find median values of 5.8 × 10⁵ M_⊙. In the outflows and streamers, we find cloud masses that are significantly lower (median 5.5–7.6 × 10⁴ M_⊙) than in the disk, even when considering the case where these regions have a higher conversion factor (dashed violins; medians 1.4–1.9 × 10⁵ M_⊙). As was the case with the radius and velocity dispersion, there do not appear to be any significant statistical differences between the mass distribution in the outflows and that in the streamers.

The total gas mass of all clouds identified by Fellwalker is 6.19 × 10⁸ M_⊙, assuming a starburst conversion factor (α_CO = 1.0 M_⊙ (K km s⁻¹ pc²)⁻¹) throughout, and 1.09 × 10⁹ M_⊙, assuming two different conversion factors, as discussed above. This is ∼75% of the total mass obtained from the moment zero map and slightly lower but consistent with previous interferometric imaging (CO(1–0); Walter et al. 2002) and agrees well with single-dish observations (e.g., Salak et al. 2013 for CO(1–0) or Leroy et al. 2015 for CO(2–1)). In those studies, the masses were based on the integrated maps and not individual clouds. The emission that is not recovered in clouds is located in the outflow and streamers and is coincident with structures that are below the cloud limits (see Section 4.1). The bright disk is recovered completely by the cloud decomposition. The recovery fraction in the low-S/N regime could only be increased by lowering the detection thresholds. However, this would in turn lead to an increasing number of noise peaks being detected as clouds, which would affect the cloud statistics, which we aim to avoid.

We calculate the cloud surface density Σ = M/A, which is the surface brightness scaled by a conversion factor, and show the resulting distributions in the bottom panel of Figure 7. Following from cloud mass, surface density is also derived from luminosity and depends on the adopted CO-to-H₂ conversion factor in the same way as discussed above. We further note that surface density is less dependent on observational limits than radius or mass because they partially divide out in the calculation.

We find typical surface densities of 12–140 M_⊙ pc⁻² in the disk (median 42 M_⊙ pc⁻²). In the outflows and streamers, the inferred surface densities fall in the range of 8–25 M_⊙ pc⁻² (medians 11–13 M_⊙ pc⁻²). As for the previously discussed parameters, we again find no differences between the outflows and the streamers.

4.3. Galactocentric Distance Dependence of Cloud Properties

We now explore the dependence of the cloud properties on projected distance from the galaxy center. To do so, we define the projected galactocentric radius as its distance to M82's center at (α, δ) = (09^h55^m52.72^s, $69^\circ 40^{\prime} 45\buildrel{\prime\prime}\over{.} 7$ ). Note that we make no correction for inclination or orientation, and our "radius" is the physical distance along the plane of the sky.

Figure 8 shows the radial trends of each quantity discussed above using sliding medians over bins with a width of 0.1 kpc. As the disk is viewed in edge-on projection and highly asymmetric, the change of disk cloud properties as a function of distance is not meaningful using this radius definition and is not discussed further. We suggest viewing the disk properties primarily as interesting in contrast with the properties in the other regions. However, it is interesting to look at the trends for both the outflows and the streamers.

Radii. We find that the radii for the clouds in the streamers stay approximately constant as a function of radii. The radii for the clouds in the outflow, however—in particular toward the south—appear to decrease with increasing distance.

Dispersions. The radial trends in σ show that, on average, the velocity dispersion decreases with increasing distance. The effect is strongest in the outflow over the range 0.5 kpc < R_center < 1.5 kpc, where the line width appears to drop by almost a factor of 2 from ∼10 to 5 km s⁻¹. At a larger distance, in the streamers, the gradient becomes much shallower, or σ stays approximately constant.

Masses. In the case of the streamers, the individual cloud masses stay approximately constant or show at most a mild decline as a function of distance. The masses of the clouds in both the southern and northern outflows show a clear, rapid decrease with radius. It is unknown how the conversion factor α_CO, as the primary source of uncertainty, varies with galactocentric distance. It could be expected that α_CO increases with distance from the disk value to the outflow/streamer value proposed by Leroy et al. (2015). In this case, the masses in the streamers would, to first order, be constant with distance. The cloud masses in the outflows would still decline with distance but with a shallower slope.

Surface density. The surface density in the outflows drops steeply by a factor of ∼2 within only ΔR_center ∼ 0.5 kpc and stays approximately constant beyond ∼1 kpc. In the streamers, the surface density is only slowly decreasing with distance. Taking the uncertainty in the conversion factor into account would allow for a constant or mildly increasing surface density as a function of distance in both the outflows and the streamers. If the conversion factor were to vary in a clumpy manner instead of a smooth gradient, more complex radial profiles of surface density could be possible.

4.4. Physical State of the Clouds

We now compare the physical state of the individual molecular clouds in the different environments probed by the observations. We do so by first looking at the size–line width relation in Section 4.4.1, followed by assessing the effects of external pressure between the different environments in Section 4.4.2.

4.4.1. Size–Line Width Relation

In a turbulent medium, the measured line width will be larger when considering a larger size scale (e.g., Larson 1981). To compare the different regions independent of this overall scaling, Figure 9 shows the size (radius)—line width (rms velocity dispersion) relation for each of our five regions. In all regions, the line widths of the clouds scale with radius, suggestive of a turbulent medium, and follow a power-law relation (Table 1). As shown in the last panel and Table 1, the relations between the outflow and the streamer regions appear indistinguishable within the uncertainties. Moreover, despite finding a larger variation of R and σ between the disk and other regions, the size–line width relation in the disk is also very similar to that of the outflow and streamers.

**Figure 9.** Size–line width relation of the clouds, comparing between disk, outflow, and tidal streamers. Each panel shows all clouds in that region (black dots), their density (colors), and a power-law fit (black line) as a representative line that follows the distribution. Gray bands show the resolution limits. The bottom right panel compares the lines and allows one to relate the panels with each other.
Download figure:
Standard image High-resolution image

Table 1. Fits to the Size–Line Width Relations $\sigma ={{aR}}^{b}={10}^{{R}_{10}}{R}^{b}$

Region	a	R₁₀	b
Disk	0.20	1.57	0.96
Outflow north	0.29	1.95	0.82
Outflow south	0.18	1.50	0.95
Streamer east	0.36	2.28	0.73
Streamer west	0.32	2.07	0.78

Note. Here R₁₀ is the line width at a representative size scale of 10 pc. The statistical fitting errors are negligible and dominated by the systematic errors. We note that these values should only be used for a relative comparison of the five regions presented in this study.

Download table as: ASCII Typeset image

In other words, the bottom right panel in Figure 9 shows that the clouds we identify in M82 appear to populate a common size–line width relation regardless of region. The primary difference between the regions appears to be what part of the relation is populated. The disk clouds show high R and σ, whereas the clouds in the tidal streamers and the outflows populate the comparatively low R and σ.

4.4.2. Size, Line Width, and Surface Brightness or Surface Density

For a self-gravitating cloud, meaning a marginally bound cloud, a cloud in virial equilibrium, or even a free-falling cloud, we expect a relationship between size, line width, and mass surface density such that σ²/R ∝ Σ (e.g., Keto & Myers 1986; Field et al. 2011). Deviations from this relationship can give insight into the dynamical state of clouds, or they can highlight uncertainties with the physical parameter estimation. Although the mass, and thus surface density, estimates for our clouds are uncertain, we can still gain insight into the dynamic state of the clouds and the origins of the observed line widths. We plot σ²/R versus Σ in Figure 10.

**Figure 10.** Relationship between cloud surface density Σ and the size–line width coefficient σ²/R in different regions of M82. Each panel shows all clouds in this region (black dots), their density (colors), and a power-law fit (black line) as a representative line that follows the distribution. The bottom right panel compares the lines and allows one to relate the panels with each other.
Download figure:
Standard image High-resolution image

The figure shows significant scatter but a similar scaling between Σ and σ²/R for the outflow and streamer regions. The outflow regions tend to show more high surface brightness and high line width points than the streamers, reflecting mostly bright emission near the disk. In all four regions, the slope is slightly steeper than ∼1, but the data appear consistent with a line of slope ∼1. We do not trust our adopted conversion factor enough to know if the clouds are bound, virialized, or in another state. But their combined line width–radius–surface density scaling is consistent with self-gravity playing an important role.

The situation appears different in the central disk, where we observe a shallower slope and a very wide range of surface brightness or surface density. This could reflect that the bright, dense regions are self-gravitating, while the lower surface density regions form part of a more extended, diffuse molecular medium. In this case, the lower surface brightness clouds might be more affected by "external pressure" or simply more dominated by turbulent motions, representing temporary fluctuations in a turbulent medium. Alternatively, it might indicate significant changes in the conversion factor across the disk region. In any case, the differences may be influenced by chance superpositions of clouds that are geometrically more likely in the (close to) edge-on disk.

5. Summary and Conclusion

We present a 154 pointing NOEMA mosaic of the CO(1–0) line emission in and around the nearby starburst galaxy M82. The observations reach a spatial resolution of ∼30 pc, sufficient to spatially resolve the molecular gas in the central starburst disk, the molecular outflow, and the tidal streamers. We obtained a striking CO(1–0) data cube in which we have identified 1891 individual molecular clouds. These data enable us to study the properties of the molecular gas in greatly distinct environments, in particular the following.

Tidal features. Also known as "streamers," they are thought to be due to a recent interaction with M81. The molecular features coincide with the tidal features seen in H i emission (e.g., Yun et al. 1993; de Blok et al. 2018) and are not thought to harbor any significant star formation activity. The expectation is that these regions are "cold," i.e., not highly excited by an external radiation field.

Outflows. Molecular gas is clearly associated with M82's prominent outflow. This outflow is known to harbor hot X-ray gas, with temperatures exceeding 10⁶ K (e.g., Lehnert et al. 1999; Ranalli et al. 2008; Lopez et al. 2020). Even though it is currently unknown how the molecular gas is associated within the hot outflow (mixed gas versus entrained emission), the environment of the gas in the outflow is quite different from that of the cold tidal features.

Disk. The molecular gas in M82's disk has been the focus of many investigations and is not the main target of this study. This environment is known to be very extreme, yielding the high star formation rate density in the central disk that gives rise to the prominent outflow.

These observations thus allow us to compare molecular cloud properties in distinctly different environments. In particular, the molecular outflow and tidal streamers have never been observed at such high spatial and spectral resolution. This allows us to separate clouds that are close in (projected) position but have significant offsets in velocities. In the tidal arms, this indicates overlapping velocity components due to previous galaxy interactions. In the outflows, these offsets are likely due to projection effects of the front and backside of the outflow cone.

There are tendencies for the clouds in the outflow to have smaller sizes (∼0.2 dex kpc⁻¹) and lower velocity dispersions (∼0.3–0.4 dex kpc⁻¹), masses (computed from their luminosity, ∼0.7 dex kpc⁻¹), and surface densities (surface brightness, ∼0.3–0.4 dex kpc⁻¹) at larger distances from the galaxy (Figure 8). The decrease in size, dispersion, and mass (luminosity) is particularly clear toward the more prominent southern outflow. The reduction in the molecular gas surface density (surface brightness) is consistent with that observed in previous work at significantly lower spatial resolution (Leroy et al. 2015). Together, these trends could be indicative of evaporation of the clouds embedded in the hot outflow. To further study the fate of the molecular clouds, multiline high-resolution and sensitivity measurements will be needed alongside simulations that can be directly compared to the observations.

Unlike the clouds in the outflow, we find that the clouds in the streamers stay approximately constant in size and mass, irrespective of distance to M82. They slowly decrease in line width (∼0.1 dex kpc⁻¹) and surface density (∼0.1 dex kpc⁻¹). Also, the size–line width relations for the clouds show indistinguishable behavior between the extragalaxy regions and the galaxy disk itself (Figure 9). The distribution of individual clouds in the σ²/R versus Σ space does not show obvious effects due to changing external pressure in the outflow or streamer clouds, although the disk clouds show a different behavior (Figure 10).

This is a first analysis of this extremely rich and complex data set of M82. In the future, we expect to include ¹³CO information to have an independent estimator of surface density that would enable us to more thoroughly assess the dynamical state of the clouds in the outflow and streamers, as well as multitransition CO spectroscopy. These data will enable comparisons with increasingly realistic simulations of starburst-driven, multiphase outflows in galaxies (e.g., Sparre et al. 2020; Kim et al. 2020).

We thank the referee for a very constructive report that helped to improve the analysis presented in this paper. This work is based on observations carried out under project Nos. w18by and 107-19 with the IRAM NOEMA Interferometer and the IRAM 30 m telescope, respectively. IRAM is supported by INSU/CNRS (France), MPG (Germany), and IGN (Spain).

Facilities: IRAM NOEMA - , IRAM 30 m. -

Software: Gildas,¹⁵ CASA (McMullin et al. 2007), Astropy (Astropy Collaboration et al. 2013, 2018), NumPy (Harris et al. 2020), SciPy (Virtanen et al. 2020), pyCupid,¹⁶ spectral-cube.¹⁷

Appendix: Cloud Catalog

The cloud catalog derived in Section 4 is available in machine-readable format. It includes position, region attribution, distance from M82's center, size, line width, and mass and surface density derived from surface brightness for the 1891 clouds. Table 2 shows an abbreviated overview.

Table 2. Cloud Catalog

Index	R.A.	Decl.	V	Region	R_center	R	σ	M	Σ	μ
	(deg)	(deg)	(km s⁻¹)		(kpc)	(pc)	(km s⁻¹)	(M_⊙)	(M_⊙ pc⁻²)	(K km s⁻¹ pc⁻²)
	(a)	(b)	(c)	(d)	(e)	(f)	(g)	(h)	(i)	(j)
1	148.955	69.677	110.0	Disk	0.32	134.7	23.5	3.23 × 10⁷	568.57	11.49
2	148.971	69.680	250.2	Disk	0.04	122.6	17.8	1.50 × 10⁷	317.91	6.40
3	148.975	69.681	315.3	Disk	0.16	177.5	28.0	5.03 × 10⁷	508.83	10.35
⋯	⋯	⋯	⋯	⋯	⋯	⋯	⋯	⋯	⋯	⋯
1889	148.948	69.688	160.1	Outflow north	0.73	16.2	3.3	9.22 × 10³	11.21	0.114
1890	148.980	69.660	85.0	Outflow south	1.24	19.6	3.5	1.40 × 10⁴	11.55	0.140
1891	148.934	69.669	130.0	Streamer west	1.02	49.2	12.6	1.07 × 10⁵	14.13	0.261

Note. The cloud catalog is published in its entirety in machine-readable format. A portion is shown here for guidance regarding its form and content. (a)–(c) Location of the peak within the cloud. (d) Region according to the definition in Section 3.2. (e) Distance to M82's center at (α, δ) = (09^h55^m5272, $69^\circ 40^{\prime} 45\buildrel{\prime\prime}\over{.} 7$ ). (f) Cloud radius according to Section 4.2.1. (g) Cloud velocity dispersion according to Section 4.2.2. (h) Cloud mass assuming a starburst conversion factor of α_CO = 1.0 M_⊙ (K km s⁻¹ pc²)⁻¹. (i) Cloud surface density assuming a starburst conversion factor of α_CO = 1.0 M_⊙ (K km s⁻¹ pc²)⁻¹. (j) Cloud surface brightness.

Only a portion of this table is shown here to demonstrate its form and content. A machine-readable version of the full table is available.

Download table as: Data Typeset image

Please note that comparisons to other data sets using different methodologies or targets need to be considered carefully. The definitions and parameters (e.g., size, line width, or conversion factor) used here may not match other works, which can lead to differences in normalization or scaling of cloud properties.

NOEMA High-fidelity Imaging of the Molecular Gas in and around M82

Article metrics

Permissions

Author e-mails

Author affiliations

Author notes

ORCID iDs

Dates

Abstract

1. Introduction