ALMA [C i]³P₁–³P₀ Observations of NGC 6240: A Puzzling Molecular Outflow, and the Role of Outflows in the Global α_CO Factor of (U)LIRGs

With the aim of investigating the extent and morphology of the outflow, we produced interferometric maps of the CO(1–0), CO(2–1), and [C i](1–0) high-velocity emissions. The maps, shown in Figures 1(a)–(c), were generated by merging together and imaging the uv visibilities corresponding to the blue- and redshifted wings of the molecular lines, integrated respectively within v ∈ (−650, −200) km s⁻¹ and v ∈(250, 800) km s⁻¹. These are the velocity ranges that, following from the identification of the outflow components performed in Section 3.3 (and further discussed in Section 4.1 and Appendix B), are completely dominated by the emission from outflowing gas. The data displayed in panels (a)–(c) have a matched spatial resolution equal to $\sim 1\buildrel{\prime\prime}\over{.} 2$ (details in Section 2). Panels (d), (e) of Figure 1 show the maps of the blue and red [C i](1–0) line wings at the native spatial resolution of the ALMA Band 8 observations ($0\buildrel{\prime\prime}\over{.} 24,$ Section 2).

Zoom In Zoom Out Reset image size

The extended (>5 kpc) components of the outflow are best seen in Figures 1(a)–(c), whereas panels (d), (e) provide a zoomed view of the molecular wind in the inner 1–2 kpc. The bulk of the outflow extends eastward of the two AGNs, as already pointed out in the previous analysis of the CO(1–0) data done by Feruglio et al. (2013a). In addition, we identify for the first time a western extension of the molecular outflow, roughly aligned along the same east–west axis as the eastern component. At the sensitivity allowed by our data, we detect at a S/N > 5 CO emission features associated with the outflow up to a maximum distance of $13\buildrel{\prime\prime}\over{.} 3$ (6.8 kpc) and $7\buildrel{\prime\prime}\over{.} 2$ (3.7 kpc) from the nucleus, respectively, in the east and west directions (Figure 1(a)–(c)). The [C i](1–0) map in Figure 1(c) shows extended structures similar to CO, although the limited field of view of the Band 8 observations does not allow us to probe emission beyond a radius of $\sim 7\buildrel{\prime\prime}\over{.} 5$.

Feruglio et al. (2013a) hinted at the possibility that the redshifted CO detected in NGC 6240 could be involved in the feedback process, but did not explicitly ascribe it to the outflow because of the smaller spatial extent of the red wing with respect to the blue wing. In Figure 2 we directly compare the blue and red line wings using the ALMA CO(2–1) data. Based on their close spatial correspondence, whereby the red wing overlaps with the blue one across more than 7 kpc along the east–west direction, we conclude that both the redshifted and blueshifted velocity components trace the same massive molecular outflow. It follows that the eastern and western sides of the outflow are detected in both their approaching and receding components. At east, the blueshifted emission is brighter than the redshifted one and dominant beyond a 3.5 kpc radius.

Zoom In Zoom Out Reset image size

The ALMA [C i](1–0) data can be used to identify with high precision the location of the inner portion of the molecular outflow. Figures 1(d), (e) clearly show that the red wing peaks in the midpoint between the two AGNs, and that the blue wing has a maximum of intensity closer to the southern AGN, as already noted by Feruglio et al. (2013a). However, these data reveal for the first time that neither the blue nor the redshifted high-velocity [C i](1–0) emissions peak exactly at the AGN positions. The blue wing has a maximum of intensity at R.A.(J2000) = 16:52:58.8946 ± 0.0011s, decl.(J2000) = $+\mathrm{02.24.03.52}\pm 0\buildrel{\prime\prime}\over{.} 02$, offset by $0\buildrel{\prime\prime}\over{.} 18\,\pm 0\buildrel{\prime\prime}\over{.} 02$ to the north-east with respect to the southern AGN. The red wing instead peaks at R.A.(J2000) = 16:52:58.9224 ± 0.0007s, decl.(J2000) = $+\mathrm{02.24.04.0158}\pm 0\buildrel{\prime\prime}\over{.} 007$ (i.e., at an approximately equal distance of $0\buildrel{\prime\prime}\over{.} 8$ from the two AGNs). The [C i](1–0) red and blue wing peaks are separated by $0\buildrel{\prime\prime}\over{.} 65\pm 0\buildrel{\prime\prime}\over{.} 02$. This separation is consistent with the distance between the CO(2–1) peaks reported by Tacconi et al. (1999), although in that work they were interpreted as the signature of a rotating molecular gas disk.

The presence of such nuclear rotating H₂ structure has been largely debated in the literature, especially due to the very high CO velocity dispersion in this region (σ > 300 km s⁻¹), and to the mismatch between the dynamics of H₂ gas and stars (Gerssen et al. 2004; Engel et al. 2010). Following Tacconi et al. (1999) and Bryant and Scoville (1999), if a rotating disk is present, its signature should appear at lower projected velocities than those imaged in Figures 1(d), (e). In Figure 3 we show the high resolution intensity-weighted moment maps of the [C i](1–0) line emission within $-200\,\lt v[\mathrm{km}\,{{\rm{s}}}^{-1}]\lt 250$. The velocity field does not exhibit the characteristic butterfly pattern of a rotating disk, but it presents a highly asymmetric gradient whereby blueshifted velocities dominate the southern emission, whereas near-systemic and redshifted velocities characterize the northern emission. These features in the velocity field are correlated in both velocity and position with the high-velocity wings (Figures 1(d), (e)). Furthermore, the right panel of Figure 3 shows that the velocity dispersion is uniform ($50\lesssim {\sigma }_{v}[\mathrm{km}\,{{\rm{s}}}^{-1}]\lesssim 80$) throughout the entire source and enhanced (σ_v ≥ 100 km s⁻¹) in a central hourglass-shaped structure extending east–west, which is the same direction of expansion of the larger-scale outflow. This structure has a high-σ_v peak with σ_v ≥ 150 km s⁻¹ to the east and another possible peak to the west with σ_v ≥ 130 km s⁻¹. Such high-σ_v points coincide with the blue and redshifted velocity peaks detected in the moment 1 map. Therefore, based on Figure 3, we conclude that the molecular gas emission between the two AGNs is dominated by a nuclear outflow expanding east–west and connected to the larger-scale outflow shown in Figure 1. The regions of enhanced turbulence may represent the places where the outflow opening angle widens up, hence increasing the line of sight velocity dispersion and velocity of the molecular gas.

