Dust and Cold Gas Properties of Starburst HyLIRG Quasars at z ∼ 2.5

Some high-z active galactic nuclei (AGNs) are found to reside in extreme star-forming galaxies, such as hyperluminous infrared galaxies (HyLIRGs), with AGN-removed L IR of >1013 L ⊙. In this paper, we report NOEMA observations of six apparent starburst HyLIRGs associated with optical quasars at z ∼ 2–3 in the Stripe 82 field, to study their dust and molecular CO properties. Five out of the six candidates are detected with CO(4–3) or CO(5–4) emission, and four in the 2 mm dust continuum. Based on the linewidth– LCO(1−0)′ diagnostics, we find that four galaxies are likely unlensed or weakly lensed sources. The molecular gas mass is in the range of μMH2∼0.8–9.7×1010M⊙ (with α=0.8M⊙Kkms−1pc2−1 , where μ is the unknown possible gravitational magnification factor). We fit their spectral energy distributions, after including the observed 2 mm fluxes and upper limits, and estimate their apparent (uncorrected for possible lensing effects) star formation rates (μSFRs) to be ∼400–2500 M ⊙ yr−1, with a depletion time of ∼20–110 Myr. We notice interesting offsets, of ∼10–40 kpc spatially or ∼1000–2000 km s−1 spectroscopically, between the optical quasar and the millimeter continuum or CO emissions. The observed velocity shift is likely related to the blueshifted broad-emission-line region of quasars, though mergers or recoiling black holes are also possible causes, which can explain the spatial offsets and the high intrinsic star formation rates in the HyLIRG quasar systems.


