H i Kinematics along the Minor Axis of M82

We set up the VLA correlator to observe the full, polarized L-band continuum with 1 MHz channels across the full range ∼1–2 GHz. We also configured spectral windows to target the H i and OH lines at higher spectral resolution. For the H i line, which is the focus of this paper, we used 2048 3.906 kHz channels centered on the H i line at the velocity of M82. For the OH lines, we used a slightly coarser channel width of 7.812 kHz. We calibrated the amplitude, bandpass, and phase with measurements of 1331 + 305 (3C 286; amplitude and bandpass; flux density = 15.0 Jy at 1.4 GHz) and J0921 + 6215 (B-configuration phase; flux density = 1.186 ± 0.003 Jy/beam at 1.55 GHz), and J0841 + 7053 (C- and D-configuration phase; flux density = 3.340 ± 0.018 Jy/beam at 1.57 GHz), respectively, using standard techniques. This observation strategy mirrored that used successfully for THINGS, VLA ANGST, and LITTLE THINGS (Walter et al. 2008; Hunter et al. 2012; Ott et al. 2012). While M82 was extensively observed by the VLA in the 1980s (Yun et al. 1993, 1994), the current observations are more sensitive, have higher velocity resolution, and have much better dynamic range due to upgrades from the early 1990s to late 2000s that improved sensitivity and bandwidth. The latter is particularly important, given M82's substantial continuum emission.

As the continuum emission from M82 is quite strong, we self-calibrated all measurements. The reduction in phase residual following self-calibration resulted in a factor of three to five improvement in peak signal-to-noise for all measurements. After self-calibration, we split out the spectral window containing H i, binning the data to have 8 s integrations. We identified line-free channels from the integrated spectrum of the u–v data. We then fit and subtracted a first-order polynomial from each visibility measurement. We combined all of the continuum-subtracted visibility data for all four observations, in the process binning the data to a final channel width of 5 km s⁻¹.

We imaged the data in stages to account for the substantial absorption in the center of M82. First, we created a high-resolution image with Briggs weighting parameter robust = 0 and no u–v taper. This image had a beam size of ∼5''. The surface brightness sensitivity of this image was too poor to study the outflow in detail, but yielded a good initial model of the inner part of the galaxy. We cleaned in interactive mode until the maximum residual appeared comparable to the noise in the image. We then imaged the data again, beginning from the model output by the previous run, and this time used robust = 0.5 weighting to emphasize surface brightness sensitivity slightly more. Again, we cleaned until the maximum residual was comparable to the noise. Then, starting from the model output by that imaging run, we imaged the data with an 18'' u–v taper. This dramatically improved the surface brightness sensitivity of the data. Because this step began with the previous, higher-resolution model, we achieved a substantial improvement in surface brightness sensitivity without significant, negative impact due to the strong absorption in the center. We cleaned this image until the maximum residual again resembled the noise of the image and smoothed the output image to have a round, 24'' beam (FWHM).

The M81/M82 and NGC 2403 area was observed by moving the telescope in declination and sampling every 3', with an integration time of 1–3 s per sample (the integration time varied with the observation session). Strips of constant declination were spaced by 3', and the telescope was moved in right ascension to form a "basket weave" pattern over the region. The GBT spectrometer was used with a bandwidth of either 12.5 or 50 MHz, depending on the observation session. The combined bandwidth for the final map is 10.5 MHz, and has a velocity range from −890 to 1320 $\mathrm{km}\,{{\rm{s}}}^{-1}$. The typical system temperature for each channel of the dual-polarization receiver was ≈20 K.

The GBT data were reduced in the standard manner using the GBTIDL and AIPS data reduction packages. Spectra were smoothed to a channel spacing of 24.4 kHz, corresponding to a velocity resolution of $5.2\,\mathrm{km}\,{{\rm{s}}}^{-1}$. A reference spectrum for each of the nine observation sessions was made from an observation of an emission-free region, usually from the edges of the maps. The reference spectrum was used to perform a (signal-reference)/reference calibration of each pixel. These calibrated spectra were scaled by the system temperature and corrected for atmospheric opacity and GBT efficiency. We adopted the GBT efficiency from equation one of Langston & Turner (2007) for a zenith atmospheric opacity of τ₀ = 0.009. The frequency range observed was relatively free of radio frequency interference (RFI) and less than 0.3% of all spectra were adversely affected. The spectra exhibiting RFI were identified by tabulating the root-mean-square (rms) noise level in channels free of neutral hydrogen emission. Spectra that showed unusually high noise across many channels were flagged and removed. The observations were then gridded using the AIPS task SDIMG, which also averages polarizations. After amplitude calibration and gridding, a first-order polynomial was fit to line-free regions of the spectra and subtracted from the gridded spectra using the AIPS task IMLIN. The effective angular resolution, determined from maps of 3C286, is 915 ± 005. To convert to units of flux density, we observed the calibration source 3C286. The calibration from Kelvin to Janskys was derived by mapping 3C286 in the same way that the HI maps were produced, and the scale factor from K/Beam to Jy/Beam is 0.43 ± 0.03. Due to the patchwork nature of the observations, the rms noise varies considerably across the data cube, ranging between 6 and 30 mJy/beam. The average rms noise in the final data cube is 20 mJy per 24.4 kHz channel.