Zoom In Zoom Out Reset image size

Figure 4 shows the CO(1–0), CO(2–1), and [C i](1–0) spectra extracted from the $12^{\prime\prime} \times 6^{\prime\prime}$-size rectangular aperture reported in Figure 1(a), encompassing both the nucleus and the extended molecular outflow of NGC 6240. The analysis of these integrated spectra, described later, is aimed at deriving a source-averaged ${\alpha }_{\mathrm{CO}}$ for the quiescent and outflowing molecular ISM in NGC 6240. Such analysis is included here because it can be useful as a reference for unresolved observations (e.g., high redshift analogues of this merger). We stress, however, that the quality of our data, the proximity of the source, and its large spatial extent allow us to perform a much more detailed, spatially resolved analysis. The latter will be presented in Section 3.3 and delivers the most reliable ${\alpha }_{\mathrm{CO}}$ values for the outflow and the quiescent gas.

Zoom In Zoom Out Reset image size

The spectra in Figure 4 were fitted simultaneously using two Gaussians to account for the narrow core and broad wings of the emission lines, by constraining the central velocity (v) and velocity dispersion (σ_v) of each Gaussian to be equal in the three transitions. Table 1 reports the best fit results and the corresponding line luminosities calculated from the integrated fluxes following Solomon et al. (1997). The [C i](1–0) line luminosities listed in Table 1 are employed to measure the molecular gas mass (M_mol, including the contribution from helium) associated with the narrow and broad line components. The expression for local thermodynamic equilibrium (LTE; i.e., uniform T_ex) and optically thin emission (${\tau }_{[{\rm{C}}{\rm{I}}](1\mbox{--}0)}\ll 1$), assuming a negligible background (CMB temperature, T_CMB ≪ T_ex) and the Rayleigh–Jeans approximation (hν_{[C i](1–0)} ≪ kT_ex), is

$$\begin{eqnarray}{M}_{\mathrm{mol}}[{M}_{\odot }] & = & (4.31\cdot {10}^{-5})\cdot {X}_{{\rm{C}}\,{\rm{I}}}^{-1}\cdot (1+\\ & & 3{e}^{-23.6/{T}_{\mathrm{ex}}[{\rm{K}}]}+5{e}^{-62.5/{T}_{\mathrm{ex}}[{\rm{K}}]})\cdot \\ & & {e}^{23.6/{T}_{\mathrm{ex}}[{\rm{K}}]}\cdot {L}_{[{\rm{C}}\,{\rm{}}{\rm{I}}](1\mbox{--}0)}^{{\prime} }[{\rm{K}}\,\mathrm{km}\,{{\rm{s}}}^{-1}\,{\mathrm{pc}}^{2}],\end{eqnarray} \tag{ 1 }$$

where X_{C i} is the $[{\rm{C}}\,{\rm{I}}/{{\rm{H}}}_{2}]$ abundance ratio and T_ex is the excitation temperature of the gas (see detailed explanations by Papadopoulos et al. 2004 and Mangum & Shirley 2015). We adopt T_ex = 30 K and X_{C i} = (3.0 ± 1.5) × 10⁻⁵, which are appropriate for (U)LIRGs (Weiß et al. 2003, 2005; Papadopoulos et al. 2004; Walter et al. 2011; Jiao et al.2017). These assumptions will be further discussed in Section 4.1.2.

Table 1.  Results of the Simultaneous Fit to the Total Spectra^a

	CO(1–0)	CO(2–1)	[C i](1–0)
	Narrow component
vb [km s−1 ]	−9.1 ± 1.0	−9.1 ± 1.0	−9.1 ± 1.0
${\sigma }_{v}$ [km s−1 ]	101.0 ± 1.0	101.0 ± 1.0	101.0 ± 1.0
Speak [mJy]	211 ± 4	1058 ± 12	1360 ± 70
$\int {S}_{v}{dv}$ [Jy km s−1 ]	53.4 ± 1.0	268 ± 4	344 ± 18
${L}^{{\prime} }$ [109 K km s−1 pc2]	1.55 ± 0.03	1.94 ± 0.03	0.55 ± 0.03
	Broad component
vb [km s−1 ]	78.0 ± 1.2	78.0 ± 1.2	78.0 ± 1.2
σv [km s−1 ]	264.4 ± 1.0	264.4 ± 1.0	264.4 ± 1.0
Speak [mJy]	301 ± 3	1458 ± 12	1040 ± 40
$\int {S}_{v}{dv}$ [Jy km s−1 ]	199 ± 2	966 ± 9	690 ± 30
L' [109 K km s−1 pc2]	5.78 ± 0.06	7.01 ± 0.06	1.09 ± 0.05
	Total line
$\int {S}_{v}{dv}$ [Jy km s−1]	253 ± 2	1234 ± 10	1030 ± 30
L' [109 K km s−1 pc2]	7.33 ± 0.07	8.96 ± 0.07	1.64 ± 0.05

Notes.

