<figure><img src="image_1.png" /></figure>

Deep Observations of CO and Free-Free Emission in Ultraluminous Infrared QSO IRAS F07599+6508

Qing-Hua Tan¹ Ⓓ, Yu Gao²,¹ Ⓓ, Emanuele Daddi³ Ⓓ, Xiao-Yang Xia⁴, Cai-Na Hao⁴ Ⓓ, Alain Omont⁵ Ⓓ, and Kotaro Kohno⁶,⁷ Ⓓ

¹Purple Mountain Observatory & Key Laboratory for Radio Astronomy, Chinese Academy of Sciences, 10 Yuanhua Road, Nanjing 210023, People's Republic of China; qhtan@pmo.ac.cn

²Department of Astronomy, Xiamen University, 422 Siming South Road, Xiamen 361005, People's Republic of China; yugao@xmu.edu.cn

³CEA, IRFU, DAp, AIM, Université Paris-Saclay, Université de Paris, CNRS, F-91191 Gif-sur-Yvette, France

⁴Tianjin Astrophysics Center, Tianjin Normal University, Tianjin 300387, People's Republic of China

⁵Sorbonne Université, UPMC Université Paris 6 and CNRS, UMR 7095, Institut d'Astrophysique de Paris, 98bis boulevard Arago, F-75014 Paris, France

⁶Institute of Astronomy, Graduate School of Science, The University of Tokyo, Osawa, Mitaka, Tokyo 181-0015, Japan

⁷Research Center for the Early Universe, Graduate School of Science, The University of Tokyo, 7-3-1 Hongo, Bunkyo-ku, Tokyo 113-0033, Japan

*Received 2020 December 9; revised 2021 March 29; accepted 2021 March 31; published 2021 May 27*

**Abstract**

Infrared quasi-stellar objects (IR QSOs) are a rare subpopulation selected from ultraluminous infrared galaxies (ULIRGs) and have been regarded as promising candidates of ULIRG-to-optical QSO transition objects. Here we present NOEMA observations of the CO (1–0) line and 3 mm continuum emission in IR QSO IRAS F07599+6508 at $z=0.1486$, which has many properties in common with Mrk 231. The CO emission is found to be resolved with a major axis of ~6.1 kpc that is larger than the size of ~4.0 kpc derived for 3 mm continuum. We identify two faint CO features located at a projected distance of ~11.4 and 19.1 kpc from the galaxy nucleus, both of which are found to have counterparts in the optical and radio bands and may have a merger origin. A systematic velocity gradient is found in the CO main component, suggesting that the bulk of molecular gas is likely rotationally supported. Based on the radio-to-millimeter spectral energy distribution and IR data, we estimate that about 30% of the flux at 3 mm arises from free-free emission and infer a free-free-derived star formation rate of $77 M_{\odot} \text{ yr}^{-1}$, close to the IR estimate corrected for the AGN contribution. We find a high-velocity CO emission feature at the velocity range of about –1300 to –2000 km s$^{-1}$. Additional deep CO observations are needed to confirm the presence of a possible very high velocity CO extension of the OH outflow in this IR QSO.

*Unified Astronomy Thesaurus concepts: Starburst galaxies (1570); Galaxy mergers (608); AGN host galaxies (2017); Interstellar molecules (849); Galaxy evolution (594)*

**1. Introduction**

Mergers of gas-rich galaxies may trigger both intense star formation and nuclear activity, leading to a violent transition from an obscured accretion and starburst stage to an optically bright quasi-stellar object (QSO) and eventually an elliptical galaxy (e.g., Sanders et al. 1988a, 1988b; Barnes & Hernquist 1992; Di Matteo et al. 2005; Hopkins et al. 2006). In this merger-driven evolutionary scenario, the active galactic nucleus (AGN) was thought to be energetically dominant but obscured in the ultraluminous infrared galaxy (ULIRG; $L_{\text{IR}} \ge 10^{12} L_{\odot}$) phase, where the most luminous sources that reach a QSO-like luminosity represent the rare ULIRG-to-QSO transition objects. These systems are very likely experiencing the so-called feedback or blowout phase that expels the gas of the host galaxy in the form of outflowing winds, which have been invoked as one of the main drivers for regulating the growth of the stellar mass and black hole (BH) accretion and producing the tight BH-bulge mass relation (see Silk & Rees 1998; Hopkins et al. 2008; Fabian 2012; Kormendy & Ho 2013; King & Pounds 2015, and references therein), and therefore represent a key phase in the lifetime of a quasar.

The large gas masses discovered in local ULIRGs via CO observations have established the presence of abundant fuel for the intense star formation activity (e.g., Solomon et al. 1997). Along with the clear evidence that many (U)LIRGs are strong dynamical interactions and mergers (Sanders & Mirabel 1996), all the properties make these systems unique local examples of dust-enshrouded galaxies at high z. Xia et al. (2012) have

investigated for the first time the CO survey in 19 ultraluminous infrared type 1 QSOs, selected from the local ULIRGs and referred to as IR QSOs (see Zheng et al. 2002). These IR-selected type 1 QSOs (IR QSOs) were originally compiled from the ULIRGs in the QDOT redshift survey (Lawrence et al. 1999), the 1 Jy ULIRG survey (Kim & Sanders 1998), and the cross-correlation of the IRAS Point Source Catalog with the ROSAT All-sky Survey Catalog. Studies of luminosity function for large samples of IRAS galaxies reveal that the space densities of ULIRGs and QSOs are comparable in the local universe (Soifer et al. 1987), while the fraction of IR QSOs is less than 10% of ULIRGs (Lawrence et al. 1999), indicating that the IR QSO phase may last only a few times $10^7$ yr if the space density of objects is simply related to the timescale of different phases (Hao et al. 2005).

The CO observations of IR QSOs show that the molecular gas masses are a few times $10^9-10^{10} M_{\odot}$ (Xia et al. 2012), similar to that found in local ULIRGs but significantly higher than that of optically selected PG QSOs (e.g., Scoville et al. 2003; Evans et al. 2001, 2006; Shangguan et al. 2020). This supports the scenario that IR QSOs are rich in fuel to feed both the star formation and the accretion of AGNs. A similar conclusion has also been reached recently for type 2 QSOs, of which a sample of QSO hosts are found to be gas-rich based on the CO observations (Jarvis et al. 2020). Together with the far-IR (FIR) excess observed in IR QSOs (Hao et al. 2005) and the multiwavelength properties (i.e., large blueshifts in optical emission lines, warm IR color, intermediate slope of mid-IR spectral energy distribution (SED) between ULIRGs and PG
QSOs, and high accretion rates; Zheng et al. 2002; Cao et al. 2008), all of these properties imply that IR QSOs are likely the objects caught in the short-lived transition phase between ULIRG and QSO stages (Canalizo & Stockton 2001).

Unlike local (U)LIRGs, a large number of which have been spatially resolved in molecular gas emission (e.g., Downes & Solomon 1998; Wilson et al. 2008; Ueda et al. 2014; Xu et al. 2014, 2015; Imanishi et al. 2019), only the nearest IR QSO, Mrk 231 at $z=0.04$, has been well studied in multiwavelength and multiscale (e.g., Bryant & Scoville 1996; Taylor et al. 1999; Fischer et al. 2010; Rupke & Veilleux 2011; Feruglio et al. 2015; Cicone et al. 2020). Until recently, the CO emission in a handful of IR QSOs has been spatially resolved on subkiloparsec scales with ALMA (Tan et al. 2019). It was found that the host galaxies of IR QSOs resolved with ALMA have a diverse variety of CO morphology and kinematics (i.e., rotating disk, disturbed structure, and multiple interacting sources), indicative of more complicated evolutionary stages from merging (U)LIRGs to QSOs. This is somewhat unexpected since the objects in IR QSO phase are typically characterized by the final coalescence of the galaxies, according to the models of the merger-driven evolutionary scenario (Barnes & Hernquist 1992; Hopkins et al. 2008). However, such a classification suffers from small number statistics, and hence spatially resolved CO observations in more IR-luminous QSOs are required for a meaningful statistical study.

In this paper, we present CO (1–0) imaging of IRAS F07599+6508 obtained with the IRAM NOrthern Extended Millimeter Array (NOEMA) interferometer. IRAS F07599+6508 is an IR QSO at $z=0.149$ with a total IR luminosity of $3.5 \times 10^{12} L_{\odot}$ and an AGN bolometric luminosity of $L_{\text{AGN,bol}} = 3.9 \times 10^{12} L_{\odot}$ (Hao et al. 2005). It is an active galaxy of particular interest because it is highly similar to Mrk 231, the closest and most prominent template of local IR QSOs, but at a higher redshift and with a lack of detailed studies on the cold molecular gas of the host galaxy, especially the spatial distribution and kinematics. Similar to Mrk 231, IRAS F07599+6508 is a low-ionization broad absorption line (BAL) QSO with strong Fe II emission (e.g., Kwan et al. 1995; Véron-Cetty et al. 2006), which is a very rare subclass of BAL QSOs that may represent young, heavily enshrouded AGNs, where the cocoon of gas and dust has a large covering factor (Boroson & Meyers 1992; Zheng et al. 2002).

The IRAM 30 m CO observations show massive amounts of molecular gas residing in IRAS F07599+6508 and a broad FWHM line width of $380 \pm 14 \text{ km s}^{-1}$ and $490 \pm 40 \text{ km s}^{-1}$ for the CO (1–0) and CO (2–1) lines, respectively (Xia et al. 2012). Based on six independent methods, Veilleux et al. (2009) estimate an average fraction of AGN contribution to the bolometric luminosity of $\alpha_{\text{AGN}} = 0.75$, consistent with other IR QSOs (e.g., ~70% for Mrk 231), where the $\alpha_{\text{AGN}}$ is found to be typically higher than that (~35%–40%) of local cooler ULIRGs. The weakness of both soft ($L_{0.1-2.4 \text{ keV}} = 6.9 \times 10^{42} \text{ erg s}^{-1}$) and hard X-ray emission ($L_{2-10 \text{ keV}} = 1.3 \times 10^{42} \text{ erg s}^{-1}$) seen in IRAS F07599+6508 (e.g., Xia et al. 2001; Brightman & Nandra 2011; Luo et al. 2014; La Caria et al. 2019), i.e., the X-ray emission is not at the level expected from its optical/UV emission compared to typical QSOs, implies that its central AGN is very likely highly obscured.

The existence of multiphase outflow winds in IRAS F07599+6508 has been demonstrated by the detection of large blueshifts