INTRODUCTION
Active galactic nuclei (AGNs) represent a crucial phase in the evolution of supermassive black holes (SMBHs), and may strongly influence the evolution of their host galaxies (Fabian 2012;Kormendy & Ho 2013).Observations have shown that the powerful phase of AGN, i.e. quasars, influences host galaxies through either radiation pressure (e.g.Laor & Draine 1993;Scoville & Norman 1995) or AGN wind (e.g.Weymann et al. 1991;Pounds et al. 2003;Tombesi et al. 2012) during the so-called "quasar mode" feedback.This is consistent with the evolutionary model of quasars by Sanders et al. (1988), where quasars develop from dusty ultra-luminous infrared galaxies (ULIRGs) with infrared (IR) luminosity L IR > 10 12 L ⊙ .As the AGN feedback swepdf the gas and dust in the core region, the central AGN gets exposed in the line of sight and appears to be a type-I broad-emission-line quasar.The star formation (SF) in the host galaxy is suppressed during the process.
This model indicates a transitional stage in quasar evolution where a quasar coexists with a large amount of IRluminous galactic dust.In observation, there are 10-30% quasars with bright submillimeter/far-infrared (SFR) excesses (e.g.Dai et al. 2012Dai et al. , 2018;;Ma & Yan 2015;Dong & Wu 2016).Some of them have the most IR-luminous host galaxies, i.e., starburst hyper-luminous infrared galaxies (HyLIRGs) with starburst-dominated L IR > 10 13 L ⊙ .With IR-traced SF rate (SFR) of ≳ 10 3 M ⊙ yr −1 (e.g.Casey et al. 2012;Ivison et al. 2013;Banerji et al. 2013), they link the most powerful AGNs and the most extreme SF activities in the host galaxy.
However, the luminosity of these galaxies is questionable because of their potential gravitational magnification.Many apparent HyLIRGs have been found to be the results of lensing by large-area millimeter/submillimeter surveys (e.g.Vieira et al. 2010;Negrello et al. 2010;Wardlow et al. 2013;Bussmann et al. 2013;Cañameras et al. 2018) and follow-up studies (e.g.Yang et al. 2017;Zhang et al. 2018).For instance, the Herschel Astrophysical Terahertz Large Area Survey (H-ATLAS, Eales et al. 2010) revealed ∼1000 strongly lensed sources (González-Nuevo et al. 2012).After correcting for the lensing effect, many sources turn out to be of lower IR luminosities, thus no longer HyLIRGs.As a case, Timmons et al. (2016) found that the apparent HyLIRG HATLAS J132427 is an intrinsic ULIRG, after correcting for the magnification factor of five.Another apparent HyLIRG SDP.81 has a magnification factor of ∼18, reconstructed with ALMA data, and is indeed a ULIRG with an SFR of ∼ 100M ⊙ yr −1 (Rybak et al. 2020).
The high apparent IR luminosity can also be an effect of a collection of IR-bright sources.Luminous infrared galaxies (LIRGs, L IR > 10 11 L ⊙ ) have been extensively observed in the merging process (e.g., Chapman et al. 2003;Tacconi et al. 2006;Bothwell et al. 2010;Engel et al. 2010;Ivison et al. 2011Ivison et al. , 2013;;Riechers et al. 2011), which is also predicted by simulations (e.g., Swinbank et al. 2008;Narayanan et al. 2010;Hayward et al. 2011Hayward et al. , 2012;;McAlpine et al. 2019).Stateof-the-art telescopes with sub-arcsecond resolution powers, such as the Atacama Large Millimeter/submillimeter Array (ALMA) and the Very Large Array (VLA), have resolved some apparent HyLIRGs into multiple sources, confirming the resolved sources to be LIRGs or ULIRGs instead.For example, Fu et al. (2013) resolved two merging ULIRGs in the apparent HyLIRG 1HERMES S250 J022016.5-060143, which used to be considered as an unusually bright HyLIRG in the Herschel Multi-tiered Extragalactic Survey (HerMES, Oliver et al. 2012).
CO observations are crucial to study the physical properties of these galaxies.Firstly, it traces the immediate star-forming material in the host galaxy.This is connected to the feedback of AGN, as controversial results have been reported on whether they drive out gas or accelerate star formation efficiency (SFE, e.g., Kirkpatrick et al. 2019;Bischetti et al. 2021).Besides, it can reveal the direct feedback from quasars in the form of galactic scale outflows, which are exhibited as broad line wings exceeding a velocity of 500 km s −1 (e.g.Feruglio et al. 2010;Cicone et al. 2012).Finally, Harris et al. (2012) has found that strongly lensed galaxies can be distinguished by CO line emission.Thus, we are able to estimate the lensing property of galaxies under limited resolution.
In this work, we conduct NOrthern Extended Millimeter Array (NOEMA) observations of the mid-J CO rotational emission (J=4-3 or J=5-4) and 2mm dust continuum in six starburst apparent HyLIRG-quasars at z ∼ 2.5.The sample is selected from the quasar catalog of the Sloan Digital Sky Survey (SDSS) in the Stripe 82 field, with apparent HyLIRG-level IR luminosity from Herschel observations (Dong & Wu 2016).We use the millimeter observations to probe the molecular gas and dust properties in these sources and to identify if they are intrinsic or lensed starburst HyLIRG-quasars.
The paper is organized as follows: In section 2, we describe the sample selection, the observations, and the data reduction process.In section 3, the observational results are presented, including the continuum and CO emission morphology and properties, and the optical-to-mm spectral energy distributions (SEDs) of our sample.In section 4, we discuss the sample's location in the HyLIRG diagnostics, the spatial and velocity offsets of the sample between various tracers, and the estimated SFR and depletion time, followed by a summary in section 5.  , 58107, 58098, 52933, 58079, and 56979, respectively). b The Herschel coordinates and flux densities at 250µm (S250µm), 350µm (S350µm), and 500µm (S500µm) are taken from the HerS catalog (Viero et al. 2014).
c Integrated between 8-1000 µm derived from SED fitting by Dong & Wu (2016) considering only the starburst grey-body (i.e. cold dust) component.Note that we have derived improved values for these luminosities, µLIR,SB, in Table 7.
The virial BH masses estimated from broad C IV lines (Shen et al. 2011;Dong & Wu 2016).
e SDSS spectroscopic redshifts, which are dominated by broad C IV lines (see Shen et al. 2011;Lyke et al. 2020, and discussion in Section 4.3).
We select six apparent starburst HyLIRG-quasars from the catalog in Dong & Wu (2016).This sample was selected in the Stripe 82 field from the SDSS quasar catalogs (Schneider et al. 2010;Shen et al. 2011;Pâris et al. 2014) , which were pre-selected to be brighter than M i = −22.0 and have at least one optical line with full width at half-maximum (FWHM) larger than 1000 km s −1 (Type 1).These quasars were then cross-matched with the Herschel Stripe 82 Survey (HerS, Viero et al. 2014) by Dong & Wu (2016), and 207 showed Herschel SPIRE detections at 250, 350, and 500 µm.These wavelengths cover the spectral regions close to the peak of the cold dust emission at z ∼ 2.5 and thus can better constrain the cold dust properties in the SED fitting.AGN-subtracted IR luminosity L IR,8−1000µm was then calculated based on a grey-body dust component (Dong & Wu 2016).
We then selected sources with L IR,8−1000µm > 10 13 L ⊙ (Table 1), i.e., apparent starburst HyLIRGs.To avoid possible gravitational lensing and blending issues, we further required the targets to be point sources in the SDSS images, without any close companion within 5 ′′ , which is slightly larger than the NOEMA resolution (∼ 4 ′′ ) (D configuration).This way, we selected six starburst HyLIRG-quasars from the Dong & Wu (2016) catalog, namely DW001 to DW006.All of the selected sources have spectroscopic redshifts between 2 and 3, corresponding to the peaks of cosmic evolution for both star formation and AGN accretion (Förster Schreiber & Wuyts 2020).Some of them have multi-epoch observations by SDSS, including those from the Baryon Oscillation Spectroscopic Survey (BOSS, Eisenstein et al. 2011) and the Extended Baryon Oscillation Spectroscopic Survey (eBOSS, Dawson et al. 2016).For the convenience of later comparison with millimeter observations, we use optical properties derived from the spectra taken at the closest time to our NOEMA observations.The physical properties of the six quasars are listed in Table 1.The virial black-hole mass listed in Table 1 was based on the broad CIV lines (Shen et al. 2011), with a typical of M BH ≳ 10 9 M ⊙ (except DW004), placing them among the most massive quasars.The bolometric luminosities (L bol ) of the selected sources are 10 46.4−47.4erg s −1 .Figure 1 shows their positions in the redshift-L bol plane.Our sources have comparable bolometric luminosities with AGNs at similar redshifts from the literature.