The VLA-only cube still exhibits "bowling" and other large-scale artifacts, in part due to the absence of short and zero spacing data. To deal with this, we combined the VLA cube with the GBT cube. First, we extracted a subcube centered on M82 from the larger M81 GBT survey. We then masked channels dominated by Galactic emission, carried out one additional round of first-order baseline fitting, applied the primary beam taper of the VLA data to the GBT data, and combined the VLA and GBT data with the CASA task feather. This combination removed or suppressed most of the large-scale artifacts that were present in the VLA-only cube. After the combination, we corrected the cube for the primary beam response of the VLA. The final cube has a channel width of $5\,\mathrm{km}\,{{\rm{s}}}^{-1}$, beam size of 24'' (FWHM), and covers the primary beam of the VLA out to its half-power point (≈32'). Before the primary beam correction, the cube had an rms noise of ≈0.4 mJy beam⁻¹ at this resolution and channel width, equivalent to ≈4 K noise for ≈1050 K Jy⁻¹.

We converted the cube to have units of Kelvin, and also made a version of our H i cube with individual channels in units of column density. To do this we used the 5 km s⁻¹ width of each channel and assumed optically thin H i emission so that N ≈ 1.823 × 10¹⁸ cm⁻² I_{H i} with I_{H i} in K $\mathrm{km}\,{{\rm{s}}}^{-1}$. To highlight only bright, positive emission, we also created a mask that included all regions of the cube where the intensity exceeded a signal-to-noise ratio of five over two successive channels. This corresponds to a column density limit of ≈7 × 10¹⁹ cm² across a 10 $\mathrm{km}\,{{\rm{s}}}^{-1}$ velocity range. We used this mask to create integrated intensity and intensity-weighted mean velocity maps. Note that most of our analysis focuses on direct analysis of the cubes.

The integrated H i intensity map and velocity field are shown in Figure 1, and channel maps of H i intensity in select velocity channels are shown in Figure 2. The H i intensity map shows substantial H i extending along the major axis for over ±5  kpc. The H i emission to the left (east) shows a significant spur that extends up (north) from the disk starting at about 4–5 kpc. This is the tidal stream seen in previous H i observations (Cottrell 1977; van der Hulst 1979; Appleton et al. 1981; Yun et al. 1993, 1994). This feature is labeled "NE Spur" in Figures 1 and 2. Our data also show some of the extensive pair of streams that extend west of the major axis, and then toward the south (Yun et al. 1993). This feature is labeled "SW Spur" in Figures 1 and 2. The H i emission along the minor axis extends approximately 5 kpc toward the north and 10 kpc toward the south. M82 is to the north of M81, so the southern part of M82 is closer to M81.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The H i velocity field shows some evidence of rotation, similar to the molecular gas (Young & Scoville 1984; Yun et al. 1993; Leroy et al. 2015), although not as regular as the K-band stellar rotation curve measured by Greco et al. (2012). The H i rotation curve exhibits several drops in velocity (Sofue et al. 1992), which are most likely due to the influence of the bar (Telesco et al. 1991) and the tidal interaction (Yun et al. 1993). The velocity field along the minor axis has a steep gradient within about 1 kpc, where it decreases rapidly by over 100 $\mathrm{km}\,{{\rm{s}}}^{-1}$ toward the south and increases rapidly by over 100 $\mathrm{km}\,{{\rm{s}}}^{-1}$ toward the north. This spatial extent is where the start of the outflow is observed at other wavelengths (McKeith et al. 1995; Westmoquette et al. 2009). The amplitude of the velocity field then decreases toward the systemic velocity of M82. This is particularly prominent toward the south, where continuous, albeit filamentary, H i is detected out to 10 kpc from the midplane. At this point the velocity has decreased to approximately half the value at 1 kpc. Figure 3 shows a position–velocity diagram along the minor axis within a synthesized slit of width ±0.5 kpc relative to the center of M82. The boundaries of this slit are shown in Figure 1.

Zoom In Zoom Out Reset image size

The most remarkable feature of the minor axis H i kinematics is the decrease of the typical velocity amplitude from about 100 to 50 $\mathrm{km}\,{{\rm{s}}}^{-1}$ over the projected distance range from about 1.5 kpc to 10 kpc (Figure 3). This decline is very obvious along the minor axis toward the south. The kinematics to the north are broadly consistent with this same velocity profile, although the H i intensity is substantially less and the gas is only detected to about 5 kpc in projection.

The velocity amplitude decrease from 1.5 to 10 kpc is not consistent with H i clouds that have been accelerated by the wind and then launched on ballistic trajectories, and consequently demonstrate that an additional force is required. This is because the inferred outflow velocities of these clouds is on the order of 500–600 $\mathrm{km}\,{{\rm{s}}}^{-1}$ (for 80° disk inclination) at a projected height of 1.5 kpc and the mass of M82 is insufficient to slow the clouds by 50% from 1.5 kpc to 10 kpc. This point is illustrated by the dashed lines in Figure 3. The lines show the projected velocity of a test particle calculated with galpy and the mass model for M82 developed in Section 3. The middle line shows the trajectory for a test particle on the symmetry axis of the bicone (at an angle of 10° relative to the plane of the sky) and the other lines show test particles at ±8° on the nominal edges of a cone with a 16° opening angle. The H i kinematics clearly decrease much faster than expected from gravitational forces alone. They are also inconsistent with cooling of the H i gas from the hot wind. H i produced from cooling in a hot wind should have the velocity of the hot wind, which Strickland & Heckman (2009) estimate to be as high as ∼2000 $\mathrm{km}\,{{\rm{s}}}^{-1}$. The H i kinematics also demonstrate that the cold H i clouds do not experience continual acceleration, nor even free expansion into the halo of M82.

One possible explanation for this decrease in projected velocity is that the inclination of the bicone is substantially different from the stellar disk, as this would decrease the H i outflow velocity. We find that the H i velocities do diminish by the required amount from 1.5–5 kpc if the symmetry axis of the cone is tilted toward us by 30° relative to the plane of the sky. However, in this case the H i velocities continue to decrease to the systemic velocity, and do not reproduce the relatively flat velocity profile seen from 5–10 kpc. Other possible explanations include that the H i experiences some additional drag force, the H i does not represent a continuous stream of material, but instead is forming in situ (although not from the hot wind), and that the H i kinematics have been substantially affected by tidal interactions.

If the H i consists of clouds flowing outward along the edge of the bicone, they could experience a drag force if they encounter ambient material outside of the bicone. This material could be halo gas of M82 that is not part of the wind, or tidal material that is approximately stationary relative to the wind. We model the drag force with the drag equation ${f}_{\mathrm{drag}}\,={C}_{d}{\rho }_{\mathrm{amb}}{v}^{2}A$, where C_d is the drag coefficient, ρ_amb is the density of the ambient medium, v is the velocity of the H i relative to this medium, and A is the cross section of a typical H i cloud. We then make the assumption that C_dρ_amb A ∝ r⁻². This would be true if C_dA were constant and ρ_amb ∝ r⁻². The coefficient C_d is approximately constant for spheres (of constant cross section A) in a supersonic flow, while ρ_amb ∝ r⁻² is characteristic of mass outflow at a constant rate at constant velocity. Based on measurements by Leroy et al. (2015), we expect the density profile is somewhat steeper with ${\rho }_{\mathrm{amb}}\propto {r}^{-3\to -4}$ for the dust and molecular gas.⁹ To maintain our assumption in this case would require the drag coefficient to increase with r, which could happen if the Reynolds number decreased with radius. We find that the addition of a drag force with C_dρ_amb A ∝ r⁻² provides a reasonable match to the data. This model is shown in Figure 6.

Zoom In Zoom Out Reset image size

If the velocity change of the H i gas is due to a drag force imposed by the ambient medium, the non-gravitational deceleration of the test particles that best match the data constrain the product of C_dρ_ambA_cl. If we assume these are spherical clouds of H i gas with uniform density ρ_cl and radius R_cl, the non-gravitational deceleration is

$$\begin{eqnarray}&&\displaystyle \frac{{dv}}{{dt}}=\displaystyle \frac{3}{4}{C}_{d}\left(\displaystyle \frac{{\rho }_{\mathrm{amb}}}{{\rho }_{\mathrm{cl}}}\right)\displaystyle \frac{{v}_{\mathrm{cl}}^{2}}{{R}_{\mathrm{cl}}}.\end{eqnarray} \tag{ 1 }$$

At 2 kpc from the midplane, the deprojected H i velocity is about 500 $\mathrm{km}\,{{\rm{s}}}^{-1}$ and the best-fit model has a non-gravitational deceleration of $1.8\times {10}^{-7}\,\mathrm{cm}\,{{\rm{s}}}^{-2}$. For R_cl = 10 pc, the density ratio is ∼0.0044 (0.5/C_d). If the cloud particle density is 10 cm⁻³, the ambient medium at 2 kpc from the midplane must be relatively dense with 0.044 cm⁻³. This is substantially denser than the hot wind, although the wind has much larger velocities so it is unlikely to be the same medium. If drag forces decelerate the H i gas, then the ambient medium responsible for the deceleration must be distinct from the hot wind, and at higher density than the typical hot, coronal gas of galactic halos. This could be due to material stripped out of the disk by the tidal interaction with M81. The tidal interaction may also be the main origin of the decrease in the H i velocity, rather than a drag force due to the ambient medium.

The combination ρ_clR_cl in Equation (1) is the column density of the clouds in the outflow. If the majority of the material in the cloud is H i, then this quantity can be measured with H i emission maps or UV/visible absorption line studies with sufficiently high resolution and sensitivity. Such a measurement, combined with the acceleration and the velocity of the clouds, could then be used to solve for the ambient density required to produce the observed deceleration. The angular size of a ∼10–20 pc cloud at the distance of M82 is 05–1''. This is a challenge for current H i facilities, but may be feasible with UV/visible absorption line studies toward background sources.

Similar to Leroy et al. (2015), we have used the mass distribution and velocity information to derive the mass flux associated with the H i emission. Under the assumption that all of the H i is associated with an outflow, Figure 7 shows the mass flux computed from the product of the mass surface density (Figure 5), the velocity from our outflow model, and the ±0.5 kpc width of the central region of our position–velocity diagram. The main differences from the calculation and figure shown in Leroy et al. (2015) are that we used a smaller region width (1 kpc, rather than 3 kpc) and we used the velocity profile, rather than a fixed outflow speed of 450 $\mathrm{km}\,{{\rm{s}}}^{-1}$. While the mass flux we have derived near the midplane is lower than in Leroy et al. (2015), largely because of the smaller region width, the mass flux is detected out to about twice the distance from the midplane, although it drops to a fraction of a solar mass per year beyond 4 kpc. This is also approximately the point where the mass flux is comparable to the outflow rate computed by Strickland & Heckman (2009) for the hot gas. This H i outflow rate is about an order of magnitude less than the star formation rate.

Zoom In Zoom Out Reset image size

The H i velocity decrease along the minor axis may be so substantial that the H i will not escape the halo. At a projected height of 10 kpc, the line-of-sight velocity has dropped to 50 $\mathrm{km}\,{{\rm{s}}}^{-1}$ relative to the systemic velocity, which corresponds to an inclination-corrected outflow velocity of about 300 $\mathrm{km}\,{{\rm{s}}}^{-1}$. This is sufficiently comparable to the escape velocity of our mass model that the ultimate fate of the gas will depend on assumptions in our model for the drag force on the cold gas, the impact of tidal forces in the M81 group, and the halo model for M82. Within the context of our model, the H i would finally stall at 70–80 kpc after about 500 Myr. The H i gas could consequently form a cold "fountain" as proposed by Leroy et al. (2015), rather than substantially add to the intergalactic medium. We searched in our data for any evidence of material that is falling back, but do not see any. This included a careful inspection of the non-masked data, which has somewhat greater sensitivity, although at the expense of more artifacts.

A simple bicone model is not a good match to the spatial and velocity distribution of the H i gas. Figures 1 and 2 clearly show the H i gas does not have a conical shape. The opening angle derived from line splitting in the cold gas by Leroy et al. (2015) is also smaller than values derived closer to the plane from the warm, ionized gas (Heckman et al. 1990; McKeith et al. 1995), which suggests the cold gas is more columnated. Another argument against a simple bicone is the presence of outflowing material at offset positions parallel to the minor axis. Figure 1 shows five slices parallel to the minor axis that are 1 kpc wide and offset at +2, +1, 0, −1, and −2 kpc along the midplane. Figure 8 shows the position–velocity at each of those five positions. This figure clearly shows the H i velocity profile observed within ±0.5 kpc of the minor axis (center) is also seen in positions offset by as much as ±2 kpc, albeit at somewhat lower intensity. A simple bicone with an opening angle of 8° (half-width) would have no flux within about ±7 kpc of the major axis in the offset positions at ±2 kpc. While the intensity at these offset positions could be explained with a larger opening angle, a larger opening angle would also produce a much larger range of velocities at each position that would not be consistent with the observations.

Zoom In Zoom Out Reset image size

We have consequently compared the data to a biconical frustum, rather than a simple bicone. Since a frustum is a conical shape with the top removed, we use this representation to model the origin of the outflow from a broader region of the disk than the single point of a cone. This shape is a better physical match to the spatial extent of the starburst, yet maintains the relatively narrow range of velocities observed further from the midplane (see Westmoquette et al. 2011 for the application of a frustum to NGC 253). We find that a frustum radius of 1 kpc at the midplane and an opening angle half-width of 8° is a reasonable match to the data. The outline of this biconical frustum is shown in Figure 1, and position–velocity diagrams for this model are shown on the right side of Figure 8 at the same five offset positions as the data on the left. The choice of a 1 kpc radius at the midplane was motivated by the physical extent of the starburst region (∼500 pc; Förster Schreiber et al. 2003), but increased to provide a better match to the observed position–velocity diagrams. The velocity of the gas in the model was calculated at each point based on the radial distance from the center of the model. The velocity includes gravity and drag forces as in Figure 6, and was normalized such that the initial velocity of the gas matches the observed H i kinematics at 1.5 kpc from the midplane. The model does not include any material closer than 1.5 kpc from the center, as we assumed this region is substantially affected by the disk kinematics. The main disagreement between the model and data are that the model has a somewhat broader velocity distribution. This could be because our simple model adopts a uniformly filled frustum. The model would better match the data if the material were more concentrated on the symmetry axis, although physically that would be in conflict with the expected location of the hot wind. At larger distances, the density is assumed to decrease as n ∝ r⁻².

It is also possible that the H i gas does not represent a continuous population of clouds flowing from the disk, but is instead partially or completely produced in situ. In this case, one potential origin is the dissociation of molecular gas. The molecular gas measurements from Leroy et al. (2015) within a few kiloparsecs of the midplane show similar kinematics as the atomic gas; however, the mass density of molecular gas is substantially lower than the H i off the midplane, and decreases much more rapidly with distance from M82. Deeper molecular gas observations would help to further test this hypothesis, along with calculations of the disassociation rate.

Another possibility is that the H i gas has radiatively cooled from the hot phase (Wang 1995a, 1995b). The basic physical picture, which was developed in detail by Thompson et al. (2016), starts with cool gas clouds in the disk that are initially accelerated and shredded by the hot wind. The addition of this cloud material then increases the mass loading in the wind, and then this additional mass seeds thermal and radiative instabilities that precipitate cool gas out in the wind at larger distances. Following the wind model of Chevalier & Clegg (1985), a key dimensionless parameter is the mass loading parameter β, which describes the ratio of the mass outflow rate ${\dot{M}}_{W}$ in the wind to the star formation rate ${\dot{M}}_{* }$, that is, ${\dot{M}}_{W}=\beta {\dot{M}}_{* }$. The radiative cooling requires some minimum mass loading due to the shape of the cooling function. For the case of M82, Hoopes et al. (2003) find weak evidence of cooling on large scales and Strickland & Heckman (2009) estimate a fairly low mass loading factor of β = 0.2–0.5 in the core, although the wind could be strongly mass loaded on larger scales. Thompson et al. (2016) therefore conclude that their model may not apply to M82, and note the additional complication of tidal material. The monotonic decline in the characteristic H i velocity with projected distance also does not seem consistent with radiative cooling because the velocity profile is inconsistent with the inferred velocity profile for the hot wind.

There is some phenomenological similarity between the velocity decrease in M82 on kiloparsec scales and the ionized gas kinematics observed in nearby AGNs on scales approximately an order of magnitude smaller. Studies of Seyfert galaxies such as NGC 1068 (Crenshaw & Kraemer 2000), Mrk 3 (Ruiz et al. 2001), and NGC 4151 (Das et al. 2005) are reasonably well matched with constant acceleration to 100–300 pc, followed by constant deceleration to the systemic velocity. Everett & Murray (2007) explored Parker wind models to explain the ∼100 pc outflows in AGNs, but found no self-consistent temperature profile that could explain all of the kinematics. They did note that the generic hydromagnetic winds explored by Matzner & McKee (1999) could explain the acceleration, and there are morphological similarities between that model and the observations of May & Steiner (2017), but it was not clear how to explain the deceleration without interaction with the surrounding medium.