(~2000 km s⁻¹) of optical emission lines, Na I D interstellar absorption lines, and an OH 119 µm P Cygni/blueshifted absorption profile, closely resembling Mrk 231 (Boroson & Meyers 1992; Zheng et al. 2002; Spoon et al. 2013; Veilleux et al. 2013; Rupke et al. 2017). The sensitive CO (1–0) mapping observations will help constrain the distribution and kinematics of the molecular gas reservoir in IRAS F07599+6508 and offer the possibility of detecting the molecular outflow based on the presence of broad wings in the CO line emission. A tremendous amount of observing time was devoted to Mrk 231 on NOEMA, revealing the outstanding molecular outflows, including in dense gas tracers (Aalto et al. 2012), and likely the tentative detection of large-scale molecular oxygen emission (Wang et al. 2020).

Similar to Mrk 231, IRAS F07599+6508 has typical features expected for a system transiting from the heavily obscured starburst phase to the unobscured QSO phase and could also be regarded as a representative of IR QSOs that have low number density in the local universe but are expected to be common at high redshift (Veilleux et al. 1999; Le Floc'h et al. 2005), and thus it provides a nearby laboratory to probe the ULIRG-to-optical QSO transition phase that may hold the key to understanding the galaxy evolution in a merger-driven scenario.

The paper is organized as follows. In Section 2, we describe the NOEMA observations and data reduction process. Section 3 presents the results from the analysis of 3 mm continuum and CO (1–0) line data, including morphology, spatially decomposed CO components, kinematics, and an analysis of archival radio continuum data from VLA. In Section 4, we discuss the nature of 3 mm continuum emission, the origin of the off-center CO components, and make a comparison of the molecular gas properties between IRAS F07599+6508 and Mrk 231. We summarize our main findings in Section 5. Throughout this work, we assume a standard $\Lambda$CDM cosmology with $H_0 = 70 \text{ km s}^{-1} \text{ Mpc}^{-1}$, $\Omega^{m} = 0.3$, and $\Omega^{\Lambda} = 0.7$.

## 2. Observations and Data Reduction

### 2.1. NOEMA CO (1–0) Observations

We have observed the CO (1–0) line and the continuum emission in IRAS F07599+6508 with the IRAM NOEMA interferometer between 2018 May and July. The observations were obtained with a total on-source time of 12.7 hr in the compact array D configuration (with eight antennas covering baselines of 15.2–175.9 m). We used the PolyFix correlator in 3 mm Band 1, which provides $2 \times 8 \text{ GHz}$ of instantaneous dual-polarization bandwidth with a spectral resolution of 2 MHz (corresponding to $6.0 \text{ km s}^{-1}$ at the redshifted CO (1–0) frequency of 100.353 GHz). The tuning frequency was set to 98.2 GHz, with frequency ranging from 94.2 to 102.2 GHz in the upper sideband (USB) and from 78.8 to 86.8 GHz in the lower sideband (LSB). In this paper, we analyze the CO (1–0) line and the continuum data in the USB; the rest of the detected lines (i.e., CN (1–0) lines at 113.191 and 113.491 GHz, CO (1–0) isotopologue lines, and the 109.174 GHz $\text{HC}_3\text{N}$ (12–11) line) covered in the USB and the LSB line data will be presented in a future paper.

The data calibration was performed using the GILDAS<fnref n="8" /> package CLIC with help from the staff in Grenoble. The typical uncertainty of flux calibration in the 3 mm band is about 10%. We binned the u-v tables by $20 \text{ km s}^{-1}$ and redefined the rest

<footnote n="8">⁸ http://iram.fr/IRAMFR/GILDAS; for more information, see Guilloteau & Lucas (2000).</footnote>
<figure><img src="image_2.png" /><figcaption>**Figure 1.** NOEMA cleaned map of 3 mm continuum emission (left), the residual map of 3 mm emission after subtracting a circular Gaussian model (middle), and the corresponding *u*-v plot showing the visibility amplitude as a function of the *u*-v distance (right) for IRAS F07599+6508. For the 3 mm continuum map, contours start at ±2σ (1σ rms noise is 10 µJy beam⁻¹) and increase by a factor of 1.5, with positive (negative) contours shown as solid (dotted) lines, while for the residual map, contours start at ±2σ and increase in steps of 1σ. The synthesized beam of 6.25 × 4.22 with position angle 106.2° is shown in the lower left corner. The plus sign indicates the centroid of the 3 mm continuum emission derived from a circular Gaussian fit to the *u*-v data. In the right panel, the black circles with error bars shown in the *u*-v plot are the data binned in *u*-v radius steps of 8 kλ. Error bars show the statistical photon noise on the average amplitude in each bin. The red curve indicates the best-fitting model with a circular Gaussian function to the *u*-v data.</figcaption></figure>

frequency to the sky frequency of CO (1–0) transition. We then used the software GILDAS/MAPPING for image cleaning and analysis in the *u*-*v* plane. We also exported the data to a CASA <fnref n="9" />-compatible format for a comparison analysis. The continuum *u*-*v* table was created from the calibrated visibilities by averaging the line-free channels using the MAPPING task *uv.continuum* with line emission channels filtered (*uv_filter* task). The continuum emission was subtracted from the line *u*-*v* data by fitting a first-order polynomial baseline for each visibility (*uv_baseline* task). The USB and LSB continuum *u*-*v* data were combined using the *uv_merge* task to form a single output table centered at 90.5 GHz. The synthesized beam size using natural weighting is 5.67 × 3.81, and the NOEMA primary beam at the sky frequency of CO (1–0) is 50.2 × 50.2. We obtain a sensitivity of ~0.44 mJy beam⁻¹ in the spectral velocity resolution of 20 km s⁻¹ and ~10 µJy beam⁻¹ for the continuum image.

## 2.2. Archival VLA Radio Data

In this work we also make use of ancillary radio continuum data observed with the Very Large Array (VLA) for a comprehensive analysis. IRAS F07599+6508 has been observed at 1.49, 4.86, 8.44, 14.94, and 22.46 GHz in CD configuration (project ID: AB0783). Higher-resolution observations in AB configuration (project ID: AH0406) were available at 1.49, 4.86, and 8.44 GHz, as well as in A configuration at 14.94 GHz (project ID: AN0104). We obtained the raw data sets from the VLA archive. <fnref n="10" /> The data flagging and reduction were performed by using the AIPS package (van Moorsel et al. 1996). The full image at each frequency and array was processed individually with the AIPS CLEAN task *imagr* and then transferred to FITS format for further analysis with the CASA package. The synthesized beam using Briggs weighting with a robustness parameter of 0.5 ranges from 65.2 × 28.3 at 1.49 GHz to 1.67 × 1.30 at 43.34 GHz for

observations in CD configuration and from 5.33 × 2.08 at 1.49 GHz to 1.8 × 1.3 at 14.94 GHz for higher-resolution observations in AB and A configurations.

# 3. Results and Analysis

## 3.1. 3 mm Continuum Emission

The continuum emission at 90.5 GHz (rest-frame 103.9 GHz; ~3 mm) is clearly detected toward the nucleus of IRAS F07599+6508. Figure 1 shows the 3 mm continuum map and the corresponding *u*-*v* plot that shows the amplitude of the visibilities as a function of the *u*-*v* distance. It can be seen from the *u*-*v* data that the 3 mm continuum source is resolved in our observations, as an unresolved source would give a flat distribution of visibilities.

Assuming that the surface brightness distribution of the continuum follows a symmetric Gaussian profile, we performed a fit to the visibilities in GILDAS using task *uv_fit* and found that the continuum source is best fitted by a 2D circular Gaussian model with a deconvolved FWHM size of 1.53 ± 0.19 (~4.0 ± 0.5 kpc) and a flux density of 560 ± 20 µJy (see Table 1). A 2D elliptical Gaussian model fit to the continuum visibilities returns a deconvolved FWHM size of (1.54 ± 0.22) × (1.50 ± 0.39) and a flux density of 560 ± 20 µJy, which are in good agreement with those derived from the fit with a circular Gaussian model, but with slightly larger uncertainties on the size measurement. In addition, the position angle (0° ± 293°) derived from the elliptical Gaussian fit is rather uncertain owing to the large uncertainty.

In the middle panel of Figure 1, we show the residual after removing the fit with a circular Gaussian model. No significant peaks or negative regions with signal-to-noise ratio (S/N) ≥ 3 are found in the residual map. The red curve in the *u*-*v* plot (see the right panel of Figure 1) shows the result of the circular Gaussian model fit. Moreover, as can be seen from the right panel of Figure 1, most data points are consistent with the model fit (see the red curve) with a circular Gaussian function within 1σ. All of these suggest that a single circular Gaussian model may provide a good description to the 3 mm continuum emission based on the current data.

<footnote n="9">⁹ http://casa.nrao.edu; see McMullin et al. (2007).</footnote>

<footnote n="10">¹⁰ https://science.nrao.edu/facilities/vla/archive/index</footnote>
<table><caption>**Table 1**
Measured Properties of CO (1–0) Line and 3 mm Continuum Emission of IRAS F07599+6508</caption><thead><tr><th>Emission<br>(1)</th><th>Component<br>(2)</th><th>R.A.<br>(J2000)<br>(3)</th><th>Decl.<br>(J2000)<br>(4)</th><th>Size FWHM<br>(arcsec)<br>(5)</th><th>FWHM CO (1-0)<br>(km s<sup>-1</sup>)<br>(6)</th><th>v<sub>offset</sub><br>(km s<sup>-1</sup>)<br>(7)</th><th>Integrated Flux<br>(8)</th></tr></thead><tbody><tr><td rowspan="3">CO (1-0)</td><td>F07599C</td><td>08:04:30.468 ± 0.004</td><td>64:59:52.72 ± 0.02</td><td>(2.35 ± 0.11)×(1.66 ± 0.09)</td><td>430 ± 10</td><td>0</td><td>11.24 ± 0.09<br>Jy km s<sup>-1</sup></td></tr><tr><td>F07599S</td><td>08:04:30.41 ± 0.03</td><td>64:59:48.36 ± 0.12</td><td>...</td><td>100 ± 10 (narrow)</td><td>100 ± 10</td><td>1.12 ± 0.07<br>Jy km s<sup>-1</sup></td></tr><tr><td>F07599NE</td><td>08:04:31.56 ± 0.07</td><td>64:59:55.2 ± 0.3</td><td>...</td><td>500 ± 80 (broad)<br>390 ± 70</td><td>40 ± 30</td><td>0.45 ± 0.07<br>Jy km s<sup>-1</sup></td></tr><tr><td>3 mm</td><td>total</td><td>08:04:30.432 ± 0.010</td><td>64:59:52.71 ± 0.04</td><td>1.53 ± 0.19</td><td></td><td></td><td>560 ± 20 μJy</td></tr></tbody></table>

**Note.** Column (1): CO (1–0) line and 3 mm continuum emission. Column (2): components of CO (1–0) emission from source spatial decomposition in the *u-v* plane. Column (3) and (4): source position derived from model fit to the *u-v* data. Column (5): FWHM source size derived from model fitting to the visibilities using the GILDAS task *uv_fit*. Ellipses in the column indicate an unresolved component. The physical scale of 1'' corresponds to ~2.6 kpc. Column (6): FWHM CO line width derived from a Gaussian profile fit to the spectrum (except for F07599S, where the CO spectrum is fitted simultaneously with two Gaussians consisting of one narrow and one broad component centered at the same velocity). Column (7): CO (1–0) line centroid velocity offset from the systematic velocity. The centroid velocity is derived from a Gaussian fit to the spectrum. Column (8): CO (1–0) velocity-integrated intensity and 3 mm continuum emission flux. The CO integrated fluxes are measured from the spectra shown in Figure 4 by integrating the intensity over the line-emitting region in each channel. The continuum flux is measured from a circular Gaussian fit to the *u-v* data. Errors in fluxes quoted in this table are statistical and do not include the absolute flux calibration uncertainty.

As an independent check on the results from our analysis of *u-v* data, we also measured the size and flux of the continuum emission by fitting the data in the image plane using the *imfit* task in CASA. A 2D Gaussian fit to the 3 mm continuum map of Figure 1 finds a deconvolved FWHM size of (1''84 ± 0''45) × (1''40 ± 0''49) and an integrated flux of 550 ± 20 μJy, consistent with the estimates from the *uv_fit* within 1σ uncertainties.

## 3.2. CO (1–0) Line Emission

### 3.2.1. CO Morphology and Spatial Decomposition

Figure 2 shows the continuum-subtracted CO (1–0) spectrum and the integrated intensity map of IRAS F07599+6508. The CO spectrum is extracted by fitting a point-source model to the *u-v* table using the *uv_fit* task. By fitting a Gaussian profile to the spectrum, we measure a redshift of $z_{CO} = 0.1486$, being used to set the systematic velocity. The velocity width (FWHM) derived from the Gaussian fit to the spectrum is 430 ± 10 km s⁻¹. We note that the CO (1–0) line profile is asymmetrical and a single Gaussian function does not provide a good fit to the line peak and wings. A double-peak profile is detected at the line center, albeit not very prominent, which is often taken to indicate orbital motion in a merger or a rotating disk (e.g., Downes & Solomon 1998; Greve et al. 2005; Daddi et al. 2010). In addition, the line profile also shows a slightly weaker peak at a velocity of ~100–150 km s⁻¹, indicating that the motion of molecular gas in the host galaxy is complex.

The irregular morphology of the CO emission with extended structure is shown by the map in the right panel of Figure 2, obtained by collapsing the channels in the [-630, +450] km s⁻¹ velocity range where the flux is positive in each channel when the CO spectrum is rebinned to a channel width of ~60 km s⁻¹. The CO emission appears to be resolved with an extension to the northeast and to the south of the nucleus, respectively. To characterize the distribution of CO emission, we fitted an elliptical Gaussian model to the *u-v* data with channels averaged in the [-630, +450] km s⁻¹ velocity range using the task

*uv_fit*. The best fit gives a source size FWHM of (2''33 ± 0''11) × (1''91 ± 0''07). The peak of CO emission is found to be cospacial with the 3 mm continuum emission.

In the residual map after subtraction of an elliptical Gaussian component, we detect two off-center CO subcomponents at above 5σ, which we denote as northeast (NE) and south (S), and we refer to the main component fitted with an elliptical Gaussian in the galaxy center as F07599C (see Figure 3). To take into account the off-center components, we remake the fitting in the *u-v* plane using combinations of functions by allowing all parameters to be free in the modeling, and the details of the model fitting are described in Appendix A. The best fit shows that the CO emission can be spatially decomposed into three components, an elliptical Gaussian source (F07599C) and two point sources (F07599S and F07599NE).

For the main component F07599C, the model fit returns a source size FWHM of (2''35 ± 0''11) × (1''66 ± 0''09), corresponding to 6.1 × 4.3 kpc, which is about two times the median source size (~3.2 kpc) measured for the local IR QSOs with well-resolved CO morphology (Tan et al. 2019). The three components identified in the CO emission and their parameters are listed in Table 1.

### 3.2.2. Properties of the Decomposed CO Components

Figure 4 shows the CO (1–0) spectra of the three components and the corresponding velocity-integrated CO maps. The CO spectra of the three components are extracted simultaneously from the line cube in the *u-v* plane with fitting models identical to the best-fit models summarized in Table 1 (see Section 3.2.1). The spectrum of F07599C is obtained with the *u-v* modeling fit using an elliptical Gaussian function by fixing the parameters of centroid and major/minor axes, while for the off-center components, F07599S and F07599NE, which are offset from the nucleus of the main component with a projected distance of 4''4 (~11.4 kpc) and 7''4 (~19.1 kpc), respectively, the spectra are extracted by a point-source model with centroids fixed.
<figure><img src="image_3.png" /><figcaption>**Figure 2.** NOEMA continuum-subtracted spectrum (left) and cleaned map (right) of CO (1–0) emission of IRAS F07599+6508. The CO (1–0) spectrum is extracted from the central beam, and the spectrum is binned to a channel width of ~20 km s⁻¹. The red dotted line represents a Gaussian fit to the spectrum. The zero velocity corresponds to the redshift of zCO = 0.1486 derived from a Gaussian fit to the spectrum. The map of CO (1–0) emission is integrated in the [-630, +450] km s⁻¹ velocity range. Contours start at ±2σ (1σ rms noise is 0.092 Jy beam⁻¹ km s⁻¹) and increase by a factor of 1.5, with positive (negative) contours shown as solid (dotted) lines. The synthesized beam (5''61 × 3''77) is shown in the lower left corner of the map. The plus sign indicates the peak of 3 mm continuum emission.</figcaption></figure>

<figure><img src="image_4.png" /><figcaption>**Figure 3.** Residual maps of CO (1–0) emission after subtracting a single elliptical Gaussian model (left) and a combination of an elliptical Gaussian model and two point-source models (right). The black plus sign marks the peak of 3 mm continuum emission, while the two red crosses indicate the peak positions of the two point-like CO features. Contours start at ±2σ and increase in steps of 1σ for both maps. The two point-like subcomponents are labeled as F07599S and F07599NE, respectively.</figcaption></figure>

A single Gaussian fit to the line profile of F07599C and F07599NE shown in Figure 4 yields a peak flux density of 26.0 ± 0.3 mJy and 1.0 ± 0.2 mJy and an FWHM of 430 ± 10 km s⁻¹ and 390 ± 70 km s⁻¹, respectively. F07599NE has a velocity offset of 40 ± 30 km s⁻¹, consistent within the uncertainties with the systematic velocity. For F07599S, the spectrum is best fitted by two nested Gaussians (with a shared center), e.g., a narrow Gaussian to fit the central core of the line, and a broad Gaussian to fit the broad emission. The best fit gives a peak flux density of 3.9 ± 0.4 mJy and 1.3 ± 0.3 mJy and an FWHM of 100 ± 10 km s⁻¹ and 500 ± 80 km s⁻¹ with velocity centroid redshifted at a velocity of about 100 ± 10 km s⁻¹ for the narrow and broad emission, respectively. The results of Gaussian fits derived for the three components are reported in Table 1.

Similar to the CO spectrum extracted from the central beam, the CO line profile of F07599C shows characteristic double peak at the line center and a weaker peak at a velocity of ~100–150 km s⁻¹, which is found to be coincident with the CO line centroid velocity of F07599S (see Figure 4, top panel). The velocity coincidence is likely to suggest that a significant portion of molecular gas is detected in the overlap region between F07599C and F07599S, where the perturbation of molecular gas makes an imprint on the CO line profile. Higher spatial resolution observations are needed to resolve the emission and clarify the nature of the source.

The CO morphology of each component is shown by the map in the bottom right panel of Figure 4, obtained by averaging the velocity range of [-330, +330] km s⁻¹, [-270, +330] km s⁻¹,
<figure><img src="image_5.png" /><figcaption>**Figure 4.** Top: CO (1–0) spectra of the three components extracted from GILDAS task uv_fit using an elliptical Gaussian and two point-source functions for F07599C (orange), F07599S (purple), and F07599NE (green), respectively. The CO spectrum of F07599C is divided by 5 for comparison purposes. All spectra are smoothed to a velocity resolution of ~20 km s⁻¹. Bottom left: zoomed-in views of CO (1–0) spectra for F07599S and F07599NE; the spectra are rebinned in channels of ~60 km s⁻¹. The black lines show the Gaussian fits to the spectra. For F07599S, we modeled the line profile with two Gaussians (one narrow and one broad), while for F07599NE the spectrum is fitted by using a single Gaussian. Bottom right: corresponding continuum-subtracted CO (1–0) maps (colored contours) of the three components overlaid on the CO (1–0) map of the whole IRAS F07599+6508 (gray scale; same as the map in the right panel of Figure 2). The colored contours show the CO (1–0) emission integrated in the [−330, +330] km s⁻¹ (orange), [−270, +330] km s⁻¹ (purple), and [−330, +270] km s⁻¹ (green) velocity range, respectively, while the black contour indicates the global CO (1–0) emission detected at a 3σ significance in the intensity map integrated in the [−630, +450] km s⁻¹ velocity range. The colored contours start at 2σ and increase by a factor of 1.5. The 1σ rms noise is 0.075, 0.095, and 0.082 Jy beam⁻¹ km s⁻¹ for F07599C, F07599S, and F07599NE, respectively. The CO map of F07599C was derived from a cleaned data cube, while the maps of F07599S and F07599NE were not cleaned given the relatively low S/N.</figcaption></figure>

and [−330, +270] km s⁻¹, where the emission is detected in consecutive channels at above 1σ, for F07599C, F07599S, and F07599NE, respectively. The velocity-integrated CO emission of each component was imaged separately by collapsing the velocity channels over the line-emitting region. It can be seen from Figure 4 that the superimposed image of the three components agrees well with the global CO map that is integrated in the [−630, +450] km s⁻¹ velocity range with coincident image contours (e.g., contours at a 3σ level). We derive an integrated flux density of 11.24 ± 0.09 Jy km s⁻¹, 1.12 ± 0.07 Jy km s⁻¹, and 0.45 ± 0.07 Jy km s⁻¹, by integrating the intensity over the line-emitting region in each channel of the spectra shown in Figure 4, for F07599C, F07599S, and F07599NE, respectively. From this analysis, the two subcomponents (F07599S and F07599NE) are identified at 16.0σ and 6.4σ in the intensity maps of the residual data cube, respectively. These weak features are also clearly present in the channel maps, where extended structures are found in the S and NE direction (see Appendix B).

Figure 5 shows the multiwavelength images of IRAS F07599+6508 in Sloan Digital Sky Survey (SDSS) *u*, *g*, *r*, *i*, and *z* and Pan-STARRS1 *y* bands and an overlay of the CO (1–0) emission of the three decomposed components on the Hubble Space Telescope (HST)/WFPC2 F702W optical image. The optical morphologies shown in *gri* and HST images are found to be complex with tidal-like features, indicative of an asymmetric stellar distribution, while the ultraviolet (UV; *u* band at ~0.35 µm) and near-IR (*z* and *y* bands at ~1 µm) emission shows more regular and compact distribution. A consistent morphology structure in the UV band has been presented by Surace & Sanders (2000). The tidal features seen in the optical are probably connected with an ongoing or past merging event in IRAS F07599+6508. Moreover, the absence of tidal tails observed in the near-IR bands is likely to suggest that the emission from optical tidal tails may trace the massive, young stellar population formed at recent epochs, since the near-IR emission is mostly sensitive to the longer-lived, low-mass stars. If this is indeed the case, it will also explain the morphology observed in the UV band, where the emission is largely attenuated in dust-obscured regions.
<figure><img src="image_6.png" /><figcaption>**Figure 5.** Left: optical images of IRAS F07599+6508. From left to right and top to bottom, we show the SDSS *u*-, *g*-, *r*-, *i*-, and *z*-band and Pan-STARRS1 *y*-band images. Right: overlay of CO (1–0) contours (colored) of the three components on the HST/WFPC2 F702W band image. The CO (1–0) contours are at the same levels as in Figure 4. The optical diffraction spikes in the middle of the object reflect the central luminous point source of an AGN that is significantly brighter than its host galaxy. Both the SDSS *gri* and HST optical images show clear signs of extended structure to the southwest and southeast with respect to the main source.</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>**Figure 6.** CO (1–0) velocity field (top) and velocity dispersion (bottom) maps for the whole galaxy IRAS F07599+6508 (first column) and the three subcomponents, F07599C (second column), F07599S (third column), and F07599NE (fourth column). The CO velocity and velocity dispersion maps are calculated using pixels detected at ≥3σ and are masked where the integrated intensity is <3σ in each component. The plus sign marks the peak of 3 mm continuum emission, while the two crosses indicate the CO peaks of F07599S and F07599NE, respectively. The synthesized beam is shown in the lower left corner of the first column panel.</figcaption></figure>

Interestingly, by comparing the distribution of CO molecular gas and stellar emission (see the right panel of Figure 5), we find that the tidal feature in the optical is coincident with the CO substructure F07599S. This provides evidence that the CO feature observed in the S region may be a real structure. Moreover, the CO emission also shows an extension to the SE, the same direction as the other tidal-tail-like feature observed in the optical images, although the CO extended structure toward the SE appears in only a few adjacent channels (e.g., channel maps at a velocity of about −178 and −58 km s⁻¹; see Appendix B). Meanwhile, a faint peak coincident with the location of the CO peak in F07599NE is identified in the HST image, indicative of a possible association.

### 3.2.3. CO Kinematics

In Figure 6 we show the velocity and velocity dispersion maps of CO (1–0) emission within the [−330, +330] km s⁻¹ velocity range for the whole galaxy IRAS F07599+6508 and the three decomposed components. These maps were created by applying a flux threshold of 3σ. The velocity map of the whole galaxy IRAS F07599+6508 shows irregular structure, indicative of complex velocity distribution. For the main component F07599C that is decomposed from the whole system, a systematic velocity gradient seems to be present in the north–south direction, from about −70 to +70 km s⁻¹, although the gas at the edge to the southeast and west appears
to be disturbed. This implies that the bulk of molecular gas is likely rotationally supported. However, the spatial resolution of our data is not adequate to unambiguously interpret the observed velocity gradient as due to disk rotation. An alternative possibility could be that a velocity gradient and broad-line profile may originate from an ongoing merger, as the velocity gradient is observed roughly along the direction connecting F07599C and F07599S. For F07599S, a velocity gradient of ~180 km s⁻¹ (from about −70 to +110 km s⁻¹) is also observed in the north–south direction, while no well-defined structure is visible in the velocity map of F07599NE, in which the velocity is positive around 60–90 km s⁻¹ with respect to the QSO redshift.

If we assume a rotating disk geometry for the main component F07599C, we can estimate the dynamical mass within the CO-emitting region as $M_{dyn}/M_{\odot} \approx 2.3 \times 10^5 v_{cir}^2 R$ (e.g., Wang et al. 2010), where R is the CO disk radius of 4.6 kpc (0.75 times the CO FWHM major axis from visibility fitting; see Table 1) and $v_{cir}$ is the maximum circular velocity of the gas disk in km s⁻¹. The $v_{cir}$ is estimated as $v_{cir} = 0.75\Delta v_{FWHM}/\sin i$ (Ho 2007), where i is the inclination angle between the polar axis of the disk and the line of sight and can be roughly calculated from the CO minor-to-major-axis ratio ($a_{min}/a_{maj}$; see Table 1), i.e., $i = \cos^{-1}(a_{min}/a_{maj}) = 45^\circ$. The derived $M_{dyn}$ is $2.1 \times 10^{11} M_{\odot}$, with a $v_{cir}$ of 450 km s⁻¹ for the primary disk F07599C. We stress that this estimate of dynamical mass can be significantly affected by the large uncertainties of source size, inclination angle, and the assumption of an inclined disk geometry. The total amount of molecular gas in the host galaxy is $M_{tot,mol} \sim 1.1 \times 10^{10} M_{\odot}$ by adopting a classic ULIRG-like CO-to-H₂ conversion factor of $\alpha_{CO} = 0.8 M_{\odot}$ (K km s⁻¹ pc⁻²)⁻¹. We obtain a molecular gas fraction of ~0.05 for IRAS F07599+6508 by assuming that the gas fraction can be approximated as the ratio of gas to dynamical mass, slightly lower than those (median value of 0.22 ± 0.04) of IR QSOs measured by Tan et al. (2019), and the local ULIRGs as well (Downes & Solomon 1998). Here we note the large uncertainties involved in the estimate of molecular gas mass, due to the assumption of conversion factor $\alpha_{CO}$, which is uncertain and varies significantly from source to source (Bolatto et al. 2013). Recent studies show evidence for higher $\alpha_{CO}$ values in (U)LIRGs that are affected by outflows if a substantial amount of molecular gas is in dense, gravitationally bound states (e.g., Cicone et al. 2018, and references therein).

The velocity dispersion map shows $\sigma_v \sim 160$ km s⁻¹ at the center of IRAS F07599+6508, which is significantly higher than the typical peak dispersion of ~80 km s⁻¹ observed in local IR QSOs (Tan et al. 2019), and $z \gtrsim 6$ QSOs (~130 km s⁻¹) as well (Neeleman et al. 2021). A plausible explanation could be related to the beam smearing effect that increases the observed dispersion, since the spatial angular resolution of our NOEMA observations is about an order of magnitude lower than the ALMA observations (~0.′5) shown by Tan et al. (2019). Assuming a uniform intrinsic dispersion across the galaxy disk of IRAS F07599+6508, we then estimate a dispersion of $\sigma_v \sim 70–100$ km s⁻¹, which is measured in the outer region of F07599C, where the beam smearing is considered to be less severe (Davies et al. 2011). The velocity dispersion map of F07599S shows an increase in the north direction with a peak dispersion of ~110 km s⁻¹, which is in the overlap region between F07599C and F07599S, while the gas in F07599NE appears less perturbed with a dispersion of ~30–60 km s⁻¹. The increased velocity dispersion observed in the overlap region

<figure><img src="image_8.png" /><figcaption>**Figure 7.** VLA 1.49 GHz contour (black) maps of IRAS F07599+6508 overlaid on the CO (1–0) contour (colored) maps of the three components. The CO (1–0) contours are at the same levels as in Figure 4. The contours of 1.49 GHz continuum emission start at $2\sigma$ (1σ noise level is 74.9 µJy beam⁻¹) and increase by a factor of 1.5. The solid orange and black ellipses in the lower left and right corners are the beam size of the CO and 1.49 GHz maps, with angular resolutions of 5.′′61 × 3.′′77 and 5.′′33 × 2.′′08, respectively.</figcaption></figure>

could be interpreted as a turbulence-dominated shock region induced by a merger or a possible CO outflow component. Higher-resolution observations are needed to better constrain the intrinsic velocity dispersion.

### 3.3. Radio Continuum Emission

In star-forming galaxies, the radio emission is dominated by synchrotron radiation from relativistic electrons accelerated in supernova remnants of massive stars. It is largely optically thin and unaffected by dust extinction and therefore is an excellent probe of very recent star formation activity in galaxies (Condon 1992). For the radio-quiet (RQ) AGNs, which represent the majority of the AGN population, the origin of radio emission could be connected with a wide range of possible mechanisms, namely, star formation, AGN-driven wind, jets, and accretion disk coronae (Panessa et al. 2019).

Figure 7 shows the VLA 1.49 GHz contour maps overlaid on the NOEMA CO (1–0) data that are decomposed into three components. Along the direction of the minor axis of the synthesized beam of 1.49 GHz observations, two clumpy structures outside the central bright compact radio source are potentially detected at the ~3σ and ~4σ level to the S and NE direction, respectively. The radio peaks of both the central compact source and the two clumps are well aligned with the CO peaks of the three decomposed components, reinforcing the fact that the faint CO features observed in the S and NE direction are very likely to be real.

A 2D Gaussian fit to the 1.49 GHz image returns a peak flux density of 32.33 ± 0.10 mJy beam⁻¹ and a consistent integrated flux density of 32.43 ± 0.19 mJy, indicative of a relatively compact structure of radio emission. This measurement agrees well with the 1.4 GHz flux density of 39.5 mJy derived from the NRAO VLA Sky Survey (Condon et al. 1998). IRAS F07599+6508 is additionally detected at 4.86, 8.44, 14.94, and 22.46 GHz by VLA. The maps and the flux measurements at these radio frequencies are shown in Appendix C. At redshift
<figure><img src="image_9.png" /><figcaption>**Figure 8.** Radio to millimeter SED of IRAS F07599+6508 based on the data from the VLA (circle) archive and NOEMA (square).</figcaption></figure>

$z=0.1486$, observations at the radio frequencies listed in Appendix C, Table 2 sample ~1.7, 5.6, 9.7, 17.2, and 25.8 GHz in the rest frame, respectively. In the following section, we adopt the rest-frame frequencies for clarity.

Figure 8 shows the SED of the continuum emission at radio and millimeter bands in IRAS F07599+6508. At GHz frequencies, the flux density $S_ν$ in terms of the observed frequency $ν$ is generally defined as $S_ν \propto ν^α$, with $α$ being the radio spectral index. The spectral indices measured from 1.7–5.6 GHz, 5.6–9.7 GHz, 9.7–17.2 GHz, and 17.2–25.8 GHz are $α_{5.6}^{1.7} = -1.10 \pm 0.09$, $α_{9.7}^{5.6} = -1.16 \pm 0.22$, $α_{17.2}^{9.7} = -1.18 \pm 0.26$, and $α_{25.8}^{17.2} = -0.34 \pm 0.78$, respectively. It is clear that the spectrum turns over from steep to flat at 17.2–25.8 GHz, which is consistent with the expectation that the contribution of thermal emission increases at higher frequencies. At low frequencies, the spectral index $α_{5.6}^{1.7}$ is slightly flatter than those at higher frequencies, $α_{9.7}^{5.6}$ and $α_{17.2}^{9.7}$, both of which show almost identical values. The flatter index observed at lower frequencies has been interpreted as a result of optical depth effects, as the radio spectrum at much lower frequencies typically shows a sharp decline due to free-free absorption (e.g., Condon 1992; Clemens et al. 2010). In contrast, the flux at larger frequency should be much less affected (e.g., Clemens et al. 2008). Therefore, we assume a synchrotron (nonthermal) spectral index of $α^{NT} = -1.18$, equivalent to the steepest part of the radio spectrum measured between 9.7 and 17.2 GHz (see the dashed line in Figure 8). This slope is slightly steeper than a normal synchrotron spectral index of $α = -0.8$ derived for a normal star-forming galaxy (Condon 1992), but consistent with the slopes measured for RQ AGNs (e.g., Behar et al. 2015).

We extrapolate the contribution of nonthermal emission from $α_{17.2}^{9.7}$ at the NOEMA frequency of 90.5 GHz (~103.9 GHz in the rest frame) and subtract it from the observed flux density. In addition, subtracting off an estimate of the thermal dust contribution to the 103.9 GHz, which is derived from a graybody model fit to the Hershel/SPIRE data with an emissivity spectral index $β = 1.6$ (Clemens et al. 2018), we obtain an estimate of the free-free emission component of

0.17 mJy. This implies a free-free radio fraction at 103.9 GHz of ~30%. A similar procedure applied to 25.8 GHz data leads to a free-free fraction of ~30% estimated at this frequency. This is consistent with the median thermal fraction of 20%–30% at 22 GHz measured for local ULIRGs (Vega et al. 2008).

We derive an extrapolated 1.4 GHz flux density of 43.4 mJy from the observed rest-frame 1.7 GHz flux density by assuming a radio power-law index of -1.18. Using the flux densities measured at IRAS 60 and 100 µm, we calculate an FIR-to-radio flux ratio<fnref n="11" /> $q = 1.67$. This implies a linear ratio about 4 times smaller than the average value of $q = 2.34 \pm 0.01$ found for local star-forming galaxies (Yun et al. 2001), possibly suggestive of an additional contribution from a compact AGN source that boosts the radio flux.

## 4. Discussion

The main results of our NOEMA observations are the detection of 3 mm continuum emission and the multiple spatially decomposed CO components in the host galaxy of IRAS F07599+6508. By combining with multiwavelength data, we begin the discussion by looking into the properties and the possible origin of these components.

### 4.1. Nature of the 3 mm Emission

The 3 mm continuum emission of IRAS F07599+6508 is resolved with a deconvolved size of $1''5 (~4.0 kpc)$. Assuming a nonthermal spectral index of $α^{NT} = -1.18$, we find that about 30% of the emission at 3 mm arises from a thermal free-free component. This is significantly lower than the thermal fraction of ~70% at 3 mm expected for normal star-forming galaxies (Condon 1992), indicating that the nonthermal synchrotron component from AGNs may contribute substantially to the 3 mm emission in IRAS F07599+6508. This could be expected since it does not follow the FIR-radio correlation characteristic of star-forming galaxies, with a ratio of $L_{FIR}/L_{1.4 GHz}$ being about four times lower. A similar situation has been observed in Mrk 231, the nearest IR QSO with a compact radio core (e.g., Lonsdale et al. 2003), where a rather low thermal fraction is observed at 33 GHz (Barcos-Muñoz et al. 2017). As mentioned above, an additional nonthermal component from AGNs could contribute to the emission at high radio frequencies. An alternative explanation for the observed low thermal fraction is connected with the possibility that the production of free-free emission might be suppressed owing to the absorption of ionizing stellar photons by dust that concentrates in the dense starburst regions.

Taking the 2–10 keV X-ray luminosity from Brightman & Nandra (2011), we find a ratio of radio (~3 mm) luminosity to X-ray luminosity of $L_{mm}/L_X \sim 3 \times 10^{-2}$ for IRAS F07599+6508. This value is about two orders of magnitude lower than those of radio-loud (RL) AGNs, which typically have powerful relativistic jets that shine in the radio band, but close to the ratios $(10^{-2}-10^{-4})$ measured for RQ AGNs (Behar et al. 2015; Panessa et al. 2019), suggestive of a different origin of radio emission compared with RL QSOs. We note that the ratio of $L_{mm}/L_X \sim 3 \times 10^{-2}$ derived here could be an overestimated value, given that the AGN embedded in the galaxy IRAS F07599+6508 is known to be heavily obscured and

<footnote n="11">¹¹ The FIR-to-radio flux ratio $q$ is defined as the ratio of the FIR luminosity ($40–120 \ \mu\text{m}$) to the radio power (Condon 1992): $q \equiv \log_{10}\left(\frac{L_{\text{FIR}}}{3.75 \times 10^{12} \ \text{W}}\right) - \log_{10}\left(\frac{L_{\text{1.4 GHz}}}{\text{W} \ \text{Hz}^{-1}}\right)$.</footnote>
thus a portion of X-ray photons may be absorbed (see Luo et al. 2014; La Caria et al. 2019). A connection between the radio emission and the nuclear accretion disk has been proposed to explain the correlation of the radio luminosity $L_{\text{mm}}$ and the X-ray luminosity $L_X$ observed in RQ AGNs (e.g., Laor & Behar 2008). More evidence for an X-ray corona (i.e., hot plasma with a temperature of ~$10^9$ K) origin of radio emission in RQ AGNs has been presented by the detections of millimeter excess (e.g., Behar et al. 2015; Inoue & Doi 2018). These are in line with the theoretical prediction for a radio synchrotron component from an X-ray corona in the vicinity of BH accretion disks contributing to the emission in the millimeter band (see Inoue & Doi 2014; Raginski & Laor 2016). Nevertheless, conclusive evidence for the coronal radio emission would be the detection of the radio/millimeter and X-ray variability correlation similar to that observed in stellar coronae (Panessa et al. 2019).

The thermal free-free emission arises directly from H II regions ionized by massive stars, and its intensity is proportional to the ionizing photon rate from young star-forming regions, thereby providing an excellent diagnostic for the current star formation rate (SFR) of galaxies (Condon 1992). Taking the free-free emission flux measured at rest frame 103.9 GHz and the free-free SFR calibration given in Murphy et al. (2012), we calculate an SFR of 77 $M_{\odot}$ yr$^{-1}$ for IRAS F07599+6508 by assuming a thermal spectral index of -0.1 and an electron temperature of $T_e = 10^4$ K. Compared to the IR-derived SFR of 280 $M_{\odot}$ yr$^{-1}$ that is estimated based on the heated dust component (Hao et al. 2008), a lower SFR measured by the free-free emission could be expected since free-free emission is sensitive to massive stars with ages $\lesssim 10$ Myr and consequently probes the very recent star formation activity. In addition, the free-free-derived SFR would be underestimated if a significant fraction of ionizing photons are absorbed by dust. Alternatively, it could be that a portion of IR luminosity arises from the dust heating by nuclear luminous AGNs, which would lead to an overestimate of SFR. Indeed, we obtain an SFR of 150 $M_{\odot}$ yr$^{-1}$ using the calibration given by Murphy et al. (2012) by correcting for the AGN contribution to the bolometric luminosity ($\alpha_{\text{AGN}} = 0.75$; Veilleux et al. 2009), which is close to the free-free SFR estimate. All of these, including the uncertainties (e.g., variations in the actual $T_e$) on the estimation of free-free-derived SFR, may explain the discrepancy in SFR estimates.

## 4.2. Origin of the Off-center CO Components

By exploiting the CO line cube of IRAS F07599+6508, we have identified two off-center CO subcomponents, F07599S and F07599NE, at a significance level of 16.0σ and 6.4σ in the intensity map, respectively. Figure 4 shows the CO (1–0) spectra and maps extracted for these two subcomponents and the main component F07599C as well. F07599S and F07599NE are located to the south and northeast of the disk component with a projected distance of ~11.4 and 19.1 kpc and with a fraction of 9% and 4% of CO flux contributed to the total CO flux, respectively. These two CO subcomponents are not detected in 3mm continuum emission. For the main component F07599C, the systematic CO velocity gradient observed may be a signature of a rotating disk structure. If confirmed, the source size of ~6.1 kpc measured would suggest that the CO disk is very extended, which is in sharp contrast with local ULIRGs, where the majority of molecular gas is

centrally concentrated within a radius of ~1 kpc (Downes & Solomon 1998; Wilson et al. 2008).

The optical images of IRAS F07599+6508 show evidence for a tidal-tail-like structure (see Figure 5), suggestive of a recent merging event in IRAS F07599+6508. It is interesting to note that the tidal feature shown in the optical image is spatially coincident with the CO feature detected in the S region. In addition, a faint peak coincident with the location of the CO peak in F07599NE is identified in the HST optical image, likely associated with the origin of the CO plume feature observed in the NE direction. It can be seen from the channel maps in Appendix B that the CO emission is extended to the S and NE direction, suggesting that the CO emission identified in the extranuclear region of F07599S and F07599NE is very likely to be a real component. Furthermore, the two potentially detected clumpy structures revealed in the 1.49 GHz radio image are found to align well with the peaks of F07599S and F07599NE (see Figure 7), as well as the faint substructures detected at 4.86 and 8.44 GHz (see Appendix C), supporting the reality of such a CO feature.

In Figure 4 we have shown that the CO (1–0) spectrum extracted for F07599S is best fitted with a two-component (narrow and broad) Gaussian. The detection of a narrow (FWHM = 100 ± 10 km s⁻¹) feature may indicate the presence of extended dynamical “quiescent” CO gas. For the broad component, the Gaussian fit gives an FWHM of 500 ± 80 km s⁻¹, slightly larger than the line FWHM (430 ± 10 km s⁻¹) of the primary component F07599C. Together with the shift in velocity (~100 km s⁻¹) and the tidal-tail-like structure observed in the optical image, it is likely that F07599S represents a separate component from the main galaxy and has a merger origin. One of the possibilities is that F07599S may trace the emission from a small obscured companion to the host galaxy of IRAS F07599+6508. The nondetection of UV emission in F07599S shows evidence for the large attenuation by dust (see Figure 5). More evidence for a merger origin is given by the increased velocity dispersion observed in the overlap region between F07599S and the main component, which could be interpreted as a turbulence-dominated shock region induced by an interaction.

For F07599NE, the best fit of the CO (1–0) line is a single Gaussian with both FWHM and peak velocity consistent within the errors with the main component. The significant lower velocity dispersion observed in this component implies that F07599NE may have a different origin from F07599S. Similar gas substructure has been observed in high-z starbursts in [C II] emission but with a much smaller physical separation (~2 kpc; Tadaki et al. 2020), where the formation of a subcomponent is expected to be associated with a gravitationally unstable disk or shocks. By contrast, the large physical separation (projected distance of ~19 kpc) observed between F07599NE and the galaxy nucleus indicates that the CO plume feature observed in the NE direction could be a separate component, likely associated with the tidal debris reminiscent of a previous merger. Similar results have been seen in the IR QSO IRAS F2341+0228, where an off-center clumpy structure is detected in both CO line and 3mm continuum maps and is found to align with a faint peak identified in the HST image, which was interpreted as either an ongoing merger or a clumpy star-forming region at large radii (Tan et al. 2019). Cicone et al. (2015) have shown a very extended (up to r ~ 30 kpc) [C II] cold gas component observed in a z > 6 QSO host galaxy, and they interpret that the extended [C II] emission may arise from a different interstellar medium phase than dense photon-dominated regions.
We note that there is a source accompanying IRAS F07599+6508 in the multiwavelength data set (from GALEX near-UV to WISE 3.4 µm bands; see Appendix D), with a distance of 17.5 (projected distance of ~45 kpc) to the northwest. This source has a cross-identification as WISEA J080428.28+650001.7 (hereafter J0804+6500) in the NASA/IPAC<sup>12</sup> Extragalactic Database. There are no optical spectroscopic data available for this object. In Appendix D, we show the NOEMA spectrum extracted from the position of J0804+6500 by fitting the visibilities with a point source. We do not detect any significant emission in either continuum or line from this source. With the data available, it is not possible to assess whether J0804+6500 is physically related to the IRAS F07599+6508 system, or is only a projected companion source. The detection of line emission from this source would help us measure the redshift and further determine whether it is an interacting companion of IRAS F07599+6508.

Together with the CO clumpy structures of local IR QSOs revealed with ALMA high-resolution observations (Tan et al. 2019), the multiple components revealed in the CO maps show evidence for the complex morphology that may be associated with the intense star formation and AGN activities. However, additional deep CO observations with higher angular resolution are needed for confirming the faint features observed in F07599S and F07599NE and testing the corresponding hypothesis.

### 4.3. Comparison of Molecular Gas with Mrk 231

As noted previously in Section 1, IRAS F07599+6508 and Mrk 231 share many common or similar properties in multi-wavelength, i.e., X-ray, UV, optical, and IR, and have been proposed as local templates of candidate ULIRG-to-optical QSO transition objects. High-resolution CO observations of Mrk 231 show that the molecular gas is concentrated in the central ~1 kpc rotating disk (Downes & Solomon 1998; Cicone et al. 2012; Feruglio et al. 2015), though O$_2$ emission might be present over ~10 kpc (Wang et al. 2020). Compared with Mrk 231, the molecular gas of IRAS F07599+6508 is also likely distributed in a disk with a systematic velocity gradient in the main component, but much more extended with a size of ~6.1 kpc. In addition, the CO line width of IRAS F07599+6508 is ~430 km s$^{-1}$ FWHM, much broader than that of 190 km s$^{-1}$ measured for the core component of Mrk 231 (Downes & Solomon 1998). The CO (2–1)/CO (1–0) line luminosity ratios derived for IRAS F07599+6508 and Mrk 231 are both ~0.8 (Xia et al. 2012; Cicone et al. 2012), indicative of modest subthermal molecular gas in their hosts.

The detection of massive, highly energetic, and kiloparsec-scale molecular outflows in Mrk 231 show direct observational evidence for AGN feedback onto the host galaxy (e.g., Feruglio et al. 2010, 2015; Aalto et al. 2012; Cicone et al. 2012), consistent with the AGN–galaxy coevolutionary models that propose that the outflow feedback is an important process in quenching the star formation activity in massive galaxies (e.g., Di Matteo et al. 2005; Hopkins et al. 2006). The CO (1–0) wing-to-peak ratio in Mrk 231 is ~2.5%, in agreement with the ratio of 1% to 8% typically measured for local ULIRGs (Feruglio et al. 2010; Cicone et al. 2012, 2014).

The presence of molecular outflow in IRAS F07599+6508 has been confirmed by the detection of an OH 119 µm P Cygni profile with velocities <1000 km s$^{-1}$ (Veilleux et al. 2013). A fast wind has also been seen in this system in the ionized gas as

traced by the broad blueshift optical spectra (Hα and Hβ) and in the neutral atomic phase with broad absorption observed in the Na I D line (e.g., Zheng et al. 2002; Rupke et al. 2017).

In the top left panel of Figure 9 we show the continuum-subtracted NOEMA CO (1–0) spectrum of the whole galaxy of IRAS F07599+6508, which is obtained by summing up the flux of the three decomposed components in each channel, compared with the CO (1–0) line observed with the IRAM 30 m (Xia et al. 2012). The CO peak flux measured by NOEMA is found to agree well with the IRAM 30 m data, suggesting that little emission is filtered out by our interferometric observations. The residual spectrum after removing a Gaussian fit to the core component of the CO spectrum is shown in the bottom left panel of Figure 9. A noise level (1σ) of 0.68 mJy in each ~60 km s$^{-1}$ channel was measured, corresponding to 2.3% of the peak flux density. Scaling to the signal strengths from Mrk 231 would lead to an estimated flux density of the CO broad emission of ~0.75 mJy for IRAS F07599+6508, comparable to the noise level of our NOEMA data.

It is interesting to note that several consecutive channels in the blueshifted side of the CO spectrum (see Figure 9, bottom left) show flux density higher than 1σ. We averaged the visibilities of the high-velocity component in the velocity range [−1340, −2000] km s$^{-1}$ and fitted a point-source model. The peak of high-velocity CO emission is located between F07599C and F07599S with a projected distance of 4.3 from the galaxy nucleus (see Figure 9, right). The flux derived from uv_fit is 0.228 ± 0.054 Jy km s$^{-1}$, corresponding to a detection significance of 4.2σ. However, this does not satisfy the criterion defined by Cicone et al. (2014), i.e., CO emission with velocities >500 km s$^{-1}$ is detected with a significance of at least 5σ, for claiming the detection of a molecular outflow. The OH outflow counterpart in CO emission is not significantly detected in our NOEMA data, although there is a hint of some flux at very high blueshifted velocities in the overlap region between F07599C and F07599S (see Figure 9), where the extracted CO spectrum shows a clear decomposition into broad plus narrow components (see the CO spectrum of F07599S in Figure 4). Nevertheless, it is difficult to rule out the existence of a weak CO outflow with the current data.

Given the much extended CO disk providing fuel for active large-scale star formation in IRAS F07599+6508, perhaps a much weaker outflow feedback in quenching star formation is expected here as compared to that of Mrk 231, where bursts of star formation only remain in the central compact core. If that is true, IRAS F07599+6508 should be in an earlier phase than that of Mrk 231 in the sense of ULIRGs’ evolution in transition into the classical QSOs. We stress that additional NOEMA CO observations are required to investigate this source in detail, particularly the presence of a possible CO outflow.

## 5. Summary

In this work, we present the deep NOEMA interferometric observations of the CO (1–0) line and 3 mm continuum emission in the ultraluminous IR QSO IRAS F07599+6508 at $z=0.1486$. This system shares many physical properties with Mrk 231, the nearest and most luminous IR QSO, in multiwavelength. Both the CO line and 3 mm continuum emission are resolved in our observations, and the CO gas is found to be spatially decomposed into three components. Combining these observations with results from multiwavelength data (i.e., radio, optical, and infrared), we study the properties of the 3 mm continuum emission and the spatially decomposed CO components, and the main results are summarized as follows:

<footnote n="12">12 http://ned.ipac.caltech.edu</footnote>
<figure><img src="image_10.png" /><figcaption>**Figure 9.** Left: continuum-subtracted spectra of the CO (1–0) emission line (top) and the residuals (bottom) after removing a Gaussian fit (red line in the top panel) to the CO spectrum of IRAS F07599+6508. In the top panel, the NOEMA CO spectrum (black) is obtained by summing up the flux of the three decomposed components in each channel (see the left panel of Figure 4), while the IRAM 30 m CO spectrum (blue) is also plotted for a comparison. In the bottom panel, the dashed lines denote the $1\sigma$ level. The spectra are binned to a velocity resolution of ~60 km s$^{-1}$. Right: contours from CO (1–0) dirty maps integrated over the velocity range of [−1340, −2000] km s$^{-1}$. Contours start at $\pm 2\sigma$ and increase in steps of $1\sigma$ (81 µJy beam$^{-1}$), with positive (negative) contours shown as solid (dotted) lines. The plus sign and cross indicate the CO peak positions of F07599C and F07599S, respectively. The synthesized beam is shown in the lower left corner.</figcaption></figure>

1. The 3 mm continuum emission is resolved with a deconvolved size of $1''5$ (~4.0 kpc). Assuming a nonthermal spectral index of $\alpha^{NT} = -1.18$, equal to the slope measured between 9.7 and 17.2 GHz representing the steepest part of the radio spectrum, we find that ~30% of the emission at 3 mm arises from a free-free component, significantly lower than the fraction (~70%) expected for a normal star-forming galaxy. This suggests that the nonthermal synchrotron component from AGNs may contribute substantially to the 3 mm emission. The relatively low FIR-to-radio flux ratio of $q = 1.67$ derived also shows evidence for an additional contribution from the central AGN that is boosting the radio flux. We estimate an SFR of $77 M_{\odot}yr^{-1}$ based on the thermal free-free emission component at 3 mm, which is close to the IR estimate of $150 M_{\odot}yr^{-1}$ corrected for the AGN contribution to the bolometric luminosity.

2. Two off-center CO subcomponents, F07599S and F07599NE, located to the south and northeast of the main component F07599C with a projected distance of ~11.4 and 19.1 kpc, respectively, were identified by model fits to the u-v data. More evidence for the reality of F07599S and F07599NE is given by the spatial coincidence of the tidal-tail-like feature and the faint peak in the optical images and the potentially detected clumpy structures in the 1.49 GHz radio image with these two components. The fraction of F07599S and F07599NE contributed to the total CO flux is 9% and 4%, respectively. Analysis of the CO kinematics, together with the multiwavelength data for the decomposed components, suggests that the gas in the dominant component F07599C is likely rotationally supported with a systematic velocity gradient observed, while F07599S may represent a separate component associated with a merger and F07599NE probably has a different origin.

3. Careful analysis of our current CO (1–0) data indicates that the significance of the high-velocity emission feature seen on the blueshifted side of the spectrum is not sufficient for claiming the detection of a CO outflow in IRAS F07599+6508. It should be noted, however, that the spatially broad CO core distribution with a

size larger than 6 kpc is in marked contrast with the compact ~kiloparsec CO concentration in Mrk 231, where strong molecular outflows were observed over wide velocity ranges.

Additional deeper observations with higher angular resolution are required to better constrain the CO morphology and kinematics on kiloparsec and subkiloparsec scales, confirm the presence of high-velocity gas, and further determine whether the features seen here are potential signs of outflowing gas that is escaping from the host galaxy.

We thank the anonymous referee for the constructive comments that helped improve the paper. We acknowledge the support from IRAM staff during the observations and data reduction. This work is based on observations carried out under project No. S18CA with the IRAM NOEMA Interferometer. IRAM is supported by INSU/CNRS (France), MPG (Germany), and IGN (Spain).

This work was supported by National Key Basic Research and Development (R&D) Program of China (grant No. 2017YFA0402704); NSFC grant Nos. 11803090, 11861131007, 12033004, and 11420101002; and Chinese Academy of Sciences (CAS) Key Research Program of Frontier Sciences (grant No. QYZDJ-SSW-SLH008). X.X. and C.H. acknowledge support from the NSFC (grant No. 11733002). The work by K.K is supported by the JSPS KAKENHI grant No. JP17H06130 and the NAOJ ALMA Scientific Research grant No. 2017-06B.

**Facilities:** IRAM:NOEMA, VLA, Sloan, HST, PS1.

**Software:** GILDAS (Guilloteau & Lucas 2000), CASA (McMullin et al. 2007), AIPS (van Moorsel et al. 1996).

## Appendix A
### Decomposition of CO Emission

To illustrate the spatial decomposition of CO emission in IRAS F07599+6508, we here describe the model-fitting procedure. We fit the model to the u-v data using combinations of functions by allowing all parameters to be free. Combinations of
models in a successive progression of complexity started with a
point source, then a circular Gaussian, and then an elliptical
Gaussian. The simultaneous fit was rerun each time as necessary to
achieve a uniform residual (i.e., no significant peak or negative
regions $\gtrsim 3\sigma$) after subtraction of the models. By verifying the
evolution of the uncertainties on the fitted parameters, we found
that the best fit was obtained with a combination of an elliptical
Gaussian (F07599C) and two point-source (F07599S and
F07599NE) functions with minimum uncertainties on the fitted
parameters. The residuals after subtracting an elliptical Gaussian
plus two point-source components are shown in the right panel of
Figure 3. It should be noted that there are a few peaks with S/

N ~ 3 (e.g., the peaks to the northwest of F07599S and F07599NE) in the residual map, indicative of the likely presence of extended structure. These weak features are, however, difficult to test with the current data. We should stress that higher spatial resolution observations are needed to confirm these findings.

Appendix B
Channel Maps of CO (1–0) Emission

The spatial and velocity distribution of CO (1–0) emission
in IRAS F07599+6508 can be further explored by inspecting
the channel maps in Figure 10, which shows the image of CO
emission within an interval of ~120 km s⁻¹ encompassing the

<figure><img src="image_11.png" /><figcaption>**Figure 10.** Channel maps showing the CO (1–0) emission of IRAS F07599+6508 in the [−2000, +1700] km s⁻¹ range with channel width of ~120 km s⁻¹. The red crosses mark the positions of CO subcomponents (F07599S and F07599NE), and the black plus sign indicates the peak of 3 mm continuum emission. Contour levels start at ±1.5σ and increase by factors of 1.5 with positive (negative) contours shown as solid (blue dotted) lines. The velocities in km s⁻¹ are marked in the upper left corner of each channel map, and the synthesized beam is shown in the lower left corner.</figcaption></figure>
velocity range from −2090 to 2093 km s⁻¹. The channel maps were produced from the continuum-subtracted data cube that was cleaned down to the 1σ noise level. The extended structures to the south and northeast of galaxy center observed in each channel map, spatially coincident with the subcomponents of F07599S and F07599NE (marked as red crosses in Figure 10), show evidence for the reality of such features.

## Appendix C
Radio Maps and Flux Measurements

Table 2 summarizes the radio continuum properties of IRAS F07599+6508 obtained with VLA. Except the 22.46 GHz data that were obtained in a single CD configuration, observations at the rest of frequencies were made with multiple arrays.

To investigate the distribution of radio emission and make a comparison with CO emission, we show the map obtained with

higher angular resolution for each frequency in Figure 11. For 14.94 and 22.46 GHz, the emission is detected with a compact morphology, while for the emission at 4.86 and 8.44 GHz, where the observations were deeper, the radio source exhibits more complex morphology with clumpy structures outside the main body of the galaxy to the east and west, respectively. The peak flux densities measured in different configurations were found to be in good agreement for each frequency (i.e., to within ~10%; see Table 2). The highest-resolution image was obtained at 17.2 GHz with a beam size of ~0.′′18 (~0.47 kpc), of which the peak flux density measured is found to be comparable to that measured from CD configuration (beam size ~4.′′58). The consistent peak flux densities measured between different beams for the different bands indicate that the bulk of radio emission probably originates in the inner compact region of the AGN.

<table><caption>**Table 2**
VLA Radio Continuum Properties of IRAS F07599+6508</caption><thead><tr><th>Obs. freq.<br>(GHz)</th><th>R.A.<sup>a</sup><br>(J2000)</th><th>Decl.<sup>a</sup><br>(J2000)</th><th>Conf.</th><th>Synt. Beam<br>(arcsec)</th><th>S<sup>b</sup><sub>peak</sub><br>(mJy beam<sup>-1</sup>)</th></tr></thead><tbody><tr><td rowspan="2">1.49 GHz</td><td rowspan="2">08:04:30.53</td><td rowspan="2">64:59:53.3</td><td>AB</td><td>5.33 × 2.08</td><td>32.33 ± 0.10</td></tr><tr><td>CD</td><td>65.2 × 28.3</td><td>36.72 ± 0.16</td></tr><tr><td rowspan="2">4.86 GHz</td><td rowspan="2">08:04:30.51</td><td rowspan="2">64:59:52.6</td><td>AB</td><td>1.84 × 0.94</td><td>9.59 ± 0.05</td></tr><tr><td>CD</td><td>20.8 × 8.8</td><td>9.15 ± 0.16</td></tr><tr><td rowspan="2">8.44 GHz</td><td rowspan="2">08:04:30.51</td><td rowspan="2">64:59:52.6</td><td>AB</td><td>1.10 × 0.47</td><td>4.67 ± 0.07</td></tr><tr><td>CD</td><td>11.0 × 5.1</td><td>5.23 ± 0.09</td></tr><tr><td rowspan="2">14.94 GHz</td><td rowspan="2">08:04:30.471</td><td rowspan="2">64:59:52.76</td><td>A</td><td>0.18 × 0.13</td><td>2.38 ± 0.06</td></tr><tr><td>CD</td><td>4.58 × 2.56</td><td>2.69 ± 0.08</td></tr><tr><td>22.46 GHz</td><td>08:04:30.46</td><td>64:59:52.8</td><td>CD</td><td>2.92 × 1.93</td><td>2.20 ± 0.10</td></tr></tbody></table>

**Notes.**

<footnote><sup>a</sup> Source position derived from a 2D Gaussian fit (using the *imfit* task in CASA) to the high-resolution images.</footnote>

<footnote><sup>b</sup> Peak flux density measured from *imfit* using a single 2D Gaussian component.</footnote>
<figure><img src="image_12.png" /><figcaption>**Figure 11.** VLA multiband radio continuum imaging of IRAS F07599+6508 at 4.86, 8.44, 14.94, and 22.46 GHz from left to right. While the observations at 22.46 GHz were conducted in a single array of CD configuration, the rest of the observations were made with arrays in both AB configuration (A configuration for 14.94 GHz) and CD configuration. Here we only show the maps with higher angular resolution. All the contours start at 3σ and increase by a factor of 1.5. The synthesized beam and beam size of each radio observation are shown at the bottom of each panel. The white plus sign indicates the peak of 3 mm continuum emission.</figcaption></figure>

# Appendix D

## Multiwavelength Data of IRAS F07599+6508 and a Companion Source to the Northwest

In this appendix we show the multiwavelength (from GALEX near-UV to WISE 3.4 µm band) data of IRAS F07599+6508 and a companion source WISEA J080428.28 +650001.7 to the northwest (Figure 12). J0804+6500 has a distance of 17"5 (~45 kpc) from IRAS F07599+6508, with coordinates of 08:04:28.140, +65:00:02.25.

Figure 13 shows the NOEMA spectra (without correction for primary beam attenuation (PBA)) extracted from the position of J0804+6500 by fitting a point-source model to the *u*-(*v* data. The tuning frequency covered in the USB and LSB corresponds to a redshift range of 0.128–0.223 and 0.238–0.462, respectively, if associated with a CO (1–0) transition. The 1σ rms noise level of the spectrum in the USB and LSB shown in Figure 13 is 0.20 and 0.21 mJy, respectively. We do not detect any significant emission from J0804+6500 in our observations.
<figure><img src="image_13.png" /><figcaption>**Figure 12.** Image cutouts (60" × 60" on a side) of IRAS F07599+6508 and a companion source WISEA J080428.28+650001.7 to the northwest. From left to right, we show the GALEX near-UV, SDSS *u*- and *r*-, Pan-STARRS1 *y*-, 2MASS *Ks*-, and WISE 3.4 µm band images.</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>**Figure 13.** NOEMA spectrum (without correction for PBA) of WISEA J080428.28+650001.7 covered in the USB (top) and the LSB (bottom). The spectra are extracted by fitting a point-source model to the *u*-*v* data and are binned to a channel width of ~120–140 km s⁻¹.</figcaption></figure>

## ORCID iDs

Qing-Hua Tan https://orcid.org/0000-0003-3032-0948

Yu Gao https://orcid.org/0000-0003-0007-2197

Emanuele Daddi https://orcid.org/0000-0002-3331-9590

Cai-Na Hao https://orcid.org/0000-0002-0901-9328

Alain Omont https://orcid.org/0000-0002-4721-3922

Kotaro Kohno https://orcid.org/0000-0002-4052-2394

## References

* Aalto, S., Garcia-Burillo, S., Muller, S., et al. 2012, A&A, 537, A44

* Barcos-Muñoz, L., Leroy, A. K., Evans, A. S., et al. 2017, ApJ, 843, 117

* Barnes, J. E., & Hernquist, L. 1992, ARA&A, 30, 705

* Behar, E., Baldi, R. D., Laor, A., et al. 2015, MNRAS, 451, 517

* Bolatto, A. D., Wolfire, M., & Leroy, A. K. 2013, ARA&A, 51, 207

* Boroson, T. A., & Meyers, K. A. 1992, ApJ, 397, 442

* Brightman, M., & Nandra, K. 2011, MNRAS, 414, 3084

* Bryant, P. M., & Scoville, N. Z. 1996, ApJ, 457, 678

* Canalizo, G., & Stockton, A. 2001, ApJ, 555, 719

* Cao, C., Xia, X. Y., Wu, H., et al. 2008, MNRAS, 390, 336

* Cicone, C., Feruglio, C., Maiolino, R., et al. 2012, A&A, 543, A99

* Cicone, C., Maiolino, R., Aalto, S., Muller, S., & Feruglio, C. 2020, A&A, 633, A163

* Cicone, C., Maiolino, R., Gallerani, S., et al. 2015, A&A, 574, A14

* Cicone, C., Maiolino, R., Sturm, E., et al. 2014, A&A, 562, A21

* Cicone, C., Severgnini, P., Papadopoulos, P. P., et al. 2018, ApJ, 863, 143

* Clemens, M. S., Scaife, A., Vega, O., & Bressan, A. 2010, MNRAS, 405, 887

* Clemens, M. S., Vega, O., Bressan, A., et al. 2008, A&A, 477, 95

* Clements, D. L., Pearson, C., Farrah, D., et al. 2018, MNRAS, 475, 2097

* Condon, J. J. 1992, ARA&A, 30, 575

* Condon, J. J., Cotton, W. D., Greisen, E. W., et al. 1998, AJ, 115, 1693

* Daddi, E., Bournaud, F., Walter, F., et al. 2010, ApJ, 713, 686

* Davies, R., Förster Schreiber, N. M., Cresci, G., et al. 2011, ApJ, 741, 69

* Di Matteo, T., Springel, V., & Hernquist, L. 2005, Natur, 433, 604

* Downes, D., & Solomon, P. M. 1998, ApJ, 507, 615

* Evans, A. S., Frayer, D. T., Surace, J. A., & Sanders, D. B. 2001, AJ, 121, 1893

* Evans, A. S., Solomon, P. M., Tacconi, L. J., Vavilkin, T., & Downes, D. 2006, AJ, 132, 2398

* Fabian, A. C. 2012, ARA&A, 50, 455

* Feruglio, C., Fiore, F., Carniani, S., et al. 2015, A&A, 583, A99

* Feruglio, C., Maiolino, R., Piconcelli, E., et al. 2010, A&A, 518, L155

* Fischer, J., Sturm, E., González-Alfonso, E., et al. 2010, A&A, 518, L41

* Greve, T. R., Bertoldi, F., Smail, I., et al. 2005, MNRAS, 359, 1165

* Guilloteau, S., & Lucas, R. 2000, in ASP Conf. Ser., 217, Imaging at Radio through Submillimeter Wavelengths, ed. J. G. Mangum & S. J. E. Radford (San Francisco, CA: ASP), 299
Hao, C. N., Xia, X. Y., Mao, S., Wu, H., & Deng, Z. G. 2005, ApJ, 625, 78

Hao, C.-N., Xia, X.-Y., Shu-DeMao, D. Z.-G., & Wu, H. 2008, ChJAA, 8, 12

Ho, L. C. 2007, ApJ, 668, 94

Hopkins, P. F., Hernquist, L., Cox, T. J., et al. 2006, ApJS, 163, 1

Hopkins, P. F., Hernquist, L., Cox, T. J., & Kereš, D. 2008, ApJS, 175, 356

Imanishi, M., Nakanishi, K., & Izumi, T. 2019, ApJS, 241, 19

Inoue, Y., & Doi, A. 2014, PASJ, 66, L8

Inoue, Y., & Doi, A. 2018, ApJ, 869, 114

Jarvis, M. E., Harrison, C. M., Mainieri, V., et al. 2020, MNRAS, 498, 1560

Kim, D. C., & Sanders, D. B. 1998, ApJS, 119, 41

King, A., & Pounds, K. 2015, ARA&A, 53, 115

Kormendy, J., & Ho, L. C. 2013, ARA&A, 51, 511

Kwan, J., Cheng, F.-Z., Fang, L.-Z., Zheng, W., & Ge, J. 1995, ApJ, 440, 628

La Caria, M. M., Vignali, C., Lanzuisi, G., Gruppioni, C., & Pozzi, F. 2019, MNRAS, 487, 1662

Laor, A., & Behar, E. 2008, MNRAS, 390, 847

Lawrence, A., Rowan-Robinson, M., Ellis, R. S., et al. 1999, MNRAS, 308, 897

Le Floc'h, E., Papovich, C., Dole, H., et al. 2005, ApJ, 632, 169

Lonsdale, C. J., Lonsdale, C. J., Smith, H. E., & Diamond, P. J. 2003, ApJ, 592, 804

Luo, B., Brandt, W. N., Alexander, D. M., et al. 2014, ApJ, 794, 70

McMullin, J. P., Waters, B., Schiebel, D., Young, W., & Golap, K. 2007, in ASP Conf. Ser., 376, Astronomical Data Analysis Software and Systems XVI, ed. R. A. Shaw, F. Hill, & D. J. Bell (San Francisco, CA: ASP), 127

Murphy, E. J., Bremseth, J., Mason, B. S., et al. 2012, ApJ, 761, 97

Neeleman, M., Novak, M., Venemans, B. P., et al. 2021, ApJ, 911, 141

Panessa, F., Baldi, R. D., Laor, A., et al. 2019, NatAs, 3, 387

Raginski, I., & Laor, A. 2016, MNRAS, 459, 2082

Rupke, D. S. N., Gültekin, K., & Veilleux, S. 2017, ApJ, 850, 40

Rupke, D. S. N., & Veilleux, S. 2011, ApJL, 729, L27

Sanders, D. B., & Mirabel, I. F. 1996, ARA&A, 34, 749

Sanders, D. B., Soifer, B. T., Elias, J. H., et al. 1988a, ApJ, 325, 74

Sanders, D. B., Soifer, B. T., Elias, J. H., Neugebauer, G., & Matthews, K. 1988b, ApJL, 328, L35

Scoville, N. Z., Frayer, D. T., Schinnerer, E., & Christopher, M. 2003, ApJL, 585, L105

Shangguan, J., Ho, L. C., Bauer, F. E., Wang, R., & Treister, E. 2020, ApJS, 247, 15

Silk, J., & Rees, M. J. 1998, A&A, 331, L1

Soifer, B. T., Sanders, D. B., Madore, B. F., et al. 1987, ApJ, 320, 238

Solomon, P. M., Downes, D., Radford, S. J. E., & Barrett, J. W. 1997, ApJ, 478, 144

Spoon, H. W. W., Farrah, D., Lebouteiller, V., et al. 2013, ApJ, 775, 127

Surace, J. A., & Sanders, D. B. 2000, AJ, 120, 604

Tadaki, K.-i., Iono, D., Yun, M. S., et al. 2020, ApJ, 889, 141

Tan, Q.-H., Gao, Y., Kohno, K., et al. 2019, ApJ, 887, 24

Taylor, G. B., Silver, C. S., Ulvestad, J. S., & Carilli, C. L. 1999, ApJ, 519, 185

Ueda, J., Iono, D., Yun, M. S., et al. 2014, ApJS, 214, 1

van Moorsel, G., Kemball, A., & Greisen, E. 1996, in ASP Conf. Ser., 101, Astronomical Data Analysis Software and Systems V, ed. G. H. Jacoby & J. Barnes (San Francisco, CA: ASP), 37

Vega, O., Clemens, M. S., Bressan, A., et al. 2008, A&A, 484, 631

Veilleux, S., Kim, D. C., & Sanders, D. B. 1999, ApJ, 522, 113

Veilleux, S., Meléndez, M., Sturm, E., et al. 2013, ApJ, 776, 27

Veilleux, S., Rupke, D. S. N., Kim, D. C., et al. 2009, ApJS, 182, 628

Véron-Cetty, M. P., Joly, M., Véron, P., et al. 2006, A&A, 451, 851

Wang, J., Li, D., Goldsmith, P. F., et al. 2020, ApJ, 889, 129

Wang, R., Carilli, C. L., Neri, R., et al. 2010, ApJ, 714, 699

Wilson, C. D., Petitpas, G. R., Iono, D., et al. 2008, ApJS, 178, 189

Xia, X.-Y., Boller, T., Deng, Z.-G., & Börner, G. 2001, ChJAA, 1, 221

Xia, X. Y., Gao, Y., Hao, C.-N., et al. 2012, ApJ, 750, 92

Xu, C. K., Cao, C., Lu, N., et al. 2014, ApJ, 787, 48

Xu, C. K., Cao, C., Lu, N., et al. 2015, ApJ, 799, 11

Yun, M. S., Reddy, N. A., & Condon, J. J. 2001, ApJ, 554, 803

Zheng, X. Z., Xia, X. Y., Mao, S., Wu, H., & Deng, Z. G. 2002, AJ, 124, 18