ALMA OBSERVATIONS OF A CANDIDATE MOLECULAR OUTFLOW IN AN OBSCURED QUASAR

SDSS J1356+1026 is a merging system at z = 0.123 classified as a luminous obscured (Type 2) quasar with a bolometric luminosity L_bol ≈ 10⁴⁶ erg s⁻¹ inferred from the [O iii] luminosity, ${L_{\mathrm{[O\, \scriptsize{III}]}}}$ = 10^42.77 erg s⁻¹ (Greene et al. 2009; Liu et al. 2009). It is also an ultra-luminous infrared galaxy (ULIRG; L_FIR = 2.68 ± 0.53 × 10⁴⁵ erg s⁻¹). The radio continuum is unresolved with a luminosity density of νL_ν = 3.4 × 10⁴⁰ erg s⁻¹ from the Very Large Aray Faint Images of the Radio Sky at Twenty-centimeters (FIRST) survey at 1.4 GHz (1.6 GHz rest frame), indicating that it is a radio-quiet quasar (Greene et al. 2012). The coordinate of SDSS J1356+1026 is (13:56:46.10, +10:26:09.09).

The two merging galaxies, which we refer to as the northern (N) and southern (S) nuclei, are separated by 2.5 kpc (11). The N nucleus is detected at 2–10 keV and is the primary AGN host in the system (Greene et al. 2014). Both the total molecular mass (Section 3.3) and the r-band light (Greene et al. 2009) are dominated by the northern nucleus, with a ratio of ∼4: 1 (N:S). The progenitor of the northern galaxy is likely a moderately massive early-type galaxy, with a stellar spectrum that is dominated by an old stellar population and a stellar velocity dispersion σ_* = 206 ± 36 km s⁻¹ (Greene et al. 2009), corresponding to a stellar mass M_* ≈ 10¹¹ M_☉ (Hyde & Bernardi 2009), and the BH mass is estimated to be M_BH ≈ 3 × 10⁸ M_☉ (±0.38 dex, McConnell & Ma 2013). The corresponding Eddington luminosity is L_Edd =3.1 × 10⁴⁶ erg s⁻¹ (±0.38 dex), which yields an Eddington ratio range of 0.1–1, energetic enough to drive a hot wind (e.g., Veilleux et al. 2005).

The most spectacular feature of this system is an [O iii]-emitting bubble extending ∼10 kpc to the south of the nuclei (green region in Figure 1), with a weaker symmetric counterpart to the north. Long-slit spectroscopy along the outflow reveals the distinctive double-peaked spectrum of an expanding shell. A simple geometric model yields a deprojected velocity v ≈ 1000 km s⁻¹, a dynamical time t_dyn ≈ 10⁷ yr, and a kinetic luminosity $\dot{E} \approx 10^{44-45}$ erg s⁻¹ (Greene et al. 2012). As the radio emission is compact and faint (Sections 3.1 and 5.2.2), and the star formation rate (SFR) is low (Sections 3.3 and 3.4), this ionized outflow provides a strong case for a quasar-driven wind (Greene et al. 2014).

Zoom In Zoom Out Reset image size

The ALMA CO (1–0) and CO (3–2) observations were conducted during Cycle 0 and Cycle 1 under project codes 2011.0.00652.S and 2012.1.00797.S, respectively. The CO (1–0), at sky frequency 102.6 GHz, was observed in Band 3 in two blocks on 2012 May 9 and July 30, using 16 and 23 12 m antennae for 68 and 63 minutes, respectively (27 and 24 minutes on-source). The CO (3–2) at 307.9 GHz was observed in Band 7 on 2013 July 6 with 27 12 m antennae in 24 minutes (19 minutes on-source). In the following, CO (1–0) properties are followed by CO (3–2) in parenthesis. Both lines use a single pointing covering the whole system with a field of view of 62'' (21''), and dual polarizations. They use the same spectral set-up with a total bandwidth of 2 GHz divided into 128 channels, each 15.625 MHz wide, corresponding to a channel width of 46 km s⁻¹ (15 km s⁻¹). The CO (1–0) data has an additional spectral window at 93 GHz. The channels are Hanning-smoothed within the correlator. QSO 3C279 (J1516+0015), Mars (Titan), and J1415+133 (J1347+1217) were used for the bandpass, flux, and phase calibrators, respectively. The phase calibrator is observed for 30 (90) seconds every 11 (9) minutes. The absolute flux calibration accuracy is 5% for CO (1–0) (Band 3) and 10% for CO (3–2) (Band 7).