^aThe errors quoted in this table are purely statistical and do not include the absolute flux calibration uncertainty. ^bWe employ the optical Doppler definition. The fit allows for a global velocity shift of CO(2–1) and [C i](1–0) with respect to CO(1–0) to take into account the different spectral binning. The best fit returns: ${v}_{\mathrm{CO}(2-1)}-{v}_{\mathrm{CO}(1-0)}=10.2\,\pm 0.9$ km s⁻¹ and ${v}_{[{\rm{C}}{\rm{I}}](1-0)}-{v}_{\mathrm{CO}(1-0)}=-18\pm 4$ km s⁻¹.

Download table as:  ASCIITypeset image

By defining a [C i](1–0)-to-H₂ conversion factor (α_{[C i]}) in analogy with the commonly employed ${\alpha }_{\mathrm{CO}}$ factor (Bolatto et al. 2013), Equation (1) resolves into

$$\begin{eqnarray}{M}_{\mathrm{mol}}[{M}_{\odot }] & \equiv & {\alpha }_{[{\rm{C}}{\rm{I}}]}\,{L}_{[{\rm{C}}\,{\rm{I}}]}^{{\prime} }[{\rm{K}}\,\mathrm{km}\,{{\rm{s}}}^{-1}\,{\mathrm{pc}}^{2}]\\ & & \mathrm{with}\,{\alpha }_{[{\rm{C}}{\rm{I}}]}=9.43\,[{M}_{\odot }\,{({\rm{K}}\mathrm{km}{{\rm{s}}}^{-1}{\mathrm{pc}}^{2})}^{-1}].\end{eqnarray} \tag{ 2 }$$

Using the values in Table 1 and Equation (2), we obtain a total molecular gas mass of ${M}_{\mathrm{mol}}=(1.5\pm 0.8)\cdot {10}^{10}\,{M}_{\odot }$, of which $(5\pm 3)\cdot {10}^{9}\,{M}_{\odot }$ is in the narrow component, and $(10\pm 5)\cdot {10}^{9}\,{M}_{\odot }$ in the broad wings. We then use these M_mol values to estimate ${\alpha }_{\mathrm{CO}}$:

$$\begin{eqnarray}&&{\alpha }_{\mathrm{CO}}={M}_{\mathrm{mol}}[{M}_{\odot }]\,{({L}_{\mathrm{CO}(1\mbox{--}0)}^{{\prime} }[{\rm{K}}\mathrm{km}{{\rm{s}}}^{-1}{\mathrm{pc}}^{2}])}^{-1}.\end{eqnarray} \tag{ 3 }$$

The results are reported in the first three rows of Table 2 for the the total, narrow, and broad emissions in NGC 6240. The so-derived ${\alpha }_{\mathrm{CO}}$ factors differ between the narrow and broad line components, being a factor of 1.8 ± 0.5 lower in the latter.²⁰ In the narrow component, the ${\alpha }_{\mathrm{CO}}$ is significantly higher than the typical (U)LIRG value. As mentioned in Section 1, higher ${\alpha }_{\mathrm{CO}}$ values become possible in (U)LIRGs if a significant fraction of the mass is "hidden" in dense and bound H₂ clouds. This is probably the case of NGC 6240, in which a large study using CO SLEDs from J = 1–0 up to J = 13–12 from the Herschel Space Observatory, as well as multi-J HCN, CS, and HCO⁺ line data from ground-based observatories, finds ${\alpha }_{\mathrm{CO}}$ $\sim \,2\mbox{--}4\,{M}_{\odot }\,{({\rm{K}}\mathrm{km}{{\rm{s}}}^{-1}{\mathrm{pc}}^{2})}^{-1}$ (Papadopoulos et al. 2014), consistent with our estimates.

Table 2.  ${\alpha }_{\mathrm{CO}}$ and r₂₁ Values^a

	${\alpha }_{\mathrm{CO}}$	r21
	$[{M}_{\odot }{({\rm{K}}\mathrm{km}{{\rm{s}}}^{-1}{\mathrm{pc}}^{2})}^{-1}]$
Totalb	2.1 ± 1.1	1.22 ± 0.14
Totalb narrow comp	3.3 ± 1.8	1.25 ± 0.18
Totalb broad comp	1.8 ± 0.9	1.21 ± 0.17
Meanc global	2.5 ± 1.4	1.17 ± 0.19
Meanc systemic comp	3.2 ± 1.8	1.0 ± 0.2
Meanc outflow comp	2.1 ± 1.2	1.4 ± 0.3

Notes.

^aQuoted errors are dominated by systematic uncertainties (e.g., absolute flux calibration errors, error on X_{C i}). ^bCalculated from the simultaneous fit to the total CO(1–0), CO(2–1), and [C i](1–0) spectra shown in Figure 4, whose results are reported in Table 1 (details in Section 3.2). ^cMean values calculated from the simultaneous fit to the CO(1–0), CO(2–1), and [C i](1–0) spectra extracted from the grid of 13 boxes shown in Figure 1(a), as explained in Section 3.3. The corresponding spectral fits are shown in Appendix B (Figures 8–10).

Download table as:  ASCIITypeset image

Table 2 lists also the CO(2–1)/CO(1–0) luminosity ratios, defined as

$$\begin{eqnarray}&&{r}_{21}\equiv {L}_{\mathrm{CO}(2-1)}^{{\prime} }/{L}_{\mathrm{CO}(1-0)}^{{\prime} }.\end{eqnarray} \tag{ 4 }$$

We find r₂₁ as consistently ∼1.2—hence higher than unity (at the 1.5σ level)—for both the narrow and broad Gaussian components. Since our total CO(1–0) and CO(2–1) fluxes are consistent with previous measurements (Costagliola et al. 2011; Papadopoulos et al. 2012b; Saito et al. 2018), we exclude that spatial filtering due to an incomplete uv coverage is significantly affecting the r₂₁ values. As pointed out by Papadopoulos et al. (2012b), galaxy-averaged r₂₁ > 1 values are not uncommon in (U)LIRGs and are indicative of extreme gas conditions. In this analysis of the integrated spectra, we derived r₂₁ > 1 in both the broad and narrow components. However, as we will show in Sections 3.3 and 3.5, the spatially resolved analysis will reveal that r₂₁ ≳ 1 values are typical of the outflowing gas and in general of higher-σ_v components, while the "quiescent" ISM has r₂₁ ∼ 1.

