Molecular Oxygen in the Nearest QSO Mrk 231

The observations were carried out from 2015 August 20 to 23, with the IRAM 30 m telescope in good weather conditions with precipitable water vapor less than 4 mm (Project id: 068-15, PI: Junzhi Wang). The Eight MIxer Receiver (EMIR) with dual polarizations, fast Fourier transform spectrometer backend, and standard wobbler switching mode at 0.5 Hz with ±120'' offset beam throw, were used toward Mrk 231 pointing at R.A. = 12:56:14.2 decl. = +56:52:25.0 (J2000). The beam size of the IRAM 30 m telescope at the observing frequency (∼114 GHz) is about 23''. In order to verify that the signal is from the sky rather than being radio-frequency interference at the IF frequency or from the backend, four local oscillator tuning setups were used during the observations, which are 103.77, 103.87, 103.57, and 103.07 GHz for August 20–23, respectively. We dumped the data every 1.7 minutes as one scan, and did one calibration observation every six scans. The CLASS package GILDAS¹³ was used for data reduction. During the data reduction, we checked the spectrum of each scan and dropped scans with bad baseline and/or hot spots in some channels, which comprised about 5% of the total data. We combined the data from all four days to obtain the final spectrum weighted by the rms of each individual scan. The rms of the final spectrum is 0.38 mK in T_mb after smoothing to the velocity resolution of 98.7 km s⁻¹, while the total effective observing time (on+off) is 408 minutes.

The NOEMA observations were carried out on 2017 December 1 and 3 (Project id: d17ae001, PI: Junzhi Wang) under good weather conditions with the same pointing center as that of IRAM 30 m observation in 2015. The frequency coverage of the upper sideband is from about 108.2 to 116.2 GHz, which covers the redshifted CN 1–0, CO 1–0, O₂, and HC₃N J = 13–12 lines. The original frequency resolution was 2 MHz, which corresponds to 5.26 km s⁻¹ at 114 GHz. On December 1, 3C84 was used as flux and bandpass calibrator, while 1300+580 and 1150+497 were used as phase calibrators, with seven antennae in the array and a total of 2.6 hr time on Mrk 231. There were two observing blocks on December 3 using 3C273 as flux and bandpass calibrator while the phase calibrators were again 1300+580 and 1150+497. One block was done with eight antennae in the array and a total of 0.8 hr on Mrk 231, while the other was done with nine antennae in the array and a total of 3.8 hr on Mrk 231.

After accounting for bandpass, flux, phase, and amplitude calibrations for the visibility data with "CLIC" in GILDAS, we carried out imaging and deconvolving with both "MAPPING" in GILDAS and also CASA. The final results from both software packages are consistent. The images and spectra presented in this paper were done with CASA after smoothing to 10 MHz frequency resolution corresponding to a velocity resolution of 26.3 km s⁻¹, and giving a slightly noncircular beam of 3.63'' × 3.32'' with PA = −65°.