The data calibration was carried out by the ALMA team using CASA version 3.3.0 for CO (1–0) and 4.1.0 for CO (3–2). We further used CASA version 4.1.0 to apply continuum subtraction [CO (1–0) only], heliocentric correction, imaging, and cleaning. For CO (1–0), the continuum level was estimated from line-free channels (36/43 channels on the blue/red side) and subtracted in the uv-plane. Continuum subtraction was not applied to the CO (3–2) data cube as there are not enough line-free channels to determine the continuum level. However, we extrapolate from the 100 GHz and 1.4 GHz continuum flux densities (spectral index α = −0.86, Section 3.4) to estimate a 300 GHz flux density of 0.58 mJy. This extrapolated continuum level is lower than the noise level in each CO (3–2) channel and is less than 10% of the emission line flux density at the N nucleus. The CO (3–2) channel is further binned by two to increase the signal-to-noise ratio in each channel, giving a final channel width of 30 km s⁻¹, while the CO (1–0) channel width is unchanged (46 km s⁻¹). The imaging processing was carried out with the CASA task clean. We used the Briggs visibility weighting with a robustness parameter of 0.5 (Briggs 1995), which is a compromise between maximum resolution and maximum sensitivity. The beam size is θ_beam = 19 × 13 (035 × 029) with a P.A. = −62° (− 60°), and the maximum resolvable scale is θ_MRS = 29'' (10''). All images are cleaned to the 3σ level in each channel. The resulting rms noise level in each channel is 0.37 mJy beam⁻¹ (0.77 mJy beam ⁻¹).

In order to compare the molecular and stellar components of the galaxies, we match the ALMA data cube to the Sloan Digital Sky Survey (SDSS) spectrum in velocity and the Hubble Space Telescope (HST)/WFC3 F814W (I-band) image in position (J. M. Comerford et al., in preparation). We use the SDSS DR7 spectrum (Abazajian et al. 2009), which was taken in a 3'' aperture with a spectral resolution of FWHM ≈ 150 km s⁻¹ centered on the N nucleus. The HST/WFC3 F814W image was observed on 2012 May 19 with an integration time of 900 s and resolution of 007.

To perform positional alignment, we identify another galaxy SDSS J135646.50+102553.6 that appears in the ALMA, SDSS, and HST images, and two other fainter galaxies in both SDSS and HST. We found that while the ALMA and SDSS positions are well-matched, there is a 06 offset from the HST position. After applying this offset to the HST image, the ALMA (continuum) and HST coordinates of the northern nucleus are aligned well within 004, much smaller than the ALMA CO (3–2) beam size (035 × 029). Rotation and stretching of the HST image with respect to the SDSS coordinates are also constrained to be within 04 rotation and 0.5% stretching.

The velocities in this paper are with respect to the rest frame of the stellar absorption features of N galaxy. We define the systemic velocity to be the best-fit redshift (z = 0.1231) from the SDSS spectrum, which is dominated by the light of the N nucleus, as the N nucleus is brighter and the fiber is centered on it.

Both the stellar absorption and AGN emission lines are fitted by the SDSS template. A redshift warning flag "many outliers" is raised, indicating that the template cannot capture the multi-component emission lines, but this is not a major concern for redshift accuracy. As a sanity check, we refit the stellar continuum using Bruzual & Charlot (2003) stellar population synthesis templates and the best-fit redshift differs from that of the SDSS by −30 km s⁻¹. We therefore adopt a redshift error of δz = ±0.0001 (± 30 km s⁻¹). This stellar continuum fitting also reveals that the light in the N nucleus is dominated by old stellar populations, as discussed in Greene et al. (2009).

The radio 100 GHz continuum is measured from the line-free channels in two spectral windows at 93/103 GHz (rest-frame 104/115 GHz) in the CO (1–0) data. The moment-0 map (Figure 4) shows an unresolved point source with a beam size of 19 × 13 (4.2 × 2.9 kpc) at the N nucleus with a flux density of 1.51 ± 0.07 ± 0.07 mJy. The first error is from the rms noise, and the second is the 5% calibration error. The spectral energy distribution (SED) of SDSS J1356+1026 including this 100 GHz flux density is discussed in Section 3.4.

Zoom In Zoom Out Reset image size

We extract the CO(1–0)/(3–2) spectra (flux densities F_ν, Figure 3), and fluxes (F = F_νΔv, Table 1) for the three components (the N/S nuclei and the W arm) from the cleaned data. In addition to the rms noise errors presented in Figure 3 and Table 1, there is a 5% and 10% flux calibration error for the CO (1–0) and CO (3–2) data, respectively. We start with the CO (3–2) data as it has higher resolution to spatially separate the components. We use simple aperture photometry with apertures at least twice as large as the beam in order to enclose ${\mathrel {{\rlap{{{\sim }}}{{>}}}}}95\%$ of the flux. As the N and W apertures are adjacent to each other, we expect flux contamination between these two components at a level of 10% or less for CO (3–2). The S aperture is well-separated from the other two sources and does not overlap with side-lobes of the N nucleus.