The previous analysis (Section 3.2) was based on the spectral fit shown in Figure 4, where we decomposed the total molecular line emission into a narrow and a broad Gaussian. In first approximation, these two spectral components can be respectively identified with the quiescent and outflowing molecular gas reservoirs of NGC 6240. However, the superb S/N and spatial resolution of our data allow us to take this analysis one step further and refine the definition of quiescent and outflowing components. This is done by including the spatially resolved information provided by the interferometric data, as described later.

We divide the central $12^{\prime\prime} \times 6^{\prime\prime}$ region employed in the previous analysis into a grid of 13 squared boxes and use them as apertures to extract the corresponding CO(1–0), CO(2–1), and [C i](1–0) spectra. As shown in Figure 1(a), the central nine boxes have a size of $2^{\prime\prime} \times 2^{\prime\prime}$, while the external four boxes have a size of $3^{\prime\prime} \times 3^{\prime\prime}$. The box spectra are presented in Appendix B (Figures 8–10).

For each box, the CO(1–0), CO(2–1), and [C i](1–0) spectra are fitted simultaneously with a combination of Gaussian functions tied to have the same line centers and widths for all three transitions. In the fitting procedure, we minimize the number of spectral components required to reproduce the line profiles, up to a maximum of four Gaussians per box. The Gaussian functions employed by the simultaneous fit span a wide range in FWHM and velocity, shown in Figure 5. The next step is to classify each of these components as "systemic" or "outflow." In many local (U)LIRGs, molecular outflows can be traced through components whose kinematical and spatial features deviate from a rotating molecular structure (Cicone et al. 2014; García-Burillo et al. 2014). However, in this source we do not detect any clear velocity gradient that may indicate the presence of a rotating molecular gas disk (Figure 3). Therefore, we adopt a different method and identify as "quiescent" the gas probed by the spectral narrow line components that are detected throughout the entire source extent (Figures 8–10). Our simultaneous fit to the CO(1–0), CO(2–1), and [C i](1–0) spectra returns for these narrow components typical FWHM and central velocities in the ranges: FWHM < 400 km s⁻¹ and $-200\lt v[\mathrm{km}\,{{\rm{s}}}^{-1}]\lt 250$, consistent with what found by Feruglio et al. (2013b).²¹ Based on these results, we assume that all components with $-200\lt v[\mathrm{km}\,{{\rm{s}}}^{-1}]\lt 250$ and FWHM < 400 km s⁻¹ trace quiescent gas that is not involved in the outflow. These constraints correspond to the region of the FWHM-v parameter space delimited by the blue dashed lines in Figure 5. All components outside this rectangular area are classified as "outflow." These assumptions are discussed in detail and validated in Appendix B, whereas more general considerations about our outflow identification method are reported in Section 4.1.1.

Zoom In Zoom Out Reset image size

Using the results of the simultaneous fit, we measure, for each box and for each of the CO(1–0), CO(2–1), and [C i](1–0) transitions, the velocity-integrated fluxes apportioned in the "systemic" and "outflow" components. These are computed by summing the fluxes from the respectively classified Gaussian functions fitted to the molecular line profiles. For example, in the case of the central box (labeled as "C1" in Figures 8–10), the simultaneous fit employs three Gaussians: a narrow one classified as "systemic," and two additional ones classified as "outflow." The flux of the first Gaussian corresponds to the flux of the "systemic" component for this box, whereas the total "outflow" component flux is given by the sum of the fluxes of the other two Gaussians (the errors are added in quadrature).

The velocity-integrated fluxes (total, systemic, outflow) are then converted into line luminosities and, from these, the corresponding ${\alpha }_{\mathrm{CO}}$ and r₂₁ can be derived by following the same steps as in Section 3.2 (Equations (2)–(4)). Table 2 (bottom three rows) lists the resulting mean values of ${\alpha }_{\mathrm{CO}}$ and r₂₁ obtained from the analysis of all 13 boxes. In computing the mean, we only include the components detected at a S/N ≥ 3 in each of the transitions used to calculate ${\alpha }_{\mathrm{CO}}$ or r₂₁—that is, CO(1–0) and [C i](1–0) for the former, and CO(1–0) and CO(2–1) for the latter. The new ${\alpha }_{\mathrm{CO}}$ values derived for the systemic and outflowing components, respectively equal to 3.2 ± 1.8 and 2.1 ± 1.2, are perfectly consistent with the previous analysis based on the integrated spectra. Instead, this new analysis delivers different $\langle {r}_{21}\rangle$ values for the systemic (1.0 ± 0.2) and outflowing component (1.4 ± 0.3), although still consistent if considering the associated uncertainties (dominated by the flux calibration errors).

By using the CO(1–0) line data and summing the contribution from all boxes, including both the systemic and the outflowing components, we derive a total molecular gas mass of ${M}_{\mathrm{mol}}^{\mathrm{tot}}=(2.1\pm 0.5)\times {10}^{10}\,{M}_{\odot }$. To compute the M_mol within each box we adopt, when available, the "global" ${\alpha }_{\mathrm{CO}}$ factor estimated for that same box; otherwise, we use the mean value of ${\alpha }_{\mathrm{CO}}$ = 2.5 ± 1.4.²² Compared with previous works recovering the same amount of CO flux, our new ${M}_{\mathrm{mol}}^{\mathrm{tot}}$ estimate is higher than in Tacconi et al. (1999) and Feruglio et al. (2013a), but consistent with Papadopoulos et al. (2014).