Observations and data reduction
We observed our targets with NOEMA (S20BT, PI: Dai) in the 2 mm band with 10 antennas on June 6, 15, 19, and September 17, 2020 (Table 2).The compact D configuration was chosen to achieve the highest sensitivity.
We used the Herschel coordinates of the sources as the phase centers (Figure 2).The targets were observed with the PolyFix correlator with two sidebands of 7.744GHz bandwidths.At z = 2 -3, the equivalent velocity coverage is ∼14000-16000 km s −1 in each sideband.We adjust the ∼140-160 GHz spectral windows (Table 2) to cover the 12 CO(4-3) (rest frequency 461.040GHz) for z < 2.5 targets, and 12 CO(5-4) (rest frequency 576.267GHz) for z > 2.5 targets.The expected CO lines are set close to the center of one 3.8 GHz baseband.The native channel width was 2 MHz, and resampled to ∼20 MHz during the data calibration process, corresponding to ∼40 km/s.We took advantage of the track sharing mode, and grouped our sources into two frequency tuning set (DW001, DW002,DW003,and DW006 in S20BT001 and DW004 and DW005 in S20BT002).The observations on June 15 of S20BT001 was not used due to bad data quality.For S20BT001 and S20BT002, the phase and amplitude calibrators were J0122-003 and 0215+015, respectively; the radio frequency (RF) calibrators were 3C454.3 and 3C84, respectively; and the flux calibrators were 1749+096, MWC349, 2010+723, 0106+013, 0215+015.
We used the GILDAS CLIC and MAPPING packages to manually calibrate and reduce the data.The dust continuum was extracted from the uv tables and the flux was derived with 2-D Gaussian fit in the uv plane.A 5σ upper limit was assigned for non-detections, where σ is the 1σ dispersion of the brightness distribution within a 30 ′′ × 30 ′′ field centered at the source.The CO fluxes were extracted with 2-D Gaussian fits on the uv plane for each channel, generated from the continuum-subtracted uv tables.The CO spectra were then constructed by combining the extracted fluxes of every channel.DW003 and DW005 show two CO-emitting components exceeding 4σ on the velocity-collapsed images, so we used two Gaussian components to extract the fluxes individually onn the uv plane.On the extracted spectra, we fit Gaussian profiles (±1000 km s −1 around the channel of the peak flux) and derived the line properties such as the FWHM, the peak flux density, and the integrated flux.
To test if the sources are resolved, we also extracted the continuum fluxes with a point source model for comparison.For DW001 and DW005, the point-source extracted fluxes are consistent with the 2D-Gaussian extracted fluxes (within 1σ) and the FWHMs of the Gaussian profiles are smaller than the beam size.Therefore, we treat DW001 and DW005 as unresolved and use the point-source extracted fluxes.For DW002, DW003b and DW006, the 2D-Gaussian model gives a larger flux and cleaner residual.The typical fitted FWHM is ∼ 4 ′′ , which is marginally resolved with our beam sizes of ∼ 3 ′′ − 4 ′′ , thus fluxes measured with the Gaussian profile were used.
Figure 2 displays the dust continuum emission together with the CO emission images that were obtained by collapsing the datacubes within the fitted velocity ranges of the CO emission lines.We adopted natural weighting for the mapping process, set the cleaning threshold to be 50% of the 1σ noise of the 30 ′′ × 30 ′′ dirty map, and tried various degrees of tapering to all sources.For the CO emission of DW003, tapering creates larger synthesized beams (6.2 ′′ × 3.7 ′′ ), thus concentrating the smeared-out flux and providing a higher peak-flux signal-to-noise ratio (SNR), which just exceeds our detection criterion of 5σ (from SN R ∼ 4 to ∼ 6).For other sources where the detection is not affected by tapering, we kept the original resolution (typical beam size is 3.5 ′′ × 2.5 ′′ , see Table 4).In both maps, the solid contours represent positive values from 0, while the dashed contours are negative values from -1σ, with a 1σ spacing (see Table 3 for the σ values).The beam sizes are plot in the bottom left corners.The quasar's SDSS positions are marked with black crosses, and the white stars in the CO maps mark the peak positions of the continuum.Note that for DW003, the peak of the dust continuum (DW003b) is significantly offset from the quasar's position (DW003a), while the CO emission comes mostly from the optical quasar (DW003a), though with a 4.3 ′′ offset.The right column shows the continuum-subtracted CO spectra, with a velocity resolution of ∼80 km s −1 .Red curves are the Gaussian fits for the lines.The vertical red dotted lines mark the fitted CO line centers, with their 1σ errors in pink shadows (Table 3).The vertical blue dotted lines and shades mark the expected CO frequency and 1σ dispersion based on the optical quasars' redshifts.For DW005a, the purple lines mark the positions for the two peaks, and the red vertical line represent the flux-weighted central position (see Table 3).
3. DUST CONTINUUM, MOLECULAR GAS PROPERTY AND MORPHOLOGY