Assuming Gaussian noise correlated on the scale of the beam, we estimate the errors in the spectra F_ν based on the image noise level rms(I_ν), taken from large emission-free regions, the beam $\mathrm{B}(\boldsymbol{r})$, a two-dimensional Gaussian with peak value 1, and the aperture $T(\boldsymbol{r}$) = 1 inside the aperture and zero everywhere else, adopting the following equation:

$$\begin{equation} \mathrm{rms}(F_{\nu })=\mathrm{rms}(I_{\nu })\sqrt{\int \int \mathrm{B}(\boldsymbol{r}_2-\boldsymbol{r}_1)T(\boldsymbol{r}_1)T(\boldsymbol{r}_2) {d}^2\boldsymbol{r}_1 {d}^2\boldsymbol{r}_2 }. \end{equation} \tag{ 1 }$$

To obtain the fluxes F, we integrate the spectra F_ν over a velocity range of −300 km s⁻¹ <v < 500 km s⁻¹ for the N nucleus and W arm, which is chosen as the velocity width where the CO (3–2) flux density is detected with greater than 2σ significance from the N nucleus. A velocity range of 300 km s⁻¹ <v < 500 km s⁻¹ is used for the high-velocity component of the N nucleus. For the S nucleus, although there is no CO (3–2) detection, we use the velocity range −80 < v < 50 which covers the >1σ CO (1–0) emission. While estimating the errors in the fluxes, we take into account the correlation between adjacent channels, which has a Pearson correlation coefficient of r = 0.375, due to Hanning smoothing during the observation and 2:1 binning in the image processing.

Estimating the CO (1–0) fluxes is more complicated, as the sources, especially the N nucleus and the W arm, are spatially blended. The CO (1–0) spectrum of the N nucleus clearly shows the narrow peak of the W arm at −100 km s⁻¹ (Figure 3). However, we can still extract the uncontaminated N nucleus fluxes making use of their distinct spectral shape. We first estimate the N nucleus flux in the 100–500 km s⁻¹ channels that are free from W arm emission by fitting a model of the beam to the stacked moment-0 map of these channels. This point-source assumption should be valid as the CO (1–0) beam is much larger than the CO (3–2) N nucleus size. We then recover the total N nucleus CO (1–0) flux in the entire velocity range (−300 to 500 km s⁻¹) assuming that the CO (1–0) spectrum has a similar spectral shape to the CO (3–2) spectrum.

The W arm CO (1–0) flux is then estimated as the difference between the total flux including both sources and the decomposed N nucleus flux. The CO (1–0) flux in the high-velocity component of the N nucleus is estimated by the same point source fitting procedure over the velocity range of 300–500 km s⁻¹, as is the S nucleus over a velocity range of −80 to 50 km s⁻¹, where the emission is seen. Correlations between adjacent channels are determined from emission-free regions and are taken into account for the flux error estimation. Using the estimated fluxes of the components in both lines, we express the line ratio as $L^{\prime }_{\mathrm{CO (3{-}2)}}/L^{\prime }_{\mathrm{CO (1{-}0)}}$ in Table 2, where L' is in units of K km s⁻¹ pc², proportional to the brightness temperature.

From the measured CO (1–0) luminosity we can estimate the molecular gas mass (e.g., Bolatto et al. 2013). The masses are listed in Table 2. The CO-to-H₂ conversion factor X_CO is defined as

$$\begin{eqnarray} &&X_{\mathrm{CO}} = N_{\mathrm{mol}}/W(\mathrm{CO (1{-}0)}), \end{eqnarray} \tag{ 2 }$$

where N_mol is the molecular column density in units of cm⁻², and W(CO(1–0)) is the CO (1–0) brightness in units of K km s⁻¹. Although X_CO is roughly constant in normal Galactic molecular clouds (Solomon et al. 1987), it is found to be lower in ULIRG/starburst galaxies by a factor of two to five with large scatter (Downes & Solomon 1998). Likely because the molecular clouds blend together due to tidal forces in the dense environment of the ULIRG nucleus, the CO emission emerges from a diffuse volume-filling medium rather than discrete self-gravitating molecular clouds.