An emission feature was detected at about the 12σ level, which we identify as O₂ N_J = 1₁–1₀ emission with similar velocity range to that of the red wing of outflowing gas observed in CO J = 1–0. The velocity range of this feature is about 150–650 km s⁻¹ (see Figure 1). The frequency resolution of the spectrum shown in Figure 1 is 37.5 MHz, which corresponds to 98.7 km s⁻¹ at 113.945 GHz for the O₂ N_J = 1₁–1₀ line with z = 0.04217 (Bryant & Scoville 1996; Carilli et al. 1998). This redshift corresponds to 12,642 km s⁻¹ as the optically defined velocity for this galaxy from NASA/IPAC Extragalactic Database (http://ned.ipac.caltech.edu).

Zoom In Zoom Out Reset image size

The velocity-integrated flux of O₂ N_J = 1₁–1₀ emission is 0.88 ± 0.07 K km s⁻¹, which corresponds to 4.4 ± 0.35 Jy km s⁻¹ with the conversion factor of 5.0 from K to Jy. The error of 0.07 K km s⁻¹ is from single component Gaussian fitting for the spectrum, while there will be about 20% more uncertainties due to absolute flux calibration which should be considered if comparing with different observations. The flux of detected O₂ emission in Mrk 231 is comparable to that of the red wing of CO 1–0 previously detected (Cicone et al. 2012), which is 2.32 ± 0.17 Jy km s⁻¹. Our CO J = 1–0 data obtained simultaneously with the O₂ in our IRAM 30 m observations provided similar results for the red wing of CO 1–0, with similar flux level and velocity range (see Figure 1). In addition to O₂ N_J = 1₁–1₀ line, we also detect HC₃N J = 13–12 line, at a significance of 8σ. This line is at ∼1200 km s⁻¹ in the figures as the O₂ N_J = 1₁–1₀ frequency corrected central velocity of the galaxy was used as reference frequency. The data are shown at higher velocity resolution and wider velocity coverage in Figure 2.

Zoom In Zoom Out Reset image size

Since the on-source times on August 20 and 21 are about half of those on August 22 and August 23, and the system temperatures on August 20 and 21 are higher than those on August 22 and 23, most of which are lower than 300 K, the main contributions to the final spectrum are from the data on August 22 and 23. We present the data of August 22 and 23 as the two spectra in Figure 3. Even though the noise levels of the two spectra are much higher than the 4 day combined spectrum, the two emission features, especially the one in the range ∼100–800 km s⁻¹, can be seen in the individual spectra shown in Figure 3. On the other hand, due to very little on-source time in August 20 and 21 without any emission signature, the spectra are not shown here.

Zoom In Zoom Out Reset image size

The single-dish detections of O₂ N_J = 1₁–1₀ and HC₃N J = 13–12 emission in Mrk 231 were confirmed with observations using NOEMA. O₂ N_J = 1₁–1₀ emission detected with NOEMA extends to ∼15'' (1'' = 0.87 kpc for Mrk 231) from the center of the galaxy (see Figure 4). After the inclination correction with the average value of cos(θ) for random sample, the physical projection will be $15\times 0.87\times \tfrac{\pi }{2}=20.5$ kpc. However, since the real inclination is unknown, the distances to center in the discussion part are without such correction. The total flux of the identified emission features observed with NOEMA (see Table 1) is 1.9 ± 0.2 Jy km s⁻¹, which includes about 43% of the emission detected with the 30 m telescope, and covers a consistent velocity range, from about 180 to 570 km s⁻¹. Since the signal is not a point source located at the phase center, which leads to the difficulty of obtaining the emission features in UV plane data, the UV spectrum was not presented.

Zoom In Zoom Out Reset image size

Continuum subtraction was implemented in the cleaned image data cube based on nearby line free channels, from about −2500 to −1000 km s⁻¹ and from about 1600 to 3000 km s⁻¹. HC₃N 13–12 was seen as compact emission in the nuclear region, coincident with the AGN/extreme starburst and consistent with the spectrum taken by the IRAM 30 m telescope, while O₂ emission was observed to be extended in outer disk regions ∼10 kpc from the nucleus. The emission was detected in several different regions with different velocities (see Figures 4 and 5). The flux of the HC₃N 13–12 line detected with NOEMA is 1.55 ± 0.07 Jy km s⁻¹, which agrees with the total emission detected with the IRAM 30 m single-dish telescope within the uncertainties of the observations at the two facilities.

Zoom In Zoom Out Reset image size

We present the velocity-integrated maps for different velocity ranges in Figure 5, smoothed to about 150 km s⁻¹ velocity resolution. The noise level (1σ) is about 0.045 Jy km s⁻¹ for each map, while the contours are from 2σ with 1σ steps. Emission features can be seen in several velocity components above 3σ, especially from 180 to 570 km s⁻¹, which is consistent with the emission feature detected with the IRAM 30 m telescope. The O₂ moment zero map (Figure 4, in both gray scale and contours) includes emission features in the velocity range 180–570 km s⁻¹. The corresponding angular regions are marked in Figure 5. With information from NOEMA observations, even the velocity range of O₂ N_J = 1₁–1₀ emission is similar to that of CO red wing, which traces molecular outflows within 3'' (Feruglio et al. 2010; Cicone et al. 2012), O₂ is from different regions of outflowing gas traced by CO red wing.

The image rejection of the EMIR on the IRAM 30 m telescope is better than 13 dB, which corresponds to ≥20:1. The lower sideband image frequency corresponding to the feature of interest in the upper sideband is different on different days, since different LO frequencies were used. The result is that any signal from the lower sideband will move around in frequency rather than adding together when the data from different days are added. The image frequencies are in the range 96–97.5 GHz during the four days, and are not near the strong dense gas tracer CS J = 2–1 at the rest frequency of 97.980953 GHz. Even the strongest dense gas tracer in the 3 mm band, HCN 1–0, produced an antenna temperature of only 6.5 mK in Mrk 231 with the IRAM 30 m telescope (Jiang et al. 2011). The lower sideband image frequencies with NOEMA corresponding to the upper sideband emission feature identified as O₂ N_J = 1₁–1₀ are between 92 and 94 GHz (due to the different IF frequency), and there were no emission features detected in the lower sideband, while the detected CS J = 2–1 line with NOEMA observation is about half the intensity of HCN J = 1–0. The expected confusion from the lines in the image sideband should thus be much less than 0.2 mK for the final spectrum, which does not affect the line identification.

The rest frequency range of the emission feature from −100 to 900 km s⁻¹ in Figure 1 is from about 118.4 to 118.8 GHz. From splatalogue (http://www.cv.nrao.edu/php/splat/), which combines the JPL molecular catalog (Pickett et al. 1998) and the CDMS catalog (Müller et al. 2001), there are 350 lines in this frequency range, including the O₂ N_J = 1₁–1₀ line at 118.75034 GHz, which is the only one that has been detected in space (Larsson et al. 2007). With the criterion of lower level energy less than 100 K, the number of lines is decreased to 83, which are mainly from complex molecules, such as CH₃OCH₃ and CCCN, or isotopologues, such as ³³SO₂. The abundances of such complex molecules and isotopologues are not expected to be high enough to produce emission stronger than HC₃N J = 13–12 in Mrk 231. Otherwise, such features should be easily detected in Galactic molecular clouds, while HC₃N J = 13–12 was the strongest (and possibly the only) emission line detected in the frequency range from 117.8 GHz to 118.75 GHz toward the Sgr A complex region in the Galactic center with Odin (Sandqvist et al. 2008).

If we restrict the emission feature to fall within ±100 km s⁻¹ of the central velocity of the galaxy, with the rest frequency range from 118.54 to 118.62 GHz, there are 57 lines from splatalogue, none of which has been reported as detected. All of them, except for SO₂ at 118.57743 GHz and the vibrationally excited HC₃N J = 13–12 1e line at 118.56145 GHz, are extremely complex molecules with more than five atoms. The lower level energy of SO₂ and HC₃N transitions are from about 250 to 1600 K, which makes them unlikely candidates for the observed emission. In addition, the vibrationally excited line would almost certainly be significantly weaker than the corresponding ground vibrational state line, as found even in regions of high infrared flux that can radiatively excite the vibrational levels, e.g., Goldsmith et al. (1982).

If the origin of the emission were a low-lying transition of a complex molecule, it should have been detected in sources in the Milky Way. However, no such emission feature was reported by Sandqvist et al. (2008). Thus, the most probable identification of the emission feature from 100 to 800 km s⁻¹ is O₂ N_J = ${1}_{1}$–1₀ with the velocity of the red wing of the gas in Mrk 231.

Stellar components traced at optical wavelengths can be seen (Rupke & Veilleux 2011) in the region with high velocity redshifted O₂ emission, while studies of molecular gas traced by CO lines (Feruglio et al. 2010; Cicone et al. 2012; Feruglio et al. 2015) were focused on the central 2–3 kpc region, even though some outflow features can be seen up to 5'' (∼4 kpc) north of the central region (Cicone et al. 2012). In the region with O₂ detection extending from 450 to 570 km s⁻¹, the flux of the simultaneously obtained CO J = 1–0 line was about four times that of the O₂ N_J = 1₁–1₀ line. However, the velocity range of CO 1–0 is primarily associated with the central region of Mrk 231, which is quite distinct from the O₂ emission. Thus, even though the O₂ emission has a velocity range similar to the red wing of molecular outflow traced by CO 1–0, it is coming from a different portion of Mrk 231. Molecular line emission of CH⁺ J = 1–0 without associated CO emission has been found in several high redshift starburst galaxies (Falgarone et al. 2017), which indicates that CO emission is not necessarily detectable in the regions in which some special molecules, such as O₂ and CH⁺ can be detected.

The CO J = 1–0 fluxes in the regions with O₂ detections were comparable to, or even less than that of O₂ N_J = 1₁–1₀ line, if we consider that only the CO 1–0 emission with the same velocity range as the O₂ line is from the same gas as responsible for the O₂ emission. Using the same CO 1–0 flux to H₂ conversion factor in the O₂ emission region as that of the host galaxy and assuming that the O₂ line is optically thin with a 15 K excitation temperature (T_ex), which is the same as was assumed in NGC 6240 (Combes et al. 1991), the O₂–H₂ abundance ratio can be estimated to be higher than 10⁻⁴ (see Table 2), nearly two orders of magnitude higher than that of Galactic sources detected in O₂, and higher than the CO-to-H₂ abundance ratio. The dependence of the O₂ column density on its excitation temperature in LTE is N(O₂) ∝ ${T}_{\mathrm{ex}}^{0.67}$ (Liseau & O₂din Team 2005). Thus if the excitation temperature of O₂ were higher than 15 K, the O₂–H₂ abundance ratio would be even higher.

The nondetection of O₂ emission in the central 2 kpc of Mrk 231 with our NOEMA observations gave the best upper limit at the extragalactic O₂ N_J = 1₁–1₀ to CO 1–0 line ratio of 1 × 10⁻³ (1σ), much lower than the 1.2 × 10⁻² value for NGC 6240, which was used previously for estimating the O₂/H₂ abundance ratio (Combes et al. 1991). Using the same CO J = 1–0 flux to H₂ conversion factor as that in NGC 6240 and similar excitation conditions of the O₂ molecules (Combes et al. 1991), the O₂ to H₂ abundance ratio should be less than 8 × 10⁻⁸ in the host galaxy of Mrk 231, which is consistent with that found in the vast majority of dense molecular gas in the Milky Way (Goldsmith et al. 2000; Pagani et al. 2003). The only exceptions are Orion (Goldsmith et al. 2011) (having significantly higher abundance in 10⁻⁶ in a very small region) and ρ Oph (Larsson et al. 2007; Liseau et al. 2012) (with comparable abundance 5 × 10⁻⁸).

We compiled the sources with O₂ detections, the best upper limit of O₂–H₂ abundance ratio in NGC 6240 with emission (Combes et al. 1991) and in B0218+357 with absorption (Combes et al. 1997) in the literature, and show the results in Table 2. If we wish to detect O₂ N_J = 1₁–1₀ emission in extragalactic sources without O₂ abundance enhancement, observations of the redshifted O₂ N_J = 1₁–1₀ line with five times lower noise level than that of our observations with NOEMA toward Mrk 231 are needed. With the Atacama Large Millimeter/submillimeter array (ALMA), about 10 hours on-source time toward nearby gas-rich galaxies with bright CO emission can achieve such sensitivity.

The velocity range of our O₂ emission is mainly from redshifted gas, especially from ∼345 to 615 km s⁻¹, with projected distance of about 10 kpc north and south from the central region, while the low velocity components, from about −400 to +300 km s⁻¹, were seen as weaker emission and at smaller projected distances to the center (see Figure 4). On the other hand, high velocity blueshifted components were even weaker than the low velocity ones and were hardly detected in this observation. The asymmetry in spatial and velocity distribution of O₂ emission, including the absence of O₂ emission along the line of sight around the nuclear region, can be attributed to the geometry of molecular gas in Mrk 231 and the asymmetry of the molecular outflow from the central AGN, which had been found in CO observations up to about 3 kpc scale away from the center of Mrk 231 in projection (Feruglio et al. 2010; Cicone et al. 2012).

The absence of a blueshifted component of O₂ emission may be caused by asymmetric distribution of molecular gas in the regions about 10 kpc away from the center, because O₂ should be produced by the interaction between the outflowing gas from the AGN center and the local molecular gas there. With velocity of about 500 km s⁻¹, 20 Myr are required for gas to travel 10 kpc, which is the typical projected distance between the O₂ emission regions and the center of Mrk 231. Thus O₂ emission in Mrk 231 traces the outflow from the central AGN on a characteristic dynamic timescale of 20 Myr.

X-ray emission extended up to 25 kpc away from the center of Mrk 231 was found with half-mega second Chandra spectral imaging. A merger remnant was proposed as the main contributor to the large-scale X-ray halo of Mrk 231 (Veilleux et al. 2014). However, the central AGN is responsible for emission at smaller (≤6 kpc) distances. We suggest that large scale AGN molecular outflow is the most plausible mechanism to produce a shock that drives the enhancement of the O₂ abundance but a merger is another possibility for producing a shock (or multiple shocks) that enhanced the abundance of O₂.

Models of shock chemistry (Chen et al. 2014; Neufeld et al. 2014; Melnick & Kaufman 2015; Godard et al. 2019) suggest that essentially all of the gas-phase oxygen will be incorporated into H₂O and CO if the shock velocity exceeds 20 km s⁻¹. Thus, a single shock that provides the velocity offset of the redshifted O₂ emission of Mrk 231 could not be responsible for dramatically enhancing the abundance of molecular oxygen. Rather, a number of lower velocity shocks produced in the outflow, having different velocities relative to the general outflow, could collectively provide the observed large O₂ abundance. A high O/CO ratio up to 2 × 10⁴ can be produced in a model of a UV irradiated molecular shock (Godard et al. 2019). Such irradiated molecular clouds with shocks located far away from the center of galaxies can have physical properties similar to those of the O₂-emitting regions in Mrk 231. Chemical networks coupled with dynamic models including radiation field are necessary for understanding the high O₂/CO abundance ratio there.

For 15 K gas under LTE conditions, the 118.75034 GHz O₂ transition is the brightest among the 114 O₂ lines listed between 10 and 1000 GHz (Pickett et al. 1998). In total, only 12 transitions make meaningful contributions, which is consistent with what had been predicted with an LVG model (Liseau & O₂din Team 2005). The total O₂ flux amounts to 2.9 times the flux at 118 GHz. Similar procedures applied to 17 transitions of CO yielded a CO total flux 5.1 times that of CO 1–0. Given the comparable fluxes of the O₂ N_J = 1₁–1₀ transition and the CO 1–0 line in outer disk region of Mrk 231, the O₂ molecule is comparable to CO as an important coolant of such molecular gas.

Based on CO observations (Cicone et al. 2014; Fluetsch et al. 2019), massive molecular outflows, which have negative feedback on star formation in galaxies, are quite common in AGN. The O₂ N_J = 1₁–1₀ transition is thus a promising new tool for studying such outflows on dynamic timescales of a few tens of megayears. AGNs with molecular gas detected via strong CO emission will be good candidates for searching for O₂ emission in extragalactic sources.