Dust Continuum
We detect the dust continuum (peak flux density > 5σ) in four of the six sources: DW001, DW002, DW005, and DW006.We also detect a > 5σ dust continuum ∼90 kpc away from the DW003 optical position.The left column of Figure 2 shows the continuum contours overlaid on the optical images from SDSS.The measured dust continuum flux densities are listed in Table 3.
We note that the dust continuum of DW003 at the optical position (i.e,DW003a) is only of 3.5σ.However, ∼ 10 ′′ north-east, there is a > 6σ continuum component (DW003b) without an optical counterpart, suggesting an optically-obscured submillimeter galaxy (SMG).Given the tentative CO detection at the same redshift, DW003b is likely associated with the targeted quasar (Section 3.2).We note that the 250µm HerS position is between DW003a and DW003b, indicating blending in the FIR dust emission.For DW004, we detect nothing at 3σ-level.For DW005, although the continuum is detected at 7σ, ∼ 4.7 ′′ east of the quasar, we find a companion with significant CO emission only (DW005b, see Section 3.2 and Figure 2).To calculate the continuum flux upper limits of DW004 and DW005b, we use the 5σ upper limit (see Section 2.2) times the synthesized beam size.The slightly extended 3σ continuum contour of DW006 could be a possible lensing effect, which will be discussed in Section 4.1.We also find spatial offsets and velocity shifts between the optical and mm observations, which will be described and discussed in Section 3.3 and Section 4.3.

Molecular CO Emission
The central column of Figure 2 shows the molecular CO emissions, and the right column is the continuum-subtracted, integrated spectra.The CO(4-3) or CO(5-4) emission are detected with a peak SNR level of ∼ 6 − 13 for five of the six sources, except DW004.The measured line properties are listed in Table 3.
For DW001 and DW002, the CO(5-4) emissions both show a slightly extended structure along the major axis of the clean beam.Both have a CO emission that aligns well with the 2mm dust continuum.
For DW003, there was no >5σ detection within the original beam resolution (4.02 ′′ ×2.78 ′′ ) so tapering was adopted.After tapering, a larger beam size (6.16 ′′ ×3.73 ′′ ) yields a 6σ CO(4-3) detection (DW003a), which extends to the south-east of the optical quasar position (cross in Figure 2).This corresponds to a molecular gas reservoir up to 10 ′′ (SNR>3 region of ∼80 kpc).To the north-west of DW003a, at the location of DW003b, we also detect a weak CO(4-3) emission of ∼ 3σ.The ratio between continuum flux density and CO peak flux density of DW003b (∼ 0.3) is higher than other sources (< 0.1, Table 3), which may indicate the existence of an obscured AGN.DW003a and DW003b have almost identical emission line frequencies, suggesting that they are in a pair system.The projected distance is ∼ 110 kpc (∼ 14 ′′ ).Similar separations have been reported before (e.g., between NGC7679 and NGC7682 by Ricci et al. 2021), suggesting possible interaction between the two systems, which might have result in the extended CO morphology of DW003a.
DW004 is not detected in CO, despite integrating through a range of 2000 km s −1 around the expected frequency.The spectrum of DW004 is subtracted using a polygonal aperture covering the phase center, which is of similar sizes to the other detected sources (∼10 beams).We estimate the upper limit of CO(4-3) emission flux as 5σ noise times one beam size.
DW005 have two components.The higher SNR CO(4-3) feature aligns well with the peak of the 2mm dust continuum (DW005a).The related CO emission shows a double-peak line profile, with separation of ∼ 650 km s −1 .The doublepeaked profile may indicate either disk-like rotation in the system, or two distinct velocity components.In the following analysis, we fit the two lines separately, and also calculate an average velocity for DW005a, weighted by the fitted fluxes of the two lines.The average velocity corresponds to a redshift is 2.283±0.001.We also find a potential second CO emitter (DW005b) ∼ 5 ′′ east of DW005a (Figure 2).DW005b has an integrated CO flux of > 3σ, at an almost identical redshift of the bluer peak of DW005a, though not associated with any > 3σ dust continuum(Section 3.1).The different line profiles of the two components suggest that they are not likely caused by gravitational lensing.
In DW006, the CO(5-4) line is detected at a ∼ 10σ level with a relatively narrow width (FWHM = 250 km s −1 ).Coincident dust continuum and optical positions, the strong emission and narrow FWHM are indicative of a lensed system as will be discussed in Section 4.1.