In this section we use the results of the spatially resolved spectral analysis presented in Section 3.3 to constrain the mass (M_out), mass-loss rate (${\dot{M}}_{\mathrm{out}}$), kinetic power ($1/2{\dot{M}}_{\mathrm{out}}{v}^{2}$), and momentum rate (${\dot{M}}_{\mathrm{out}}v$) of the molecular outflow. We first select the boxes in which an outflow component is detected in the CO(1–0) spectrum with S/N ≥ 3. As described in Section 3.3, the outflow component is defined as the sum of all Gaussian functions employed by the simultaneous fit that lie outside the rectangular region of the FWHM-v parameter space shown in Figure 5. With this S/N ≥ 3 constraint, 12 boxes (that is, all except W1) are selected to have an outflow component in CO(1–0), and for each box,²³ we measure the following:

• (i)  
  The average outflow velocity (v_out,i), equal to the mean of the (moduli of the) central velocities of the individual Gaussians classified as "outflow"
• (ii)  
  The molecular gas mass in outflow (M_out,i), calculated by multiplying the ${L}_{\mathrm{CO}(1\mbox{--}0)}^{{\prime} }$ of each outflow component by an appropriate ${\alpha }_{\mathrm{CO}}$. In ten boxes, the outflow component is detected with S/N ≥ 3 also in the [C i](1–0) transition; hence for these boxes we can use their corresponding ${\alpha }_{\mathrm{CO}}$ factor (see Section 3.3). For the remaining two boxes (E1 and W2), we adopt the galaxy-averaged outflow ${\alpha }_{\mathrm{CO}}$ of 2.1 ± 1.2 (Table 2).
• (iii)  
  The dynamical timescale of the outflow, defined as τ_dyn,i = R_i/v_out,i, where R_i is the distance of the center of the box from R.A. = 16:52:58.900, decl. = 02.24.03.950. This definition cannot be applied to the central box (C1) because the so-estimated R would be zero, hence boosting the mass-loss rate to infinite. Therefore, for box C1, we conservatively assume that most of the outflow emission comes from a radius of $1^{\prime\prime}$; hence we set R = 0.5 kpc. For all boxes, we assume the uncertainty on R_i to be $0\buildrel{\prime\prime}\over{.} 6$ (0.3 kpc), which is half a beam size.
• (iv)  
  The mass-loss rate ${\dot{M}}_{\mathrm{out},i}$, equal to M_out,i/τ_dyn,i

All uncertainties are derived by error propagation.

The resulting total outflow mass and mass-loss rate, obtained by adding the contribution from all boxes, are respectively ${M}_{\mathrm{out}}=(1.2\pm 0.3)\times {10}^{10}\,{M}_{\odot }$ and dM_out/dt = 2500 ± 1200 ${M}_{\odot }\,{\mathrm{yr}}^{-1}$. As discussed in Appendix B, the largest contribution to both ${\dot{M}}_{\mathrm{out}}$ and its uncertainty is given by the central box. Indeed, box C1 has at the same time the highest estimated M_out and the smallest—and most uncertain—R, because the outflow is launched from within this region, likely close to the midpoint between the two AGNs, as suggested by Figures 1(d), (e).

Similar to the mass-loss rate, we calculate the total kinetic power and momentum rate of the outflow by summing the contribution from all boxes with a CO(1–0) outflow component, and we obtain, respectively, $1/2{\dot{M}}_{\mathrm{out}}{v}_{\mathrm{out}}^{2}\,\equiv {\sum }_{i}1/2{\dot{M}}_{\mathrm{out},i}{v}_{\mathrm{out},i}^{2}=(0.033\pm 0.019)\,{L}_{\mathrm{AGN}}$ and $v{\dot{M}}_{\mathrm{out}}\,\equiv {\sum }_{i}{v}_{i}{\dot{M}}_{\mathrm{out},i}=(80\pm 50)\,{L}_{\mathrm{AGN}}/c$. If all the gas carried by the outflow escaped the system and the mass-loss continued at the current rate, the depletion timescale of the molecular gas reservoir in NGC 6240 would be τ_dep = 8 ± 4 Myr. All the relevant numbers describing the properties of the source and of the molecular outflow are reported in Table 3 and will be discussed in Section 4.3 in the context of feedback models.

Table 3.  Summary of Source and Outflow Properties

Source properties
${L}_{\mathrm{TIR}(8\mbox{--}1000\mu {\rm{m}})}$ [erg s−1]	2.71 × 1045a
LBol [erg s−1]	3.11 × 1045b
LAGN [erg s−1]	(1.1 ± 0.4) × 1045c
αAGN ≡ LAGN/LBol	0.35 ± 0.13
SFR [M${}_{\odot }$ yr−1]	46 ± 9d
M${}_{\mathrm{mol}}^{\mathrm{tot}}$ [M${}_{\odot }]$	(2.1 ± 0.5) × 1010e
Molecular outflow properties (estimated in Section 3.4)
rmax [kpc]	2.4 ± 0.3f
$\langle {v}_{\mathrm{out}}\rangle$ [km s−1 ]	250 ± 50g
$\langle {\sigma }_{\mathrm{out}}\rangle$ [km s−1]	220 ± 20g
$\langle {\tau }_{\mathrm{dyn}}\rangle$ [Myr]	6.5 ± 1.8g
Mout [M${}_{\odot }]$	(1.2 ± 0.3) × 1010
${\dot{M}}_{\mathrm{out}}$ [M${}_{\odot }$ yr−1]	2500 ± 1200
$v{\dot{M}}_{\mathrm{out}}$ $[{\rm{g}}\,\mathrm{cm}\,{{\rm{s}}}^{-2}]$	(3.1 ± 1.2) × 1036
$1/2{\dot{M}}_{\mathrm{out}}{v}^{2}$ $[\mathrm{erg}\,{{\rm{s}}}^{-1}]$	(3.6 ± 1.6) × 1043
$\eta \equiv {\dot{M}}_{\mathrm{out}}$/SFR	50 ± 30
($v{\dot{M}}_{\mathrm{out}}$)/(LAGN/c)	80 ± 50
($1/2{\dot{M}}_{\mathrm{out}}{v}^{2}$)/LAGN	0.033 ± 0.019
${\tau }_{\mathrm{dep}}\equiv {M}_{\mathrm{mol}}^{\mathrm{tot}}/{\dot{M}}_{\mathrm{out}}$ [Myr]	8 ± 4