SDSS J1356+1026 is identified as a ULIRG and has disturbed gas dynamics, so we adopt an X_CO value for ULIRGs X_CO = 0.4 × 10²⁰ cm⁻² (K km s⁻¹)⁻¹ with an error of 0.5 dex (Downes & Solomon 1998; Bolatto et al. 2013). This corresponds to α_CO = 0.86 M_☉ (K km s⁻¹ pc²)⁻¹, matching the value used in other ULIRG/AGN CO studies (e.g., Cicone et al. 2014). This factor of three uncertainty in the X_CO factor dominates the errors in the mass estimates. We find a total molecular mass of $9^{+19}_{-6}\times 10^{8} \,M_{\odot }$. Half of the molecular mass is in the extended W arm ($5^{+11}_{-3}\times 10^{8} \,M_{\odot }$), while the N nucleus shares a third ($3^{+6}_{-2}\times 10^{8} \,M_{\odot }$). As CO (1–0) is only marginally detected at the S nucleus, a conservative mass limit of <3.8 × 10⁸ M_☉ is estimated from the 3σ detection limit and the uncertainty in X_CO.

Since molecular gas is the fuel for star formation, we can estimate the star formation that can be supported by the molecular content in SDSS J1356+1026. Using the Schmidt–Kennicutt law (Kennicutt 1998), we find the SFR at the N nucleus to be SFR = 1.2 M_☉ yr⁻¹ (±0.5 dex) assuming that all of the molecular gas at the N nucleus is distributed in a disk with a radius of 300 pc, half of the beam deconvolved source FWHM in CO (3–2). As the morphology of the molecular gas in the W arm and the S nucleus are uncertain, their SFR is poorly constrained. A conservatively high estimate of the entire system can be calculated by putting all the molecular gas, including that in the diffuse W arm, in a disk with a radius of 300 pc, which gives 5 M_☉ yr⁻¹. Taking the factor of three X_CO uncertainty into account, an absolute upper limit is placed at 16 M_☉ yr⁻¹. However, the SFR inferred from the CO luminosity is indirect. Whether the assumed X_CO factor and/or the Schmidt–Kennicutt law apply in this environment is uncertain. In the following section, we constrain the SFR directly from the infrared SED decomposition.

We use the infrared SED decomposition to constrain the SFR (Figure 5 and Table 3). It is usually thought that the dust close to the AGN is heated to much higher temperature than the diffuse dust heated by starlight. The difference in the dust temperatures is at the core of the SED-fitting method to distinguish between AGN-dominated and star formation-dominated sources. Greene et al. (2012) compiled a UV to radio SED of SDSS J1356+1026 making use of the Two Micron All Sky Survey (Skrutskie et al. 2006), IRAS (Neugebauer et al. 1984), FIRST (Becker et al. 1995), and NVSS (Condon et al. 1998) data. We further include data from the Wide-field Infrared Survey Explorer (WISE; Wright et al. 2010) and the Herschel Space Observatory (Pilbratt et al. 2010) PACS (Poglitsch et al. 2010) and SPIRE (Griffin et al. 2010) data from A. Petric et al. (in preparation), which extends to 500 μm to cover the Rayleigh–Jeans tail of the dust emission.

Zoom In Zoom Out Reset image size

Following the decomposition procedures of Kirkpatrick et al. (2012), we fit a two-temperature modified blackbody to the seven-point SED. This SED of 22 μm WISE, 60 μm IRAS, and 70–350 μm Herschel is sensitive to the temperature of the dust emission. There are three free parameters in the model: the temperature T_warm and luminosity L_warm of the warm component and the luminosity of the cold component L_cold. The temperature of the cold dust is fixed at T_cold = 35 K, as found in a wide range of systems by Kirkpatrick et al. (2012). Even if we fit the temperature of the cold component, we recover T_cold = 35 K. The emissivity of the dust is assumed to be a power-law function of frequency with index β = 1.5. All photometry errors are assumed to be 10%, which is representative of the systematic errors. The best-fit model (black line, Figure 5) matches the data well with a reduced χ² of 1. The best-fit temperature of the warm component is high T_warm = 81 K, and the luminosity is L_warm = 2.6 × 10⁴⁵ erg s⁻¹, much higher than that of the cold component L_cold = 4.3 × 10⁴⁴ erg s⁻¹.