Positional and Velocity Offsets
We find common (3 out of 5 detected) positional offsets and velocity shifts between the NOEMA continuum emission and SDSS optical observations in our sample, as listed in Table 4.However, the moment maps, generated using the GILDAS MAPPING package, do not show any significant signs of velocity components.
For DW001, DW002, DW005 and DW006, the peak positions of the 2mm continuum are offset from the quasar optical positions by ∼ 0.4 ′′ − 1.6 ′′ .We calculate the spatial uncertainties using where δθ is the positional uncertainty, θ Beam is the cleaned beam size (3 ′′ − 4 ′′ ), and SN R peak is the signal-to-noise ratio of the peak detection on the map (Reid et al. 1988).All our sources have consistent Gaia (Gaia Collaboration et Liu et al. (2) the size and (3) positional angle (PA) of the clean beam; (4) the angular distance between the CO emission peak and the 2mm continuum peak; (5) the angular distance between the 2mm continuum peak and the optical quasar position; (6) the physical projected distance between the 2mm continuum peak and the optical quasar position; (7) the velocity difference between the CO emission line and optical redshifts with both optical and CO redshift uncertainties considered (positive values indicate a redshifted CO emission compared to the optical lines).Since DW003a and DW005b do not have SNR>5 continuum detection, we only list the separation between the CO peak and the optical positions in column ( 5) and ( 6).For DW003b, the spatial uncertainty is a combination of NOEMA pointing uncertainty and positional uncertainty.
For others, a 0.5 ′′ spatial uncertainty is the pixel scale.For DW005a, the velocity offset is the average between the two peak positions, weighted by their relative fitted fluxes.
al. 2016, 2021) positions with the SDSS coordinates, except for DW003 and DW005, whose Gaia positions are 40 mas north.Since 40 mas is negligible compared to the 2mm uncertainties, in the following analysis and Table 4, all offsets are calculated based on the SDSS optical positions.The positional uncertainty δθ, the NOEMA pointing accuracy of 0.2 ′′ , the NOEMA pixel size of 0.5 ′′ , and the SDSS positional uncertainty of < 0.1 ′′ together yield a spatial uncertainty of 1.0 ′′ for DW003b, and 0.5 ′′ for the other sources.Since DW003a and DW005b have no significant dust continuum detection, in Table 4 (column 5 and 6) we present their CO flux peak position offsets relative to the optical positions.
Column (2) in Table 4 lists the offsets of the peaks between the CO and dust continuum emissions (s CO−mm ) with an error calculated from CO and continuum spatial detection uncertainties, NOEMA pointing accuracy and the pixel size.For DW001 and DW002, their offsets are consistent within the astrometric accuracy of NOEMA (0.2 ′′ ).DW005a has a 1.0 ′′ positional offset, but of different direction compared to DW005b (south vs. south-east), thus likely not due to the latter.Misaligned gas and dust components suggest that DW005a itself may be in a pair system, while DW005b is a third component outside the pair.Difference in the spatial distribution of the molecular gas and dust continuum has been observed with small separations at higher redshifts (≤ 1 ′′ , e.g.Gururajan et al. 2021;Fogasy et al. 2022;Lamperti et al. 2022), though often times larger offsets are observed between the optical and submm components for both AGN/quasars and star forming galaxies (see Section 4.3).
In general, we find spatial offsets of ∼ 1 − 4 ′′ between the Herschel and the optical quasar positions.We note that the positional errors are at 6 ′′ level, which propagate from the pointing accuracy (2 ′′ , Pilbratt et al. 2010), spatial detection uncertainty (∼ 2 ′′ ), and the pixel size (6 ′′ at 250µm, Viero et al. 2014).Given the large positional errors, the optical-FIR offsets are not discussed later.
At z ∼ 2 − 3, 1 ′′ corresponds to ∼ 8 kpc.Thus for our targets, the observed projected distance correspond to 3 to 13 kpc.This offset is significant compared to the typical galaxy sizes at the Cosmic Noon (a few kpc, Förster Schreiber & Wuyts 2020).For targets with a companion 2mm continuum component (DW003 and DW005), their offsets to the quasar positions range from 35 to 90 kpc, suggesting that the second mm source is likely another galaxy.
Velocity shifts are also observed between the optical and mm spectroscopic redshifts (Table 4).The red dashed lines in Figure 2 (right) mark the expected frequencies and ranges based on the SDSS spectroscopic redshifts and associated uncertainties.The differences between the optical and CO redshifts (i.e.δv) are all positive, corresponding to redshifted CO lines relative to the optical lines.In DW001 and DW002, the velocity difference is relatively small with a large uncertainty and can be treated as consistent, despite the relatively large spatial offset in DW002 (Column ( 5) and ( 6) in Table 4).The velocity difference is large (1100-1800 km s −1 ) for DW003, DW005, and DW006.We find significant velocity shifts, as compared to the redshift uncertainties, in DW003a and DW005a.These two systems may undergo volatile kinematic activities, possibly related to the secondary components.The origin and nature of these offsets will be further discussed in Section 4.3.
To estimate the molecular gas mass, we adopt a linear relation between L ′ CO and H 2 masses (i.e., M H2 = αL ′ CO(1−0) ), assuming α = 0.8 M ⊙ (K km/s pc 2 ) −1 (Dunne et al. 2022).For all detected targets (Table 5), this yields M H2 of the order of 10 10 M ⊙ .Note that these values are not corrected for possible magnification.The uncertainties listed in the table are based on the line measurement, without considering systematic uncertainties on the conversion factors.For instance, if we adopted an α = 4.0 M ⊙ (K km/s pc 2 ) −1 , the resulting M H2 would have be 5 times larger.In addition, we note that the conversion factors between different CO lines from Carilli & Walter (2013) are average values, with a large scatter (e.g., a scatter of 0.5 dex have been reported in CO(4-3) to CO(1-0), see e.g.Banerji et al. 2018;Bothwell et al. 2014;Brusa et al. 2018).