Notes.

^aFrom the IRAS Revised Bright Galaxy Sample (Sanders et al. 2003). ^bL_Bol = 1.15 L_TIR, following Veilleux et al. (2009). ^cTotal bolometric luminosity of the dual AGN system estimated from X-ray data by Puccetti et al. (2016). ^d $\mathrm{SFR}=(1-{\alpha }_{\mathrm{AGN}})\times {10}^{-10}$ ${L}_{\mathrm{TIR}}$, following Sturm et al. (2011). ^eTotal molecular gas mass in the $12^{\prime\prime} \times 6^{\prime\prime}$ region encompassing the nucleus and the outflow, derived in Section 3.4. ^fMaximum distance at which we detect [C i](1–0) in the outflow at a S/N > 3; hence the quoted r_max should be considered a lower limit constraint allowed by current data. ^gMean values obtained from the analysis of all boxes.

Download table as:  ASCIITypeset image

Very stringent lower limits on the outflow energetics can be derived by assuming that its CO(1–0) emission is fully optically thin. For optically thin gas and T_ex = 30 K, the ${\alpha }_{\mathrm{CO}}$ factor would be ∼0.34 (Bolatto et al. 2013), and the outflow mass and mass-loss rate would be ${M}_{\mathrm{out}}=(1.98\pm 0.09)\times {10}^{9}\,{M}_{\odot }$ and dM_out/dt = 430 ± 160 ${M}_{\odot }\,{\mathrm{yr}}^{-1}$. However, we stress that the assumption of fully optically thin CO(1–0) emission in the outflow is not supported by our data, which instead favor an ${\alpha }_{\mathrm{CO}}$ factor for outflowing gas that is intermediate between the optically thin and the optically thick values (for solar metallicities).

Using the results of the spatially resolved analysis presented in Section 3.3, we now study how the ${\alpha }_{\mathrm{CO}}$ and r₂₁ parameters vary as a function of velocity dispersion (σ_v) and projected distance (d) from the nucleus of NGC 6240. The relevant plots are shown in Figure 6. To investigate possible statistical correlations, we conduct a Bayesian linear regression analysis of the relations in Figure 6, following Kelly (2007).²⁴

Zoom In Zoom Out Reset image size

The left panels of Figure 6 do not indicate any statistically significant relation between ${\alpha }_{\mathrm{CO}}$ and either σ_v or d. Instead, they show that the ${\alpha }_{\mathrm{CO}}$ factor is systematically higher—although formally only at a significance of 1.2σ (Table 2)—in the quiescent gas than in the outflow, regardless of the velocity dispersion of the clouds, or of their position with respect to the merger nucleus. For the non-outflowing components, the ${\alpha }_{\mathrm{CO}}$ factors are at least twice the so-called (U)LIRG value (Downes & Solomon 1998), and reach up to Galactic values. This result is consistent with the multi-transition analysis by Papadopoulos et al. (2014), and is likely due to the state of the dense gas phase that low-J CO lines alone cannot constrain, but which instead is accounted for when using [C i](1–0) as a molecular mass tracer. Nevertheless, the outflowing H₂ gas has lower ${\alpha }_{\mathrm{CO}}$ values than the quiescent ISM. This is indeed expected from the ISM physics behind ${\alpha }_{\mathrm{CO}}$ for warm and strongly unbound gas states (Papadopoulos et al. 2012a)—that is, the type of gas that we expect to be embedded in outflows. In particular, in the case that molecular outflows are ubiquitous in (U)LIRGs, as suggested by observations (Sturm et al. 2011; Spoon et al. 2013; Veilleux et al. 2013; Cicone et al. 2014), the outflow may be the location of the diffuse and warm molecular gas phase that is not contained in self-gravitating cooler clouds—a sort of "intercloud" medium advocated by some of the previous analyses based solely on low-J CO, ¹³CO line observations (Aalto et al. 1995; Downes & Solomon 1998). Furthermore, the flat trend between ${\alpha }_{\mathrm{CO}}$ and d observed in Figure 6 does not support the hypothesis that the lower ${\alpha }_{\mathrm{CO}}$ values in (U)LIRGs are related to the collision of the progenitors' disks, since in this case we would naively expect the lower ${\alpha }_{\mathrm{CO}}$ clouds to be concentrated in the central regions of the merger. The ${\alpha }_{\mathrm{CO}}$ values measured for the outflow components are, however, significantly higher than the optically thin value, suggesting that not all of the outflowing material is diffuse and warm, but there may still be a significant amount of dense gas. These results are further discussed and contextualized in Section 4.2.

The right panels of Figure 6 show a weak correlation between the r₂₁ and σ_v (correlation coefficient, ρ = 0.4 ± 0.2) and an anti-correlation with the distance, although only for the systemic/quiescent components (ρ = −0.7 ± 0.2). The corresponding best fit relations, of the form r₂₁ = α + βx, plotted in Figure 6, have (α, β) = (0.8 ± 0.2, 2.1 ± 1.4 × 10⁻³) for x = σ_v, and (α, β) = (1.5 ± 0.2, −0.33 ± 0.13) for x = d. The systemic ISM shows 0.8 ≲ r₂₁ ≲ 1.4, whereas the outflow is characterized by higher ratios, with most components in the range 1.2 ≲ r₂₁ ≲ 2.5, although we observe a large spread in r₂₁ values at d > 2 kpc.