The luminosity-weighted temperature T_eff = 75 K is higher than that seen in Kirkpatrick et al. (2012) even in those of their objects that have an AGN-dominated SED (T_eff = 65 K). Therefore, as the warm component is consistent with being AGN-heated, the AGN is the major heating source for the bulk of the dust emission. However, star-formation cannot be ruled out as the heating source for the 35 K cold component. We therefore use the luminosity of this 35 K component as an upper limit on the luminosity of dust emission associated with star formation, and use calibrations from Bell (2003) to calculate an upper limit on the SFR of 21 M_☉ yr⁻¹. If the SFR were much higher, at the same T_cold, the dust emission would exceed the measured 250 μm to 350 μm flux. This SFR upper limit is insensitive to the temperature and luminosity of the warm component because the 250 μm and 350 μm luminosities are dominated by the cold component. Harrison et al. (2014) estimated a higher SFR of $63^{+7}_{-17}$ M_☉ yr⁻¹, inconsistent with our limit. This high SFR is inferred from the WISE and IRAS data alone covering 4.6 μm–100 μm. In this case, the Herschel photometry at 250 μm and 350 μm are critical for the SED decomposition and SFR estimation.

We plot our measured 100 GHz flux (Section 3.1) in the SED (right panel of Figure 5). The flux at 100 GHz is much higher than the extrapolated Rayleigh–Jeans tail of the dust emission, and thus must have a different origin. We find a spectral slope of α = −0.86 (F_ν∝ν^α) between the 100 GHz and the 1.4 GHz fluxes, similar to the typical spectral indices α ≈ −0.7 seen from synchrotron emission for radio-loud AGNs (Zakamska et al. 2004).

Mid-infrared (5–25 μm) luminosity has also been used as a bolometric indicator for the AGN. We use the WISE photometry to infer L_bol following similar procedures as Liu et al. (2013b). As pointed out by Liu et al. (2013b), type 2 quasars have significantly redder mid-infrared colors than type 1 unobscured quasars, possibly due to dust reddening. This red color is also seen in SDSS J1356+1026. Therefore, to avoid attenuation from the dust, we choose the reddest WISE band at 22 μm (19.5 μm rest frame) to apply the bolometric correction of L_bol = νL_ν × (11 ± 5) from Richards et al. (2006), which gives L_bol ≈ 1.1 ± 0.5 × 10⁴⁶ erg s⁻¹. However, this L_bol might still be underestimated due to the strong extinction. We further use the 30 μm luminosity density extrapolated from the 3.4, 4.6, 12, and 22 μm data assuming a power-law spectral shape to set a conservative upper limit of L_bol < 2.3 × 10⁴⁶ erg s⁻¹ with a bolometric correction factor of 13 (Liu et al. 2013b). The L_bol estimate from the mid-infrared data is similar to that inferred from the [O iii] luminosity, L_bol ≈ 10⁴⁶ erg s⁻¹ (Greene et al. 2012), using a L_bol–${L_{\mathrm{[O\, \scriptsize{III}]}}}$ relation (Liu et al. 2009) calibrated with type 1 AGNs.

There is strong CO (1–0) and CO (3–2) emission at the N nucleus with a wide velocity distribution from −300 km s⁻¹ to 500 km s⁻¹ (Figure 3). From the CO (1–0) flux, the molecular gas mass is estimated to be $3^{+6}_{-2} \times 10^8$ M_☉, which can fuel star formation at an approximate rate of ∼1 M_☉ yr⁻¹, if the molecular gas is in a disk of radius 300 pc. The N nucleus has a CO line ratio $L^{\prime }_{\mathrm{CO (3{-}2)}}/L^{\prime }_{\mathrm{CO (1{-}0)}} = 1.18\pm 0.73$. As various factors affect the line ratio, including the gas temperature and density, we are unable to infer the physical conditions of the gas from this one ratio. Instead, we empirically compare our measured line ratio to other objects using the compilation in Carilli & Walter (2013), including quasars ($L^{\prime }_{\mathrm{CO (3{-}2)}}/L^{\prime }_{\mathrm{CO (1{-}0)}}\approx$ 1.0), submillimeter galaxies (∼0.7), normal star-forming galaxies (∼0.4–0.6, Bauermeister et al. 2013a), and the Milky Way (∼0.3). The N nucleus has a relatively high line ratio similar to that of a quasar.

The CO (3–2) emission is marginally resolved (FWHM = 045 > Beam = 035) with a decomposed FWHM of 0.6 ± 0.1 kpc. This spatial resolution allows us to analyze the molecular kinematics at the N nucleus. The channel maps (Figure 6) show that the emission moves toward the southeast as the velocity increases until about 300 km s⁻¹. Above 300 km s⁻¹, there is high-velocity emission extending to 500 km s⁻¹, but the position does not align with the low-velocity gas.

Zoom In Zoom Out Reset image size