SED fitting
The SEDs for the six sources include fluxes and upper limits from: SDSS (Alam et al. 2015;Ahumada et al. 2020), the UKIRT Infrared Deep Sky Survey (UKIDSS, Lawrence et al. 2007), the Wide Infrared Survey Explorer (WISE, Wright et al. 2010), and HerS (Viero et al. 2014, see Table 1).The measured 2mm flux density greater than 5σ is used, otherwise a 5σ upper limit (Section 3.1) is adopted.We fit the SEDs of our HyLIRG-quasars using the CIGALE code (Yang et al. 2022).We adopt an AGN component and a cold dust component in our SED fitting, with the CO redshifts as input.Table 6 lists the configuration for the fits.
We note that the FIR fluxes of DW003a in the HerS catalog may suffer from blending of both DW003a and DW003b (See Section 3.1).However, the current Herschel resolution (∼ 18 ′′ at 250µm, Viero et al. 2014) makes it difficult to assign the flux.Here we assign 31% of the Herschel fluxes to DW003a, which corresponds to DW003a and b's flux ratio in the 2mm.At 500µm, the beam size of SPIRE is ∼ 36 ′′ (Viero et al. 2014).Considering possible contamination from nearby sources, we manually increase the error budget for 500µm to 50%in the fit.
The properties derived from SED fitting are listed in Table 7.For DW001, DW002 and DW006, the µL IR,SB we calculate is consistent with the values in Dong & Wu (2016) (Table 1) within ∼ 1σ.Their µL IR,SB still satisfy the HyLIRG definition.However, if a magnification factor of 5-10 is applied for DW006, as derived from Figure 4 in Section 4.1, its intrinsic L IR,SB will drop to ULIRG level.For DW003a, in contrast to the 10 13 L ⊙ cold-dust IR luminosity reported by Dong & Wu (2016), the fitted µL IR,SB is also in ULIRG level.In addition, the dust mass is of 10 8 M ⊙ , which is only ∼ 1% of molecular gas mass (Table 5).This suggests a dusty companion source that contributes to the majority of the 2mm emission.The µL IR,SB of DW004 is only 0.6 × 10 13 L ⊙ , about 50% lower than in Dong & Wu (2016).For all sources except DW003a, the fitted dust mass is of the order of 10 9 M ⊙ , indicating the dust-rich nature of these HyLIRG-quasars.
In Figure 4, galaxies clearly fall in two distinct populations: strongly lensed galaxies in the upper left and unlicensed along the power-law relation.The unlensed or weakly lensed galaxies in general follow a virial relation: where ∆V is the FWHM of the CO line in km s −1 , R is the radius of the CO emission region in parsecs, α is the conversion factor from L ′ CO(1−0) to solar mass in K km s −1 pc 2 , G is the gravitational constant, and C is a constant related to the kinematics of the galaxy.We consider two extreme cases, using parameters suggested by Erb et al. (2006): C = 2.1, R = 5 kpc, α = 4.6 for a disk model; C = 5, R = 2 kpc, α = 1.0 for a spherical model.Both models are plotted with dotted lines in Figure 4, which nicely bracket the majority of the unlensed and slightly-lensed galaxies.
The solid line in Figure 4 represents the best-fit relation L ′ CO = 5.4 × ∆V 2 derived by Bothwell et al. (2013) and also applied by (e.g., Zavala et al. 2015;Neri et al. 2020).The dashed line represents the best-fit relation from Harris et al. (2012), i.e., L ′ CO(1−0) (K km s −1 pc 2 ) = (∆V/400km s −1 ) 1.7 × 10 11 /3.5.In our sample, four out of the five CO-detected sources are safely located in the unlensed or at most weakly lensed region.For DW005a, whose CO emission shows double peaks, the separation between the two peaks was used as it's line width.Despite the large error bars, the two possible companions, DW003b and DW005b, also fall in the unlensed region.
For objects with multiple CO lines observed (e.g.Bothwell et al. 2013), we calculate the low-J (filled circle) and mid-J (triangle) emissions separately, and mark both in the diagram, using L ′ CO /L ′ CO(1−0) ratios from Carilli & Walter (2013) for quasars and SMGs.Detections with SNR≤3 (tentative detections) are not plotted.
Based on the L ′ CO -∆V diagram, DW001, DW002, DW003a, DW005a and two possible companions (DW003b, DW005b) fall in the unlensed power-law region, thus likely unlensed or at most weakly lensed, intrinsic HyLIRGs.DW004 is undetected thus not plotted.DW006 is likely strongly lensed, with an estimated magnification factor of ∼ 5 − 10, based on the offset from the power-law relation, making it a ULIRG system instead.vs the CO FWHM (∆ V ) for our sample of 6 HyLIRG-quasars, and 180 galaxies from the literature, including lensed and unlensed SMGs and DSFGs (Bothwell et al. 2013;Aravena et al. 2016;Yang et al. 2017;Harris et al. 2012;Bakx et al. 2020), local/low-to-mid-z ULIRGs (Solomon et al. 1997;Combes et al. 2011Combes et al. , 2013Combes et al. , 2006)), and high-z quasars and companions (Fan et al. 2018(Fan et al. , 2019;;Riechers 2011;Riechers et al. 2011;Wang et al. 2010;Feruglio et al. 2017;Bischetti et al. 2021;Noterdaeme et al. 2021).Data points using CO high J transitions (J≥3) are marked with triangles and converted to CO(1-0) using the same factors as in this paper (i.e., Carilli & Walter 2013), while data points using CO J≤2 transitions are marked with filled circles.The solid and dashed lines represents the approximate best-fitting quadratic relationships from Bothwell et al. (2013) and Harris et al. (2012), respectively.The dotted lines represent the virial relations assuming spherical and disk models.

Depletion Time of the HyLIRGs
Based on the gas mass and SFR calculated in Section 3.4 and Section 3.5, we estimate the gas depletion time: The calculated τ dep (Table 7) ranges between 20 − 60 Myr, similar to other starburst galaxies (tens of Myrs, e.g.Daddi et al. 2010;Combes et al. 2013), for all of our sources except DW003a.This τ dep range is much shorter than the lifetime of a galaxy (∼ 10 Gyr for elliptical galaxies, De Lucia et al. 2006).This is consistent with the scenario that HyLIRG-quasars are in a short transitional phase in the early formation of the massive elliptical galaxies (Fu et al. 2013).We also note that if a Salpeter IMF is applied, the SFR would also increase by ∼ 0.15 dex (Davé 2008), resulting Liu et al.
in even shorter τ dep of ∼15-40 Myr.We also note that, as mentioned in Section 3.4, the depletion time could be up to 5 times longer if instead of 0.8, a CO-H 2 conversion factor α of 4.0 M ⊙ (K km s −1 pc 2 ) −1 is adopted.
In Figure 5, we compare the inverse integrated Kennicutt-Schmidt relation between the molecular gas mass and µSFR of our sources, and observations from the literature.Using the µM H2 as a proxy for total molecular gas mass, we find that the star formation efficiencies (SFEs) of all CO-detected sources follow the trend for starburst galaxies, similar to AGNs from literature.6 dex to represents the starburst galaxies.The pink strip mark the starburst regions, whose lower border represent the extreme starburst galaxies with a SFE ∼ 15 times higher than that of the MS galaxies (Sargent et al. 2014).Our five CO-detected quasars all fall in the starburst region, with a CO-H2 conversion factor α = 0.8M⊙(K km s −1 pc 2 ) −1 (red stars), similar to other starburst galaxies (Tacconi et al. 2020).They will fall in a transitional region between starburst and MS galaxies if an α = 4.0M⊙(K km s −1 pc 2 ) −1 is adopted (white stars).Grey dots are the star-forming galaxies in the COSMOS deep field from Liu et al. (2019) at z=2-3 plotted for comparison.Filled circles and triangles (upper limits) are AGNs from the literature (Perna et al. 2018;Kirkpatrick et al. 2019;Bischetti et al. 2021;Circosta et al. 2021).No correction of gravitational magnification is applied to the sources.

Origin of the Observed Positional and Velocity Offsets
In Section 3.3, we present the spatial offset and velocity shifts between the optical quasars and the mm dust and molecular gas components.In this section, we investigate the possible causes of these shifts.Similar offsets have been ubiquitously observed in quasars and star forming galaxies, where the location of the molecular CO or dust components deviate from either the galaxies' optical positions, or redshifts, or both (e.g.Krips et al. 2005;Clements et al. 2009;Combes et al. 2013;Iono et al. 2016;Chiaberge et al. 2017;Ikarashi et al. 2017;Magdis et al. 2017;Barthel et al. 2018;Vayner et al. 2021).
One likely explanation for the systematic redshifted CO velocity is the known blueshift in the broad C IV emission lines in the quasar system, whose peak is often blue-shifted with respect to narrow-line emissions, such as [O III], used as the representative systematic redshift of the quasar system when available (see, e.g., Richards et al. 2011;Coatman et al. 2017;Vietri et al. 2020), and to cold gas tracers such as [C II] and CO (see, e.g., Bischetti et al. 2017;Trakhtenbrot et al. 2017;Schindler et al. 2020;Circosta et al. 2021).For our sources at z∼2-3, indeed the redshift determination from SDSS is highly reliable on the broad C IV broad emission lines, partly due to the absence of [O III] or [Ne V] lines in the wavelength coverage.We double checked the spectra of DW002, DW003 and DW005, whose Mg II broad emission lines are also well detected.We find that all of the Mg II lines are ≲ 1000 km s −1 redshifted from the SDSS redshifts, fixed as the Gaussian center wavelength when fitting the Mg II line profiles.If we treat the Mg II redshift as the quasar's representative redshift, and remove any velocity difference from the CO-optical velocity offsets, which now become -700, 600, and 500 km s −1 for DW002, DW003a and DW005, respectively.We note that the intrinsic blueshift of the broad-line-region could be as large as ∼2000 km s −1 (e.g.Schindler et al. 2020;Vietri et al. 2018), which could explain the velocity shifts observed in our quasar sample.A comparison of optical and CO redshifts on broad-line spectra is demonstrated in Figure 6.
Another scenario involves unresolved mergers in the HyLIRG-quasar systems, as in the cases in Krips et al. (2005) and a number of similar targets (i.e. a powerful AGN with a CO-rich merging partner, e.g.Walter et al. 2004;De Breuck et al. 2005).Merging galaxies with a dust-free quasar and an optically-obscured SMG could explain the observed offsets, both positionally and spectroscopically, with the SMGs contributing the majority of the observed CO and dust emissions.If two galaxies are in the later stage of a merger, they cannot be resolved under the current resolution (∼ 3.5 ′′ × 2.5 ′′ ), even if the CO emission is centered at a second component slightly offset from the quasar location.The SMG + quasar scenario could also explain the observed velocity offsets up to several hundred of km s −1 .The merging of two gas-rich galaxies can trigger starburst and explain the observed HyLIRG level luminosities, as has been demonstrated in the cases of ULIRGs (e.g.Riechers 2011;Ivison et al. 2019;Huang et al. 2023).In fact, DW003a has extended CO morphology, implying a second gas clump or galaxy at almost identical redshifts at a distance of ∼90 kpc from the main component.The double-peaked CO line profile of DW005a indicates disk-like rotation, or two distinct velocity components.This offers indirect evidence to support the merger scenario.
Finally, we can not rule out the recoiling black holes scenario: the central black hole being ejected with the broad line region (BLR) during galaxy merging, while the narrow lines still lag behind in the center of the galaxy, as well as the molecular gas (Komossa 2012).Velocity offsets generated by this effect can be of the order of 100-1000 km s −1 (as discussed in Chiaberge et al. 2017), which is consistent with the observed velocity offset in our sample.However, simulation suggests the probability of such an event is low (only < 10% of recoiling black holes formed from binary black holes have velocity > 1000 km s −1 , Civano et al. 2010).1).For DW001 and DW002, the optical and CO redshifts are consistent and overlaid with each other.For DW003, DW005, DW006, the CO lines are redshifted from the peaks of optical lines.