CO(2–1)/CO(1–0) luminosity ratios of r₂₁ ∼ 0.8–1.0 are typically found in the molecular disks of normal spiral galaxies (Leroy et al. 2009) and are indicative of optically thick CO emission with T_kin ∼ 10–30 K (under LTE assumptions). Nevertheless, such low-J CO line ratios, in absence of additional transitions, are well-known to be highly degenerate tracers of the average gas physical conditions. Higher-J data of CO, molecules with larger dipole moment, and isotopologues can break such degeneracies. Such studies exist for NGC 6240 (Greve et al. 2009; Meijerink et al 2013; Papadopoulos et al. 2014), and found extraordinary states for the molecular gas, with average densities typically above 10⁴ cm⁻³ and temperatures T_kin ∼ 30–100 K.

On the contrary, global CO(2–1)/CO(1–0) ratios exceeding unity have a lower degree of degeneracy in terms of the extraordinary conditions that they imply for molecular gas, as they require warmer (T_kin ≳ 100 K) and/or strongly unbound states (Papadopoulos et al. 2012b). In NGC 6240, optical depth effects are most likely at the origin of the ${r}_{21}\gt 1$ values. More specifically, such ratios can result from highly non-virial motions (e.g., the large velocity gradients of the outflowing clouds), causing the CO lines to become partially transparent, as also supported by the tentative trend of increasing r₂₁ with σ_v (Figure 6). This finding independently strengthens our explanation for the lower ${\alpha }_{\mathrm{CO}}$ factors derived for the outflowing gas, which are intermediate between an optically thin and an optically thick value (for typical solar CO abundances).

Our results build, on the one hand, on the identification of the outflow components, and on the other hand on the assumption that CO(1–0) and [C i](1–0) trace the same molecular gas, implying that M_mol can be measured from [C i](1–0). In this section we further comment on these steps and discuss their caveats and limitations.

4.1.1. The Outflow Identification

4.1.2. Combining [C i](1–0) and CO(1–0) Data to Infer ${\alpha }_{\mathrm{CO}}$ and M_mol

The average ${\alpha }_{\mathrm{CO}}$ factors that we measure for the quiescent and outflowing components of the ISM in NGC 6240 (Table 2) are both higher than the classic (U)LIRG ${\alpha }_{\mathrm{CO}}$. How can we reconcile this result with previous works advocating for significantly lower ${\alpha }_{\mathrm{CO}}$ values in (U)LIRGs? In the case of NGC 6240, our analysis has highlighted several effects that may have plagued previous ${\alpha }_{\mathrm{CO}}$ estimates:

• 1.  
  The widespread presence of outflowing gas implies that, at any location within this merger, the molecular line emission includes a significant contribution from the outflow, with its overall lower ${\alpha }_{\mathrm{CO}}$ and higher r₂₁. As a result, an analysis of the global ISM conditions (especially if based only on low-J CO lines; see also point 3 of this list) would get contaminated by the warm unbound H₂ envelopes in the outflow, and their larger ${L}_{\mathrm{CO}}^{{\prime} }/{M}_{{\rm{H}}2}$ ratios would drive down the global ${\alpha }_{\mathrm{CO}}$ estimate (Yao et al. 2003; Papadopoulos et al. 2012a).
• 2.  
  The outflow dominates the velocity field of the H₂ gas throughout the entire source, including the central region (Figure 3). The apparent nuclear north–south velocity gradient identified in previous CO line data (Bryant & Scoville 1999; Tacconi et al. 1999) is actually not compatible with ordered rotation once observed at higher spatial resolution, but it shows several features distinctive of the outflow. Therefore, the assumption that the molecular gas in this area is dominated by ordered motions is broken, making any dynamical mass estimate unreliable (if the outflow is not properly taken into account).
• 3.  
  Previous analyses based only on low-J CO lines have probably missed a substantial fraction of the denser gas phase that is instead accounted for when using the optically thin [C i](1–0) line as a gas mass tracer, or when probing the excitation of the ISM using high-J CO transitions and high density molecular gas tracers. Indeed, the ${\alpha }_{\mathrm{CO}}$ value derived for the quiescent gas reservoir is consistent with a significant contribution from a dense gas state. Even in the outflow, the average ${\alpha }_{\mathrm{CO}}$ is still significantly higher than the optically thin value; hence the presence of dense gas may not be negligible. A conspicuous dense gas phase has already been demonstrated in a few other galaxy-wide molecular outflows (Aalto et al. 2012; García-Burillo et al. 2014; Sakamoto et al. 2014; Alatalo et al. 2015), and its presence would make more likely the formation of stars within these outflows (Maiolino et al. 2017).
• 4.  
  The molecular gas emission in NGC 6240 is clearly very extended—both spectrally and spatially. As a result, at least some of the previous interferometric observations (especially "pre ALMA") may have been severely affected by incomplete uv coverages filtering out the emission on larger scales, hence impacting on the measured line fluxes and sizes. Furthermore, as already noted by Tacconi et al. (1999), an insufficient spectral bandwidth may have hindered a correct baseline and/or continuum fitting and subtraction. The latter can be an issue for both single-dish and interferometric observations, including observations with ALMA if only one spectral window is employed to sample the line.

Since molecular outflows are a common phenomenon in local (U)LIRGs (Sturm et al. 2011; Spoon et al. 2013; Veilleux et al. 2013; Cicone et al. 2014; Fluetsch et al. 2018), at least some of the previously listed considerations may be generalized to their entire class. Therefore, it is possible that the so-called (U)LIRG ${\alpha }_{\mathrm{CO}}$ factor is an artifact resulting from modeling the molecular ISM of such sources containing massive H₂ outflows.