In each channel, we fit the emission with a two-dimensional Gaussian profile, and derive the best-fit x and y positions (Figure 6). We fit a linear model ${\boldsymbol {x}}={\boldsymbol {a}}v+{\boldsymbol {b}}$ to the best-fit positions in the velocity range of ±300 km s⁻¹ to yield a P.A. =127°. We then extract a position–velocity diagram (hereafter pv diagram, Figure 7) along this P.A. The pv diagram is extracted over a width of 05, which encloses >95% of the flux. It can be seen in the pv diagram that the emission within ±300 km s⁻¹ roughly follows a linear position–velocity relation that is symmetric about the systemic velocity and has a best-fit gradient of 1077 km s⁻¹ kpc⁻¹. The stellar velocity dispersion at the N nucleus (within a 1'' aperture) is σ_* = 206 ± 36 km s⁻¹ (Greene et al. 2009), which would correspond to a circular velocity of $V_c = \sqrt{2}\sigma _* = 291 \pm 50$ km s⁻¹ for an isothermal sphere. Therefore, the linear velocity gradient within ±300 km s⁻¹ is consistent with a rotating disk. This velocity gradient corresponds to a constant dynamical mass density of ρ_dyn ≈ 63 M_☉ pc⁻³, which gives a dynamical mass within a radius of 0.25 kpc (0.3 kpc) of M_dyn ≈ 4 × 10⁹ M_☉ (7 × 10⁹ M_☉). The molecular gas mass of $M_{\rm {mol}} = 3^{+6}_{-2} \times 10^8$ M_☉ inferred from the CO emission accounts for 1%–20% of the dynamical mass in the N nucleus.

Zoom In Zoom Out Reset image size

Beyond this rotating component, the redshifted high-velocity feature (300 < v < 500 km s⁻¹) deviates from the linear velocity gradient and has a velocity too high to be explained by rotation alone. The bottom of Figure 6 and the left panel of Figure 7 show that the position of this feature has a slight offset (002) to the northeast from the disk plane. There is no counterpart on the blueshifted side. The CO (3–2) emission is detected with 9σ significance (0.99 ± 0.10 Jy km s⁻¹), while there is only a marginal detection 0.06 ± 0.08 Jy km s⁻¹ for CO (1–0). Therefore, the line ratio $L^{\prime }_{\mathrm{CO (3{-}2)}}/L^{\prime }_{\mathrm{CO(1{-}0)}}=1.8\pm 2.5$ is poorly constrained, but hints at a high excitation level. The 3σ detection limit of CO (1–0) places a lower limit on the line ratio $L^{\prime }_{\mathrm{CO (3{-}2)}}/L^{\prime }_{\mathrm{CO(1{-}0)}} > 0.45$. We use the CO (3–2) flux to infer a molecular mass $M_{\mathrm{mol}} =7^{+15}_{-5}\times 10^7$ M_☉, assuming a line ratio of $L^{\prime }_{\mathrm{CO (3{-}2)}}/L^{\prime }_{\mathrm{CO(1{-}0)}}=1$ as observed in the bulk of the N nucleus as well as in other AGNs. This high-velocity component accounts for a significant fraction (∼20%) of the molecular mass in the N nucleus, under the assumption of constant line ratio and X_CO factor. We discuss the interpretation of this component in Section 5.2.

The extended emission to the west of the N nucleus is 2'' (4 kpc) long, and has a blueshifted narrow line at −150 km s⁻¹ with a FWHM of 200 km s⁻¹ (Figure 3). The CO (1–0) flux is 0.82 ± 0.34 Jy km s⁻¹ once the overlap with the N nucleus is accounted for, and the CO (3–2) flux is 2.55 ± 0.40 Jy km s⁻¹. Adopting a ULIRG X_CO ratio, the W arm contains $5^{+11}_{-3} \times 10^8$ M_☉ of molecular gas, about half of the total mass in the system. The actual mass could be even higher, up to 2.5 × 10⁹ M_☉, if the molecular gas is in self-gravitating clouds and the higher Milky Way X_CO ratio applies. It has a moderate line ratio of $L^{\prime }_{\mathrm{CO (3{-}2)}}/L^{\prime }_{\mathrm{CO (1{-}0)}} = 0.35\pm 0.20$, which is lower than the N nucleus and is similar to star-forming galaxies (0.4–0.6, Bauermeister et al. 2013a; Carilli & Walter 2013).

The W arm is likely a tidal feature (Section 5.1). It is only slightly offset from the position of an extended stellar plume to the north west of the N nucleus (see the CO (1–0) map in Figure 1). Decoupling between stellar and gas components in merging systems is not only predicted by simulations (e.g., Barnes & Hernquist 1996) but also commonly seen in observations of tidal tails (e.g., Hibbard et al. 2000), as the gas is subject to dissipational processes that can affect its trajectory. Furthermore, the cool and coherent kinematics of the W arm, as seen from its narrow line width (Figure 3), across its full length of 4 kpc, makes it hard to explain by either outflow or inflow.