Figure 2 .
Figure2.The 2mm continuum (left) and CO (middle) maps, and CO spectra (right) of DW001 to DW006 (from the top).In both maps, the solid contours represent positive values from 0, while the dashed contours are negative values from -1σ, with a 1σ spacing (see Table3for the σ values).The beam sizes are plot in the bottom left corners.The quasar's SDSS positions are marked with black crosses, and the white stars in the CO maps mark the peak positions of the continuum.Note that for DW003, the peak of the dust continuum (DW003b) is significantly offset from the quasar's position (DW003a), while the CO emission comes mostly from the optical quasar (DW003a), though with a 4.3 ′′ offset.The right column shows the continuum-subtracted CO spectra, with a velocity resolution of ∼80 km s −1 .Red curves are the Gaussian fits for the lines.The vertical red dotted lines mark the fitted CO line centers, with their 1σ errors in pink shadows (Table3).The vertical blue dotted lines and shades mark the expected CO frequency and 1σ dispersion based on the optical quasars' redshifts.For DW005a, the purple lines mark the positions for the two peaks, and the red vertical line represent the flux-weighted central position (see Table3).

Figure 3 .
Figure 3. SEDs of our HyLIRG-quasar sample.The FIR dust emission is plotted in red, while the AGN component is in orange and in solid black is the composite SED.Violet empty circles mark the observed photometric data points, with upper limits shown in downward triangles.Besides the NOEMA 2mm continuum fluxes, the data points are from SDSS, UKIDSS, WISE, Herschel, and ALMA (see Section 3.5).