The extreme spatial extent of its molecular outflow makes NGC 6240 one of the few sources—all powerful quasars—hosting H₂ outflows with sizes of ≳10 kpc (Veilleux et al. 2017; see also the 30 kpc size [C ii]λ158 μm outflow at z = 6.4 studied by Cicone et al. 2015). In comparison, the H₂ gas entrained in the well-studied starburst-driven winds of M82 and NGC 253 reaches maximum scales of ∼1–2 kpc (Walter et al. 2002, 2017). Furthermore, among all the large-scale molecular outflows discovered so far in quasar host galaxies, the outflow of NGC 6240 is the one that has been observed at the highest spatial resolution (∼120 pc). Indeed, the ALMA [C i](1–0) line data allowed us to probe deep into the nuclear region of the merger, close to the launching point of the molecular wind, and surprisingly revealed that the outflow emission peaks between the two AGNs rather than on either of the two. This is apparently at odds with an AGN radiative-mode feedback scenario, in which the multiphase outflow is expected to be generated close to the central engine (Costa et al. 2015).

Nevertheless, the role of the AGN(s) is certified by the extreme energetics of the molecular outflow, which have been constrained here with unprecedented accuracy. By comparing our M_out and M_mol^tot estimates (Table 3), it appears that 60 ± 20% of the molecular medium is involved in the outflow. The estimated mass-loss rate of 2500 ± 1200 ${M}_{\odot }\,{\mathrm{yr}}^{-1}$ corresponds to $\eta \equiv {\dot{M}}_{\mathrm{out}}/\mathrm{SFR}=50\pm 30$, whereas the lower limit on ${\dot{M}}_{\mathrm{out}}$, calculated using the optically thin ${\alpha }_{\mathrm{CO}}$ prescription, corresponds to $\eta \equiv {\dot{M}}_{\mathrm{out}}/\mathrm{SFR}=9\pm 4$. Such high mass loading factors are inconsistent with a purely star formation-driven wind. As a matter of fact, stellar feedback alone can hardly bear outflows with η much higher than unity. Cosmological hydrodynamical simulations incorporating realistic stellar feedback physics, by including mechanisms other than supernovae, can reach up to η ∼ 10 (Hopkins et al. 2012). However, because η in these simulations anti-correlates with the mass of the galaxy, the highest η values are generally predicted for dwarf galaxies, whereas for galaxies with baryonic masses of several 10¹⁰ M_⊙ such as NGC 6240, the η achievable by stellar feedback can be at maximum 2–3.

Based on its energetics, we can therefore rule out that the massive molecular outflow observed in NGC 6240 is the result of star formation alone. The outflow energetics can instead be fully accommodated within the predictions of AGN feedback models (Faucher-Giguère & Quataert 2012; Costa et al. 2014; Zubovas & King 2014). In addition, these models can explain the multi-wavelength properties of NGC 6240. At optical wavelengths, NGC 6240 is known to host a ionized wind (Heckman et al. 1990), with large-scale superbubbles expanding by tens of kpc toward north–west and south–east (Veilleux et al. 2003; Yoshida et al. 2016). The Hα emission from the ionized wind shows a close spatial correspondence with the soft X-ray continuum, suggesting the presence of gas cooling out of a shocked medium (Nardini et al. 2013). Furthermore, Wang et al. (2014) detected a diffuse component in hard-X-ray continuum and Fe xxv line emission, tracing T ∼ 7 × 10⁷ K gas between the two AGNs (north–west of the southern nucleus, similar to the [C i](1–0) blue wing in Figure 1(d)), as well as in kpc-scale structures that are remarkably coincident with both the strong NIR H₂ emission (van der Werf et al. 1993; Max et al. 2005) and the Hα filaments.

At radio wavelengths, Colbert et al. (1994) reported the detection of non-thermal continuum emission with a steep spectrum extending by several kpc in an arclike structure west of the AGNs, later confirmed also by Baan et al. (2007), together with a possibly similar feature on the eastern side. This structure lacks a clear spatial correspondence with optical or NIR starlight (which excludes a starburst origin), and its complex morphology suggests a connection with the Hα outflow. Theoretically, the association between an AGN-driven wind and extended non-thermal radiation (due to relativistic electrons accelerated by the forward shock) has been predicted by Nims et al. (2015). Observationally, the rough alignment of the western arclike structure discovered by Colbert et al. (1994) with the molecular outflow studied in this work would also support this hypothesis, although the current data do not allow us to probe the presence of H₂ outflowing gas at the exact position of the radio emission. The total radio power of this arclike feature is comparable with that of extended radio structures previously observed in radio-quiet AGNs with prominent outflows (Morganti et al. 2016). Future facilities like the Cherenkov Telescope Array may reveal the γ-ray counterpart of the non-thermal emission, as expected for an AGN outflow shock (Lamastra et al. 2017).

In summary, all these multi-wavelength observational evidences point to a radiative-mode AGN feedback mechanism (Faucher-Giguère & Quataert 2012; Nims et al. 2015). However, at the same time, it is difficult to reconcile a classic model with an H₂ outflow whose emission does not peak on either of the two AGNs. A complex interplay of stellar and AGN feedback must be at work in this source (see also Müller-Sánchez et al. 2018), and we cannot exclude the additional contribution from compact radio-jets (Gallimore & Beswick 2004), which may be accelerating part of the cold material (Mukherjee et al. 2016). Finally, positive feedback may also be at work in NGC 6240. There is indeed a striking correspondence between (1) the morphology of the approaching side of the outflow north–west of the southern AGN (Figure 1(d)), (2) a peak of dust extinction, and (3) a stellar population with unusually large stellar σ_v and blueshifted velocities (Engel et al. 2010). The latter, according to Engel et al. (2010), may have been formed recently as a result of crushing of molecular clouds, which could be related to the observed outflow event (Zubovas & King 2014; Maiolino et al. 2017), provided these stars are not older than a few Myr.