At the location of the southern nucleus identified from the HST images, there is only a ∼2σ CO (1–0) detection of 0.13 ± 0.06 Jy km s⁻¹, and no detection for CO (3–2). The 3σ detection limit of CO (1–0) translates into a molecular gas mass of 1.2 × 10⁸ M_☉, about 10% of the total molecular gas mass in the system. Taking the uncertainty in the X_CO factor into account, we can place an upper limit of M_mol < 3.8 × 10⁸ M_☉. Therefore, the S nucleus is not a major reservoir of molecules in SDSS J1356+1026. Although the sensitivity is not adequate to constrain the detailed line profiles, the CO (1–0) marginal detection suggests a narrow blueshifted line centered on −30 km s⁻¹ (Figure 3). The much narrower CO (1–0) line width compared to the N nucleus can result from the shallower gravitational potential or the absence of gas outflow/inflow in the S galaxy. The $L^{\prime }_{\mathrm{CO (3{-}2)}}/L^{\prime }_{\mathrm{CO (1{-}0)}}$ line ratio is constrained to be 0.10 ± 0.11, much lower than the N nucleus and the W arm.

As we describe in Section 1.1, SDSS J1356+1026 contains an extended ionized outflow discovered by Greene et al. (2012). However, there is no detection of molecular gas at the corresponding location. We identify the region of ionized line emission from the elongated luminous structure in the HST F814W image (green region in Figure 1). The integrated CO (1–0) flux of −32 ± 102 mJy km s⁻¹ places a molecular gas mass upper limit (95% c.l.) of M_mol < 8.6 × 10⁷ M_☉, while Greene et al. (2012) give a lower limit on the ionized gas mass M_ion > 5 × 10⁷ M_☉. Therefore, the molecular gas does not dominate the mass of the extended ionized outflow.

Despite the non-detection of CO in the ionized outflow, we note a tentative detection of a double-peaked CO (3–2) emission line at an optically bright spot in the HST F814W image about 4'' south of the N nucleus where the [O iii] expanding-shell signature starts to emerge (the "bubble base" shown in magenta in Figure 1). The two CO (3–2) peaks, each having 3σ–4σ significance, are at −600 and +400 km s⁻¹, respectively. Despite the intriguing location and velocity of this component, its confirmation requires further observations.

We first discuss the origin and evolution of the molecular gas in the system. A significant amount of molecular gas is detected in the N nucleus and the W arm with a total mass of $9^{+19}_{-6}\times 10^8$ M_☉ (Section 3). However, the stellar absorption spectrum of the N nucleus is dominated by an old stellar population (Greene et al. 2009), and based on the velocity dispersion we suspect that the galaxy started as an early-type with little cold gas. For local elliptical/S0 galaxies like SDSS J1356+1026 with σ_* > 200 km s⁻¹, Young et al. (2011) find that only 1 out of 30 galaxies has a molecular gas mass as high as 8 × 10⁸ M_☉, while most of the ellipticals are below the CO detection limit of M_mol < 10⁸ M_☉. Therefore, the exceptionally high molecular content in the N galaxy was likely accreted externally, perhaps from the S galaxy. The current optical spectroscopy is not adequate to identify the stellar population of the S galaxy. However, if it was a late-type galaxy, with a stellar mass of about M_* ≈ 3 × 10¹⁰ M_☉, it could bring in plenty of molecular gas, as the typical molecular to stellar mass ratio is about 5% for local spirals (Leroy et al. 2009; Bauermeister et al. 2013b). Gas transfer from disk companions to elliptical galaxies has been observed (Schweizer 1987). Furthermore, some early-type galaxies have cold gas components that are kinematically decoupled from the stellar components, suggesting an external origin of the cold gas (Knapp 1987; Davis et al. 2013). Therefore, a primary elliptical/satellite disk galaxy merger is a likely scenario to explain the molecular gas content in SDSS J1356+1026.

In simulations, when the gas collides during a merger, its orbital energy and angular momentum are dissipated, and thus the gas tunnels into the center to form a nuclear disk and triggers AGN or star formation activity (e.g., Hernquist 1989; Barnes 2002). In particular, if it is a primary elliptical and satellite disk galaxy merger, in some cases the gas in the disk galaxy is disrupted and sinks into the center of the primary elliptical (Weil & Hernquist 1993). These nuclear gas disks, growing from infalling tidal tails, usually contain about half of the cool gas in the system and have peculiar kinematics relative to the stars. While these simulations do not treat the molecular gas phase explicitly, we expect the nuclear disks to have a large molecular fraction because of the dense environment. Observationally, molecular nuclear gas disks and rings are indeed common in infrared luminous mergers (Downes & Solomon 1998). These compact nuclear disks have radii ∼0.5 kpc, as opposed to ∼1 kpc in normal ellipticals (Davis et al. 2013), and can be as massive as 10⁹ M_☉, similar to the compact nuclear disk of SDSS J1356+1026. The energetic nuclear activity of SDSS J1356+1026 could be a direct result of the nuclear inflow induced by the interaction between the two merging galaxies (Barnes & Hernquist 1996).