Figure 5 .
Figure 5.The inverse integrated Kennicutt-Schmidt relation.The blue strip represents the relation for main-sequence (MS) galaxies, with a dispersion of 0.3 dex.The solid black line is the MS relation consistent with Figure 3 in Tacconi et al. (2020) at z ∼ 2.5.The dashed orange line represents the MS relation offset by 0.6 dex to represents the starburst galaxies.The pink strip mark the starburst regions, whose lower border represent the extreme starburst galaxies with a SFE ∼ 15 times higher than that of the MS galaxies(Sargent et al. 2014).Our five CO-detected quasars all fall in the starburst region, with a CO-H2 conversion factor α = 0.8M⊙(K km s −1 pc 2 ) −1 (red stars), similar to other starburst galaxies(Tacconi et al. 2020).They will fall in a transitional region between starburst and MS galaxies if an α = 4.0M⊙(K km s −1 pc 2 ) −1 is adopted (white stars).Grey dots are the star-forming galaxies in the COSMOS deep field fromLiu et al. (2019) at z=2-3 plotted for comparison.Filled circles and triangles (upper limits) are AGNs from the literature(Perna et al. 2018;Kirkpatrick et al. 2019;Bischetti et al. 2021;Circosta et al. 2021).No correction of gravitational magnification is applied to the sources.

Figure 6 .
Figure6.Broad-emission-line spectra of DW001-DW006.The red dashed vertical lines mark the CO redshifts, and the blue dashed vertical lines mark the optical redshifts from SDSS (Table1).For DW001 and DW002, the optical and CO redshifts are consistent and overlaid with each other.For DW003, DW005, DW006, the CO lines are redshifted from the peaks of optical lines.

Table 1 .
Source properties Source name Herschel name b for our SFR estimates. 1 Figure 1.Redshifts of selected sources versus their bolometric luminosities.Red filled squares are the six HyLIRG-quasars with bolometric luminosities from

Table 2 .
Observations , the frequency coverage is 135.692GHz to 143.436GHz (the lower sideband) and 151.180GHz to 158.924GHz (the upper sideband).For DW004 and DW005, the frequency coverage is 144.576GHz to 152.320GHz (the lower sideband) and 160.064GHz to 167.808GHz (the upper sideband).b Precipitable water vapor.

Table 3 .
Observed CO and dust properties (GHz) (km s −1 ) (Jy km s −1 ) (mJy) (Jy/beam km s −1 )a With a typical error of ∼0.04-0.06GHz.For those that have fitted uncertainty of CO line frequency lower than 0.04 GHz, we

Table 5 .
Molecular gas properties

Table 7 .
SED fitting resultsNote-fIR,AGN is fraction of AGN emission in total IR luminosity estimated from SED fitting with the CIGALE code.µLIR,SB is the AGN-removed, starburst-dominant IR luminosity.µM dust is the fitted galactic dust mass.The luminosities, SFRs and masses are apparent quantities not corrected for possible gravitational lensing effect.τ dep is the depletion time calculated from Equation 3.