Observationally, although extended molecular tails can survive after being stripped from the galaxies (e.g., through ram pressure stripping, Dasyra et al. 2012; Jachym et al. 2014), it is not common to find a significant amount of molecular gas in tidal tails. Aalto et al. (2001) found 4% (M_mol ≈ 8 × 10⁷ M_☉) of the molecular gas along a 10 kpc long optical tidal tail in NGC 4194, an elliptical/disk merger remnant. The ULIRG merger Mrk 273 also shows a 5 kpc extended CO streamer with mass M_mol ≈ 1 × 10⁸ M_☉, accounting for a few percent of the CO flux (Downes & Solomon 1998). One of the nuclei of Mrk 273 has no molecular gas, so it was possibly a gas poor galaxy or the gas was depleted by transfer (Yun & Scoville 1995). However, none of these examples contain as much molecular gas as in the W arm of SDSS J1356+1026 (5 × 10⁸ M_☉, half of the total molecular mass).

Both SDSS J1356+1026 and NGC 4194 are probably elliptical/disk minor mergers. We suspect that the formation of molecular tidal features may require particular progenitor properties and orbital configurations. For example, in the primary elliptical/satellite disk merger, the elliptical with its concentrated potential may disrupt the molecular disk in the satellite galaxy. At the same time, the absence of gas in the elliptical can reduce the shock dissociation of molecular gas during the encounter and thus preserve gas in the molecular phase in tidal features.

In Section 4.1, we find a high-velocity CO (3–2) feature at the N nucleus with a radius of ∼300 pc. This feature ranges in velocity from 300 to 500 km s⁻¹, which is too high to be part of the rotating disk. Its mass is estimated to be M ≈ 7 × 10⁷ M_☉, a quarter of the mass at the N nucleus. Also, it has complex kinematics that deviate from rotation and has a slight offset with respect to the disk plane (Figures 6 and 7). There are various possible interpretations for this high-velocity component. It could be inflowing or outflowing, or could result from tidal forces or feedback mechanisms. We discuss two of the most plausible scenarios.

First, the high-velocity feature may represent material accreting onto the disk from an external source, perhaps the W tidal arm. The projected angular momentum of the two features are parallel to each other, although they lie on opposite sides of the disk. Furthermore, the angular momentum of the high-velocity gas is smaller than the W arm, as expected for gas that is accreting onto a nuclear disk.

The second possibility is a molecular outflow. Herschel OH line observations have found molecular outflow to be a common phenomenon in ULIRG/AGN like SDSS J1356+1026 (Veilleux et al. 2013). The high-velocity feature seen in SDSS J1356+1026 meets one of the outflow criterion used by the CO interferometry study by Cicone et al. (2014), in that it has a velocity above 300 km s⁻¹ and deviates from the rotation pattern, making it a likely outflow candidate. However, CO confirmed molecular outflows in luminous ULIRG/AGN are typically more massive and larger in size. The six molecular outflows in luminous ULIRG/AGN (L_bol > 10^44.5) discussed by Cicone et al. (2014), including IRAS F08572+3915, IRAS F10565+2448, IRAS 23365+3604, Mrk 273, Mrk 231, and NGC 6240, have molecular outflow masses M ≈ 3 × 10⁸, sizes r = 0.6–1.2 kpc, and velocities v = 400–1200 km s⁻¹. This discrepancy in size could be partly due to the high angular resolution of the ALMA observations. Structures of size 300 pc may not have been resolved by previous CO studies. The fact that this compact size is measured by CO (3–2), a denser tracer than CO (1–0), may also contribute to the discrepancy. In fact, very compact molecular outflows on the scale of ∼100 pc have been discovered by OH observations combined with radiative-transfer modeling (Sturm et al. 2011; González-Alfonso et al. 2013).

With the current observations we are unable to completely confirm or rule out either scenario. Given the ubiquity of nuclear outflows observed with Herschel on similar scales, in what follows we elaborate further on the outflow scenario and discuss the implications of the possible compact molecular outflow in this source.

5.2.1. Properties of the Molecular Outflow

5.2.2. AGN or Star Formation Driven Outflow?

5.2.3. Comparison with the Ionized Outflow