DDraft version January 5, 2021
Typeset using L A TEX twocolumn style in AASTeX63
Compact Molecular Gas Distribution in Quasar Host Galaxies
Juan Molina, Ran Wang,
1, 2
Jinyi Shangguan, Luis C. Ho,
1, 2
Franz E. Bauer,
4, 5, 6
Ezequiel Treister, andYali Shao Kavli Institute for Astronomy and Astrophysics, Peking University, Beijing 100871, China Department of Astronomy, School of Physics, Peking University, Beijing 100871, China Max-Planck-Institut f¨ur Extraterrestrische Physik (MPE), Giessenbachstr., D-85748 Garching, Germany Instituto de Astrof´ısica and Centro de Astroingenier´ıa, Facultad de F´ısica, Pontificia Universidad Cat´olica de Chile, Casilla 306,Santiago 22, Chile Millennium Institute of Astrophysics (MAS), Nuncio Monse˜nor S´otero Sanz 100, Providencia, Santiago, Chile Space Science Institute, 4750 Walnut Street, Suite 205, Boulder, Colorado 80301 Max-Planck-Institut f¨ur Radioastronomie, Auf dem H¨ugel 69, 53121 Bonn, Germany (Accepted ApJ)
ABSTRACTWe use Atacama Large Millimeter/submillimeter Array CO(2–1) observations of six low-redshiftPalomar-Green quasars to study the distribution and kinematics of the molecular gas of their hostgalaxies at kpc-scale resolution. While the molecular gas content, molecular gas fraction, and starformation rates are similar to those of nearby massive, star-forming galaxies, the quasar host galaxiespossess exceptionally compact, disky molecular gas distributions with a median half-light radius of1.8 kpc and molecular gas mass surface densities (cid:38) M (cid:12) pc − . While the overall velocity field of themolecular gas is dominated by regular rotation out to large radii, with rotation velocity-to-velocitydispersion ratio (cid:38)
9, the nuclear region displays substantial kinematic complexity associated withsmall-scale substructure in the gas distribution. A tilted-ring analysis reveals that the kinematicand photometric position angles are misaligned on average by ∼ ± ◦ , and provides evidence ofkinematic twisting. These observations provide tantalizing clues to the detailed physical conditions ofthe circumnuclear environments of actively accreting supermassive black holes. Keywords: galaxies: active — galaxies: kinematics and dynamics — quasars: general INTRODUCTIONCorrelations between host galaxy bulge properties andsupermassive black hole (BH) mass (Kormendy & Rich-stone 1995; Ferrarese & Merritt 2000; Gebhardt et al.2000; Tremaine et al. 2002) have led to the conclusionthat galaxies and BHs coevolve (Richstone et al. 1998;Kormendy & Ho 2013), perhaps mediated by feedbackmechanisms from active galactic nuclei (AGNs) that reg-ulate gas accretion toward the center and star formationon galactic scales (e.g., Fabian 2012; Heckman & Best2014). AGN feedback can produce strong multi-phasegas outflows (e.g., Cicone et al. 2014; Perna et al. 2015;Feruglio et al. 2015; Karouzos et al. 2016; Morganti 2017;
Corresponding author: Juan [email protected]
Cicone et al. 2018; Fluetsch et al. 2019) that are powerfulenough to heat or remove the interstellar medium fromthe host galaxy (Silk & Rees 1998; Harrison et al. 2018),thereby suppressing subsequent star formation activity(e.g., Dubois et al. 2016) and keeping the host galaxyquiescent (Fabian 2012). It is a key ingredient in numer-ical, theoretical and semi-analytic models to reproducethe lack of massive galaxies in the high-mass end of themass function (e.g. Kauffmann & Haehnelt 2000; Cro-ton et al. 2006; Schaye et al. 2015; Sijacki et al. 2015;Lacey et al. 2016).Notwithstanding these compelling arguments, seriousdoubts remain as to whether AGN feedback effectivelyremoves sufficient cold gas from the host galaxy to cur-tail its ongoing star formation activity. The interstellarmedium content of the host galaxies of nearby AGNsshows no evidence of depletion relative to star-forminggalaxies of similar stellar mass, based on gas masses de- a r X i v : . [ a s t r o - ph . GA ] J a n Molina, J. et al. rived from direct observations of neutral atomic hydro-gen (Ho et al. 2008; Fabello et al. 2011; Ger´eb et al. 2015;Zhu & Wu 2015; Ellison et al. 2019) and CO (Maiolinoet al. 1997; Evans et al. 2001; Scoville et al. 2003; Evanset al. 2006; Bertram et al. 2007; Shangguan et al. 2020a;Jarvis et al. 2020), as well as inferred indirectly fromdust emission (Shangguan et al. 2018; Shangguan & Ho2019) or dust extinction (Zhuang & Ho 2020; Yesuf &Ho 2020). Far from being quenched, stars seem to formwith even greater efficiency in the host galaxies of lu-minous AGNs (Shangguan et al. 2020b; Zhuang & Ho2020). It should be further noted that AGN host galax-ies possess not only a “normal” gas reservoir, but thegas appears largely kinematically regular, as evidencedby their global line widths in H I (Ho et al. 2008) andCO (Shangguan et al. 2020b), which are consistent withrotational support, and by the absence of significantmolecular outflows (Shangguan et al. 2020a).The above-mentioned observations suggest that AGNfeedback, if present, imparts only a modest, likely lo-calized effect on the cold gas. For example, the exist-ing molecular gas observations cannot rule out that thecold gas is in the process of being expelled from nucleibut still associated with the host galaxy. And whileAGNs can heat the molecular gas and suppress starformation through “negative” feedback (Papadopouloset al. 2010), they can also exert the opposite effect—“positive” feedback—that can compress the cold gasand enhance star formation (Cresci et al. 2015; Carnianiet al. 2016; Maiolino et al. 2017; Cresci & Maiolino 2018;Gallagher et al. 2019). From an observational pointof view, the impact of AGN feedback is still far fromsettled, and it is essential to obtain spatially resolvedinformation on the molecular gas in active galaxies togain further insight into the physical processes that gov-ern the coevolution of supermassive BHs and their hostgalaxies.Quasars, the most luminous of the active galaxies, arethe ideal sites to investigate the possible interplay be-tween AGN feedback and the molecular gas of their hostgalaxies. In the popular merger-driven evolutionary sce-nario of quasars (Sanders et al. 1988), two gas-rich galax-ies merge, gravitational torques drive the cold gas to thecenter of the merger remnant, and vigorous starburst ac-tivity and BH growth ensue. The prodigious energy re-leased by the AGN expels the enshrouding gas and dust,giving birth to an optically visible, largely unobscuredquasar (Hopkins et al. 2008).The frequent association of a young stellar populationwith quasar host galaxies supports the notion that starformation accompanies or precedes AGN activity (e.g.,Canalizo & Stockton 2001, 2013; Jahnke et al. 2007; Kim log ( λL λ (5100 Å ) / erg s − ) l og ( M B H / M fl ) This workShangguan+20aCO − Literature z < . sample Figure 1.
Black hole mass as a function of the AGNmonochromatic luminosity at 5100 ˚A for the PG quasar hostgalaxies. We show separately the host galaxies with CO datafrom literature (open circles), our ACA observations (filledcircles; Shangguan et al. 2020a) and the six targets observedby ALMA (red-filled circles). The six targets sample the M BH (cid:46) M (cid:12) limit within the PG survey. & Ho 2019; Zhao et al. 2019), in qualitative agreementwith the merger-induced evolution scenario of quasars.However, stellar morphology studies of quasar hostgalaxies suggest that only a fraction of the systemsshow tidal and/or dynamical perturbation signaturesthat may be associated to recent merger activity (e.g.Veilleux et al. 2009). Not all the quasar hosts with en-hanced star formation activity show evidence of dynam-ical perturbations (Shangguan et al. 2020b). The equalmass merger scenario may be applicable to ultralumi-nous infrared galaxies and AGNs hosted in ellipticals.But, in other cases, unequal mass mergers which do notperturb the more massive interacting galaxy may be re-quired.This work reports new, relatively high-resolution(beam ∼ . (cid:48)(cid:48) − (cid:48)(cid:48) , which corresponds to physical scalesof (cid:46) z < .
5) type 1(broad-lined) quasars (Boroson & Green 1992). Thissample has been studied extensively, enjoying a richrepository of multi-wavelength data for the AGN andhost galaxy, including optical spectra (Boroson & Green1992; Ho & Kim 2009), radio properties (Kellermannet al. 1989, 1994), X-ray constraints (Reeves & Turner2000; Bianchi et al. 2009), dust properties for both the olecular Gas in Quasar Hosts CO) J = 2 → ν rest = 230 .
538 GHz; here-inafter CO(2–1)] for a subset of 23 z < . z < . m = 0 . Λ = 0 . H = 67 . − Mpc − (Planck Collaboration et al. 2016). SAMPLE AND OBSERVATIONS2.1.
ALMA Observations and Data Reduction
We analyze Cycle 6 ALMA follow-up observations(program 2018.1.00006.S; PI: F. Bauer) of the six PGquasar host galaxies (Table 1) taken from our z < . (cid:46) . (cid:48)(cid:48) . (cid:48)(cid:48) z PG quasarswith low to moderate BH masses and luminosities (Fig-ure 1) and disk-like stellar morphology (Table 1). TheALMA observations were carried during November 2018to March 2019 in good weather conditions with pre-cipitate water vapor (PWV) (cid:46) . ∼
11 km s − . The observational setup for each sourceis described in Table 2. The ALMA flux calibration un-certainty is (cid:46)
10 % (Fomalont et al. 2014; Bonato et al.2018).We use the Common Astronomy Software Application(
CASA ; McMullin et al. 2007) to reduce the ALMAdata, employing the standard pipeline to calibrate thedata to generate the uv visibilities. To minimize miss-ing flux from possible extended source emission and toobtain more sensitive imaging of the total CO(2–1) emis-sion, we concatenate our Cycle 6 observations with theprevious Cycle 5 ACA observations using the task con-cat . For the ACA data we consider the same spec-tral channel flagging that Shangguan et al. (2020a) em-ployed in their work, while for the 12-m antenna datawe flag the spectral channels near a sky feature foundat ∼ . uvcontsub to subtract the continuum, beforeproceeding with imaging the line emission.The emission-line imaging is performed using tclean with robust weighting ( robust = 0.5) and a channelresolution of 11 km s − . The spatial pixel scale is set tosample the synthesized beam by five pixels. Visual in-spection of the masks constructed by the masking proce-dure auto-multithresh indicates that optimal resultscan be obtained by setting the noise, sidelobe, and low-noise thresholds to 4.0, 1.0, and 1.5, respectively, andthe fractional beam size to 0.4. All other parameters of auto-multithresh were left at their default values, asrecommended for combined ACA and 12-m data. TheRMS values and the sizes of the synthesized beam aregiven in Table 2.2.2.
Emission-line Characterization and MapConstruction
We implement a three-step procedure to characterizethe CO(2–1) emission-line shapes. For each pixel, wefirst average the spectrum by considering a box/squaredregion with size comparable to the synthesized beam(e.g., Swinbank et al. 2012), and we estimate the noisefrom the line-free channels. We fit a Gaussian modelto the spectrum using the least-squares minimizationprocedure implemented in the
Python package lm-fit (Newville et al. 2014). As initial guesses, we as-sume that the line centroid equals to the line peak loca-tion in the spectrum, and that the line width equals to20 km s − . The line width is restricted to a maximumvalue of 500 km s − . We use the Bayesian informationcriterion (BIC; Schwarz 1978), which penalizes by the https://casaguides.nrao.edu/index.php?title=Automasking Guide Molina, J. et al.
Table 1.
Basic Parameters of the SampleObject R.A. Decl. z D L Morph. log M (cid:63) log L IR SFR log M gas log M BH log λL λ (5100˚A)(J2000.0) (J2000.0) (Mpc) ( M (cid:12) ) (erg s − ) ( M (cid:12) yr − ) ( M (cid:12) ) ( M (cid:12) ) (erg s − )(1) (2) (3) (4) (5) (6) (7) (8) (9) (10) (11) (12)PG 0050+124 00:53:34.94 +12:41:36.2 0.061 282.3 Disk 11.12 44.94 +0 . − . +0 . − . +0 . − . +0 . − . −
040 10:14:20.69 − +0 . − . +0 . − . −
041 11:29:16.66 − +0 . − . +0 . − . +0 . − . +0 . − . +0 . − . +0 . − . Note — (1) Source name. (2) Right ascension. (3) Declination. (4) Redshift. (5) Luminosity distance. (6) Host galaxymorphology type based on HST image and taken from (Zhang et al. 2016; Kim et al. 2017; Zhao et al. in prep.) (7) Stellarmass; the 1 σ uncertainty is 0.3 dex (Shangguan et al. 2018). (8) Total infrared luminosity of the host galaxy (Shangguan et al.2018). (9) Star formation rate derived from IR luminosity by adopting Eq. 4 of (Kennicutt 1998) and a Kroupa Initial massfunction (Kroupa 2001). (10) Total gas mass inferred from dust mass measurements. The 1 σ uncertainty is 0.2 dex (Shangguanet al. 2018). (11) Black hole mass, estimated by applying the calibration of Ho & Kim (2015) using the AGN monochromaticluminosity at 5100 ˚A (Col. 12) and the H β line width (Shangguan et al. 2018). (12) AGN monochromatic luminosity at 5100 ˚A. Table 2.
ALMA Cycle 6 Observational SetupObject Observation Bandpass & Phase PWV On-source Beam Size Beam Position RMSDate Flux Calibrator Calibrator (mm) Time (min) ( (cid:48)(cid:48) × (cid:48)(cid:48) ) Angle ( ◦ ) (mJy beam − )(1) (2) (3) (4) (5) (6) (7) (8) (9)PG 0050+124 26 Nov. 2018 J0006 − . × .
36 21.1 0.55PG 0923+129 17 Mar. 2019 J1058+0133 J0853+0654 1.57 16.72 1 . × .
00 71.6 0.47PG 1011 −
040 14 Mar. 2019 J1037 − − . × .
91 72.0 0.59PG 1126 −
041 14 Mar. 2019 J1058+0133 J1131 − . × .
89 83.1 0.37PG 1244+026 19 Mar. 2019 J1256 − . × .
17 95.4 0.41PG 2130+099 22 Mar. 2019 J2000 − . × .
13 25.5 0.34 model parameter number, instead of a χ -based crite-rion to determine whether a line is detected. We esti-mate the likelihood of the best-fit Gaussian model bycomparing it with a straight line fit (i.e., no emissionline present). We consider a 5 σ probability thresholdfor detection. If this threshold is not achieved, we binthe spectrum over a larger area by increasing the sizeof the extraction box by one pixel per side, and thenrepeat the fit until either the probability threshold isachieved or the third iteration is reached. If no detec-tion is achieved after three iterations, we assume thatno emission is present, mask the pixel, and skip to thenext one. We stop at the third binning iteration in or-der to avoid large binned zones when compared to thebeam size. We also exclude highly uncertain models by masking pixels that have lines with peak signal-to-noise (S/N) lower than 3. Pixels with high S/N oftenshow asymmetric, occasionally highly complicated lineshapes. When an emission line is detected, we increasethe number of Gaussian components and we repeat thefit. We compute the multi-Gaussian model likelihoodby calculating BIC with respect to the last acceptedmodel. Again, we consider a 5 σ probability thresholdfor model acceptance. Most spectra that are symmetriccan be well fit with a single Gaussian component, thosethat are asymmetric can be described by two Gaussians,and three components are necessary for even more com-plex shapes. Figure 2 gives some examples for the caseof PG 1126 − olecular Gas in Quasar Hosts N o r m a li ze d f l ux d e n s it y Symmetric
DataGaussian N o r m a li ze d f l ux d e n s it y Asymmetric
DataTwo-Gaussian
300 150 0 150 300Velocity [kms − ]0.00.51.0 N o r m a li ze d f l ux d e n s it y Complex
DataTriple-Gaussian
Figure 2.
Example of the three types of CO(2–1) lineprofiles observed for PG 1126 − form better than a high-order Gauss-Hermite function(van der Marel & Franx 1993).The least-squares minimization technique can behighly sensitive to the given initial guesses. To mit-igate this issue, we refit the detected emission lines byconsidering a series of new initial guesses taken from thebest-fit parameters obtained from the neighboring pix-els and from the pixel itself during the first step. Theneighbor pixels are defined as those within the binningbox used to extract the averaged spectrum. We selectthe final model that gives the lowest BIC. Note that thefitting procedure can give rise to false positive detectionsfrom noisy peaks in the spectrum. We mask these noisypixels by applying a procedure that mimics the prun-ing routine employed by the auto-multithresh taskwithin tclean , by masking detected pixel groups thathave projected sizes smaller than 0.6 times the beamsize.Finally, we use Monte Carlo resampling to derive themodel parameter uncertainties. We measure the averagespectrum noise level for each pixel and assume that thenoise follows a normal distribution. We then add thesimulated noise to the observed spectrum and fit theline. We iterate 300 times to obtain a distribution foreach parameter, and we estimate the 1 σ uncertaintiesof the parameter from the 16th and 84th percentiles ofthe distribution. The intensity, line-of-sight (LOS) velocity, and veloc-ity dispersion maps are shown in Figure 3. The lattertwo maps are constructed by calculating the luminosity-weighted moment one and two from the best-fit modelsof the emission line in each pixel.We perform a sanity check to determine how muchCO(2–1) emission is missed by our line shape fittingprocedure. We calculate the total velocity-integratedline intensity ( S CO(2 − ∆ v ) from the intensity mapsand we compare these to the values obtained by sum-ming the flux densities from the channels above the 2 σ level for each target. We find that our emission-line fit-ting procedure recovers ∼
97 % of S CO(2 − ∆ v for thetwo worst cases (PG 0050+124, PG 0923+129), while forthe best case (PG 1011 − ∼ S CO(2 − ∆ v value. These estimates indicate that ouremission line procedure recovers most of the reliableCO(2–1) emission recorded in the datacubes.2.3. CO(2–1) Luminosity and Molecular GasEstimates
We estimate the S CO(2 − ∆ v values by summing allthe pixel values from the intensity maps. The luminos-ity L (cid:48) CO(2 − is calculated following (Solomon & VandenBout 2005) L (cid:48) CO(2 − = 3 . × S CO(2 − ∆ v D L ν (1 + z ) [K km s − pc ] , (1)where S CO(2 − ∆ v is in units of Jy km s − , ν obs is the ob-served frequency of the line in GHz, D L is the luminositydistance in Mpc, and z is the redshift. The CO(2–1) lu-minosities are used to estimate the CO(1–0) luminosityby adopting a luminosity ratio L (cid:48) CO(2 − /L (cid:48) CO(1 − =0 .
62, the median value found by Shangguan et al.(2020a) for eight PG quasar host galaxies. We obtainmolecular gas masses assuming a CO-to-H conversionfactor α CO = 3 . M (cid:12) (K km s − pc ) − with 0.3 dex un-certainty (Sandstrom et al. 2013), a value consistentwith dust-based gas masses independently derived forthe PG quasars (Shangguan et al. 2020a).2.4. Comparison Sample
Our analysis (Section 3.5) will compare the prop-erties of quasar hosts with those of inactive galaxies.We choose, for comparison, the EDGE-CALIFA sample(Bolatto et al. 2017), an interferometric CO(1–0) studywith the Combined Array for Millimetre-wave Astron-omy of 126 nearby (23–130 Mpc) galaxies selected fromthe CALIFA survey (S´anchez et al. 2012). With a spec-tral resolution of ∼
10 km s − and an average spatial res- Molina, J. et al.
PG0050+124HST/WFC3-F105W
Intensity
Residual
LOS Velocity v I.U -0.9 0.9
I.U -240 240 km/s
10 80 km/s
PG0923+129HST/WFC3-F814W
Intensity
Residual
LOS Velocity v I.U -0.4 0.4
I.U -200 200 km/s
10 80 km/s
PG1011-040HST/WFC3-F814W
Intensity
Residual
LOS Velocity v I.U -0.08 0.08
I.U -100 100 km/s
10 60 km/s
PG1126-041HST/NICMOS-F160W
Intensity
Residual
LOS Velocity v I.U -0.3 0.3
I.U -280 280 km/s
10 120 km/s
PG1244+026HST/WFC3-F814W
Intensity
Residual
LOS Velocity v I.U -0.03 0.03
I.U -100 100 km/s
10 40 km/s
PG2130+099HST/WFPC2-F555W
Intensity
Residual
LOS Velocity v I.U -0.4 0.4
I.U -260 260 km/s
10 100 km/s
Figure 3.
Hubble Space Telescope (HST) image taken from Zhao et al. (in prep.), CO(2–1) line intensity, best-fit intensityresiduals (model fits described in Section 3.2), LOS velocity, velocity dispersion, and map of emission-line shapes for the PGquasars. In the first and second columns, the light blue contours are taken from each HST image and correspond to the[0 . , , , , × − levels relative to the maximum brightness value (for PG 2130+099 we show the [3 , . , , , × − relative levels to improve visualization). The HST images reveal the stellar component emission of quasar host galaxies. Thebeam is plotted above the scale bar of the CO(2–1) intensity map. In the LOS velocity map, the dashed line represents thedirection of the kinematic major axis. olecular Gas in Quasar Hosts log( M /M fl ) l og ( M H / M fl ) f H = . f H = . This workShangguan+20aCO − LiteraturexCOLDGASSEDGE-CALIFA1.51.00.50.00.51.01.52.02.5 l og ( S F R / M fl y r − ) + . d e x This workShangguan+20aCO − LiteraturexCOLD GASSEDGE-CALIFAMain sequence
Figure 4.
PG quasar host galaxies compared to the EDGE-CALIFA star-forming galaxies with measured CO sizes (Bo-latto et al. 2017). (Top) The M (cid:63) –SFR plane. We showhost galaxies with CO data, color-coded following Figure 1.The EDGE-CALIFA data are presented in light blue andblue squares, the latter representing the subsample chosento compare with our observations. The grey triangles showthe xCOLD GASS sample (Saintonge et al. 2017). The solidline indicates the main sequence suggested by Saintonge et al.(2016), with the orange shaded band indicating the ± . M (cid:63) – M H plane.The data are color-coded as in the top panel, with the excep-tion that upper limits for M H for the xCOLD GASS galaxiesare indicated by downward arrows. The dashed lines indi-cate values of constant f H . The PG quasar host galaxiesand the EDGE-CALIFA survey overlap in both parameterspaces. olution of ∼ . assuming α CO = 4 . M (cid:12) (K km s − pc ) − , similar tobut somewhat higher than our preferred value of α CO =3 . M (cid:12) (K km s − pc ) − . For consistency with our con-vention, we scale the EDGE-CALIFA molecular gasmasses by a factor of ∼ .
71. We only consider theEDGE-CALIFA galaxies with measured CO sizes (69systems). We further discard five sources that are likelyto be AGN hosts based on the optical line intensity ratiodiagnostics of Baldwin et al. (1981) and H α linewidthanalysis (Lacerda et al. 2020). These five AGN hostsdo not show any particular trend when compared tothe main EDGE-CALIFA sample and only one of these(MRK 79) is classified as type-I AGN (log λL λ (5100˚A) ∼ .
16 erg s − ; Lu et al. 2019). Therefore, our compari-son sample corresponds to a total of 64 galaxies takenfrom the EDGE-CALIFA survey.Figure 4 compares the PG quasar host galaxies (thesix sources mapped by ALMA plus the larger sampleof sources with CO observations from Shangguan et al.2020a) with the EDGE-CALIFA galaxies, in terms oftheir SFR, stellar mass ( M (cid:63) ), molecular gas mass ( M H ),and molecular gas mass fraction [ f H ≡ M H / ( M (cid:63) + M H )]. To improve further the statistics of the in-active galaxies, we also include the larger sample ofnearby galaxies from xCOLD GASS (Saintonge et al.2017). As discussed in Shangguan et al. (2020b, seealso Xie et al. submitted), the PG quasars generallytrack the main sequence of star-forming galaxies as de-fined by Saintonge et al. (2016), with a non-negligiblefraction lying significantly above it, to the extent thatthey can be deemed starburst systems. The ALMAsample, by virtue of their selection, is biased towardhigher SFRs, M H , and f H compared to the overall PGsample with CO observations: all six objects lie on orabove the main sequence, half of them formally exceed-ing the main sequence 1 σ scatter of 0.4 dex. The EDGE-CALIFA subsample overlaps well with the overall sam-ple of quasar hosts in terms of M (cid:63) , SFR, M H , and f H , but for the purposes of achieving a better matchwith the ALMA-mapped quasars, we further distinguishthe EDGE-CALIFA galaxies that lie above the main se-quence (blue squares in Figure 4). ANALYSIS AND RESULTS Bolatto et al. (2017) adopt α CO = 0 . M (cid:12) (K km s − pc ) − forthe ultraluminous infrared galaxy Arp 220. We follow their con-vention. Molina, J. et al.
20 40 60 80 S CO(2 − ∆ v (ACA) [Jy km s − ] S C O ( − ) ∆ v ( m + A C A ) [ J yk m s − ] PG 0050+124PG 0923+129PG 1011-040PG 1126-041PG 1244+026PG 2130+099
Figure 5.
Comparison between the concatenated (12m+ ACA) data and the ACA velocity-integrated CO(2–1)emission-line fluxes published in Shangguan et al. (2020a).The dashed line show the unity ratio. The 1 σ errorbars ofthe data (without including flux calibration uncertainties)are smaller than the symbol sizes. We find good agreementbetween both estimates (∆[ S CO(2 − ∆ v ] (cid:46)
24 %), suggest-ing that the concatenated data can recover well both thecompact and extended CO(2–1) emission.
Sensitivity of the Concatenated Data Compared tothe ACA Data
To verify whether our process of data concatenation(12 m + ACA) properly recovers the CO line emissionon extended scales, Figure 5 compares S CO(2 − ∆ v mea-sured from our intensity maps with the ACA-only valuesreported by Shangguan et al. (2020a). The two sets ofmeasurements show good agreement, but there is a mildtendency for the concatenated data to recover somewhathigher fluxes than the ACA data alone. Four of the sixsources recover extra CO emission, ranging from ∼ ∼ (cid:46) (cid:48)(cid:48) region of the galaxy, whereas our map (Fig-ure 3) shows that the host galaxy exhibits CO spiralarms that extend up to ∼ (cid:48)(cid:48) ( ∼ ∼ ∼
6% flux deficit in the con-catenated data set, but this is insignificant compared tothe ∼ −
10% uncertainty of the absolute flux scale(Fomalont et al. 2014; Bonato et al. 2018). We con-clude that the concatenated data are sensitive enoughto trace most of the CO(2–1) line emission coming fromthe central part of the host galaxies.3.2.
CO(2–1) Intensity Map Modeling
We model the two-dimensional intensity maps using aradial profile described by a S´ersic (1963) function, I ( R ) = I e exp (cid:40) − b n (cid:34)(cid:18) RR e (cid:19) /n − (cid:35)(cid:41) , (2)where b n is a constant that sets R e as the effective (half-light) radius, I e is the intensity measured at R e , and n is the S´ersic index. An exponential function, whichoften describes a cold disk, corresponds to n = 1. Thetwo-dimensional model considers seven free parameters( I e , n , R e , PA pho , b/a minor-to-major axis ratio, [ x , y ]on-sky center location, and background) and the S´ersicmodel is convolved with the synthesized beam to recoverunbiased estimates.To model the galaxy projection on the sky, we assumethat the host galaxies follow an oblate spheroidal geom-etry (Hubble 1926), such thatcos i = ( b/a ) − q − q , (3)where q is the intrinsic thickness. We set q = 0 . z < .
05; Mosenkovet al. 2015; see Yu et al. 2020 for a more complicatedprescription).We use the
Python package emcee (Foreman-Mackey et al. 2013) to find the best-fit model. Briefly, emcee implements the affine invariant ensemble sam-pler for the Markov chain Monte Carlo sampling methodto characterize the model probability density function.This allows the algorithm to perform equally well un-der affine transformations (including linear transforma-tions) between the model parameters, and therefore isless sensitive to the possible covariances among them(see Foreman-Mackey et al. 2013 for more details). Weoptimize the log-likelihood olecular Gas in Quasar Hosts I n t e n s i t y [ I . U . ] PG 0050+124 modeldata10 0 10
Radius [kpc] -0.1000.0000.100 R e s i d u a l s I n t e n s i t y [ I . U . ] PG 0923+129 modeldata10 5 0 5 10
Radius [kpc] -0.0500.0000.050 R e s i d u a l s I n t e n s i t y [ I . U . ] PG 1011-040 modeldata5 0 5
Radius [kpc] -0.0100.0000.010 R e s i d u a l s I n t e n s i t y [ I . U . ] PG 1126-041 modeldata10 5 0 5 10
Radius [kpc] -0.0500.0000.050 R e s i d u a l s I n t e n s i t y [ I . U . ] PG 1244+026 modeldata5.0 2.5 0.0 2.5 5.0
Radius [kpc] -0.0050.0000.005 R e s i d u a l s I n t e n s i t y [ I . U . ] PG 2130+099 modeldata10 0 10
Radius [kpc] -0.0500.0000.050 R e s i d u a l s Figure 6.
Intensity profiles across the major axis of the CO(2–1) image, compared with model profiles extracted from thebest-fit two-dimensional S´ersic models. We find good agreement between the data and models, as the residual panels suggest. ln L ≡ − N (cid:88) i (cid:18) ( y i − y mi ) σ i + ln(2 πσ i ) (cid:19) , (4)where y i denotes the pixel data taken from the intensitymaps, σ i is the 1 σ uncertainty, and y mi corresponds tothe model value for that pixel.The intensity maps, along with the residuals from thebest-fit model, are shown in Figure 3, while the best-fit parameters are presented in Table 3. These best-fitparameters are also used to extract the radial intensityprofile for each system (Figure 6), which is derived bysimulating a slit with width equal to the beam FWHMand aligned with respect to the major axis of the CO(2–1) image (given by PA pho ). Across the simulated slit thedata are sampled in radial bins with size equal to halfof the beam FWHM to avoid over- and under-sampling.For each radial bin, we calculate the average intensityvalue and the 1 σ uncertainty, which is given by the stan-dard deviation of the encompassed individual pixel val-ues.The stellar component of the six quasar host galax-ies are all observed by HST with high spatial res-olution at optical or NIR wavelengths. HST im-ages of PG 0050+124, PG 0923+129, PG 1011 − I -band (PI: L. C. Ho). Images of PG 1126 − H -band; PI:S. Veilleux) and WFPC2 ( B -band; PI: S. R. Heap). We collect them in Figure 3 as a reference for the CO emis-sion. The molecular gas tends to be distributed withinthe stellar component delineated by the HST images.In four sources (PG 0050+124, PG 0923+129,PG 1126 − ∼
10 kpc from the center, although thebulk of it is confined to a more compact, disk-like mor-phology. The CO morphology of PG 1244+026 is themost compact, with a total radial extension (cid:46) −
040 exhibits a complex structure.Appendix A gives comments on the morpho-kinematicsof each source.Despite the observed variety of CO morphologies,the two-dimensional flux distributions are well fittedby S´ersic models with indexes in the range of n ≈ . − .
9. Subtracting the best-fitting global componentfrom each map reveals complex sub-structures in theresidual maps, ranging from clumps (PG 1126 −
041 andPG 1011 − − Kinemetry Molina, J. et al.
Table 3.
CO(2–1) Emission and Spatial DistributionObject S CO(2 − ∆ v log L (cid:48) CO(2 − I e n R e b/a PA pho R / (Jy km s − ) (K km s − pc ) (Jy km s − ) (kpc) ( ◦ ) (kpc)(1) (2) (3) (4) (5) (6) (7) (8) (9)PG 0050+124 71 . ± .
08 9 . ± .
001 0 . ± .
001 1 . ± .
02 1 . ± .
01 0 . ± .
01 15 ± . ± . . ± .
06 8 . ± .
001 0 . ± .
001 0 . ± .
01 1 . ± .
01 0 . ± .
01 53 ± . ± . −
040 18 . ± .
04 8 . ± .
001 0 . ± .
002 0 . ± .
01 1 . ± .
01 0 . ± .
01 102 ± . ± . −
041 17 . ± .
04 8 . ± .
001 0 . ± .
001 0 . ± .
03 3 . ± .
01 0 . ± .
01 322 ± . ± . . ± .
02 8 . ± .
001 0 . ± .
001 0 . ± .
01 1 . ± .
01 0 . ± .
01 186 ± . ± . . ± .
04 8 . ± .
001 0 . ± .
001 1 . ± .
01 1 . ± .
01 0 . ± .
01 44 ± . ± . Note — (1) Source name. (2) Velocity-integrated flux density estimated from our concatenated data. We have not consideredthe ALMA flux calibration uncertainty ( (cid:46)
10 %; Fomalont et al. 2014; Bonato et al. 2018). (3) CO(2–1) line luminosity, updatedfrom the concatenated data measurements. (4) Velocity-integrated line intensity at the effective radius. (5) S´ersic index. (6)Effective radius. (7) Ratio of the projected minor axis to to major axis. (8) Position angle (North = 0 ◦ , East = 90 ◦ ). (9)CO(2–1) half-light radius calculated by implementing a tilted-ring approach and using the best-fit PA pho and b/a values (Col.7 and 8). Radius [kpc] -3 -2 -1 k v , / k v , Krajnovic+08
PG0050+124PG0923+129 PG1011-040PG1126-041 PG2130+099
Figure 7.
Radial variation of the kinemetry coefficient k v, /k v, for five PG quasar host galaxies. The dot-dashedline corresponds to (cid:104) k v, /k v, (cid:105) = 0 .
02 threshold that sep-arates between “regular” and “non-regular” rotators (Kra-jnovi´c et al. 2008). The colored arrows indicate the hostgalaxy radius at which k v, /k v, is consistent with the thresh-old value. We find that four host galaxies are classified asregular rotators, while only one system, PG 1011 − We investigate whether the observed molecular gas ro-tational motions resemble those expected from an idealdisk by analyzing the two-dimensional LOS velocitymaps using the tilted-ring approach as implemented in kinemetry (Krajnovi´c et al. 2006) . kinemetry quan-tifies the deviations from disk-like kinematics for an ob-served velocity field by parameterizing it as a functionof the radius ( R ) and the azimuthal angle ( ψ ). For anideal disk, the velocity profile is only a function of thegalaxy radius, and the sky projection adds a variationin terms of azimuthal angle that follows a cosine law: v ( R, ψ ) = v ( R ) cos( ψ ). To test this condition, kineme-try decomposes the LOS velocity into a series of tiltedrings. Along each elliptical path, it parameterizes thevelocity profile in terms of a Fourier series that onlydepends on the azimuthal angle: v LOS ( a, ψ ) = A ( a )+ N (cid:88) n =1 k v,n ( a ) cos[ n ( ψ − φ n ( a ))] , (5)where k v,n and φ n are the amplitude and phase coeffi-cients, and a is the length of the kinematic major axis.The kinematic position angle (PA kin ) is accurately re-covered in tilted rings dominated by a single component,whereas in the ellipses between multiple components ittraces the position of maximum velocity amplitude. Thecase of an ideal disk is recovered when all the k v,n coef-ficients are zero except k v, . Hence, the high-order am-plitude coefficients quantify the kinematic deviations or http://davor.krajnovic.org/idl/ olecular Gas in Quasar Hosts kinemetry does not fit the kinematic center, and thisparameter must be given in advance. We use the centerdetermined by the best-fit S´ersic model as input. Wealso fix the axis ratio equal to the value given by ourmorphological modeling, as this parameter is not well-constrained by fitting the LOS velocity map alone. Thetilted ring thickness is set equal to half of the beamFWHM, and we skip the first tilted ring iteration ineach map because the radius is smaller than the beamFWHM. We restrict the analysis to the central zone ofthe galaxy, which, in any case, is our primary interest.We find that kinemetry does not interpolate reliablythe maps in the outer regions, where spiral arms maybe present (e.g., PG 0050+124), producing spurious re-sults such as inverted kinematic position angles at largeradii. This restriction also helps to avoiding outer zonesthat may be under-sampled due to incomplete azimuthalcoverage. For each tilted ring, we use the default 75%pixel sampling limit required by kinemetry to obtainreliable estimates (Krajnovi´c et al. 2006).To consider the uncertainties associated with the kine-matic center and the LOS velocity map, we bootstrapthese values within their 3 σ error range and repeat the kinemetry procedure 100 times. Then, the 1 σ uncer-tainties of the k n,v amplitude coefficient are estimatedfrom the corresponding distribution by calculating the16th and 84th percentiles.We quantify the non-regular motions by calculatingthe ratio k ,v /k ,v (e.g., Krajnovi´c et al. 2008, 2011; vande Sande et al. 2017). The fifth-order amplitude co-efficient is used because kinemetry adjusts the lowerorder coefficients to find the best-fit tilted ring. Kra-jnovi´c et al. (2011) used an average threshold ratioof (cid:104) k ,v /k ,v (cid:105) = 0 .
04 to classify galaxies as “regular”or “non-regular” rotators for the the ATLAS sur-vey (Cappellari et al. 2011). However, this thresholddepends on data quality measured by the typical un-certainty in k ,v /k ,v . For example, Krajnovi´c et al.(2008) used a limit of (cid:104) k ,v /k ,v (cid:105) = 0 .
02 to determine ifa galaxy presents ‘disk-like rotation’ for the SAURONproject (de Zeeuw et al. 2002). Based on the typi-cal uncertainty of k ,v /k ,v measured for our systems( ∼ . (cid:104) k ,v /k ,v (cid:105) = 0 .
02 forthis study. We apply the kinemetry analysis for five of the six PG quasars . The radial median of k v, /k v, spans 0 . − .
027 (Table 4), such that four would beclassified as regular rotators, with k v, /k v, < .
02 overalmost all radii, with no clear trend toward the inneror outer regions. PG 0050+124 and PG 1126 −
041 havelocal values of k v, /k v, ≈ .
02 (as indicated by the ar-rows in Figure 7), but they are localized galaxy sub-structures observed in the intensity residual maps. OnlyPG 1011 −
040 shows k v, /k v, > .
02 at R (cid:38) . kin ) with respect to themedian value ( (cid:102) PA kin ) and as a function of radius forthe host galaxies. We find that PA kin varies smoothlywith radius in each source, with a maximum differenceof ∆PA kin ≈ ◦ for PG 1011 − − − kin > ◦ ac-cording to the criteria of Krajnovi´c et al. (2008). Foreach galaxy, we assume that the direction of the kine-matic major axis is given by (cid:102) PA kin (Table 4). The (cid:102) PA kin values differ on average by ∼ ◦ from the PA pho val-ues determined by the S´ersic models. PG 0050+124shows the higher difference between both estimates( (cid:102) PA kin − PA pho ∼ ◦ ), whereas PG 1126 −
041 presentsthe smaller difference ( (cid:102) PA kin − PA pho ∼ ◦ ). Similardifferences are found when we compare (cid:102) PA kin with thestellar photometric position angles measured from theHST images (Veilleux et al. 2009; Zhao et al. in prep.).In this case, we find differences in the range of 0 − ◦ ,with an average difference of ∼ ◦ .To understand the difference between (cid:102) PA kin andPA pho values, we also analyze the PG quasars CO(2–1)intensity maps using kinemetry . In this case, kineme-try performs a simple ellipse fitting method. We applythe same restrictions that we assumed for the LOS ve-locity map analyzes. In the bottom panel of Figure 8we show the difference of the position angles derivedfrom the intensity maps (PA int ) with respect to (cid:102) PA kin .We measure larger radial variations of PA int (∆PA int ≈ ◦ ) when compared to PA kin (∆PA kin ≈ ◦ ), high-lighting the effect of molecular gas sub-structure onthe intensity maps, but only relatively subtle finger- We avoid analyzing PG 1244+026 because of the high degree ofcompactness of the system. With a projected major axis exten-sion of ∼ (cid:48)(cid:48) , this source is only sampled by ∼ ∼ ∼ . (cid:48)(cid:48) Molina, J. et al. P A k i n − g P A k i n [ ◦ ] PG0050+124PG0923+129PG1011-040 PG1126-041PG2130+099
Radius [kpc] P A i n t − g P A k i n [ ◦ ] PG0050+124PG0923+129 PG1011-040PG1126-041 PG2130+099
Figure 8.
Radial variation of the position angle measuredfrom the LOS velocity maps (PA kin ; Top ) and intensity maps(PA int ; Bottom ) for five PG quasars. In both panels, the PAvalues are presented with respect to the median kinematicposition angle (cid:102) PA kin . The data are color-coded followingFigure 7. The horizontal dotted lines indicate the differenceof the global PA pho value with respect to (cid:102) PA kin . prints on the LOS velocity fields. With the exception ofPG 0923+129, PA int approximates to PA pho (horizontaldotted lines in Figure 8) at longer radii. On the otherhand, PA int tends to be consistent with (cid:102) PA kin at smallradii (except for PG 2130+099). The difference between (cid:102) PA kin and PA pho may be produced by the large variationof the photometric position angle values tracing multiplemolecular gas components induced by, perhaps, stellarbars (e.g., PG 1011 − s and s , respectively) maybe indicative of elliptical streaming or radial flows. El-liptical streaming generally produces an anti-correlationbetween the s and s terms, while axisymmetric ra-dial flows can be identified in the case of significant s value but negligible s term ( | s /s | (cid:38) P r o f il e [ k m s ] PG 0050+124 v LOS v P r o f il e [ k m s ] PG 0923+129 v LOS v P r o f il e [ k m s ] PG 1011-040 v LOS v P r o f il e [ k m s ] PG 1126-041 v LOS v P r o f il e [ k m s ] PG 2130+099 v LOS v Figure 9.
LOS velocity and velocity dispersion profilesacross the kinematic major axis for five PG quasars. Theshaded zones represent the inner regions that suffer fromsignificant beam-smearing effect. PG 0050+124 shows con-siderably elevated σ v values in its central zone ( R (cid:46) . olecular Gas in Quasar Hosts
13a static bar potential can produce the same kinematicsignature of a radial flow (see Wong et al. 2004, for moredetails). We obtain the s and s radial profiles directlyfrom kinemetry for the host galaxies. The s and s coefficients absolute values are within (cid:46) − , withno clear evidence of gas elliptical streaming or axisym-metric radial flows in any host galaxy. It is worth tomention that results obtained from the harmonic de-composition analyses of the velocity fields must not beconsidered as conclusive. Kinematic features producedby non-axysimmetric sub-structures in disk galaxies lim-itate further interpretation (Wong et al. 2004). More-over, radial flows signatures suggested by gravitationaltorque modeling do not necessarily coincide with theexpectation from the purely kinematic decompositionanalysis (Haan et al. 2009). Any interpretation of theobserved kinematics requires a detailed knowledge of thegalaxy potential (e.g. Alonso-Herrero et al. 2018).3.4. Regular Motions and Velocity Dispersion
Rotational patterns can be seen in the LOS velocitymaps, while velocity dispersion rises toward the center(Figure 3). To further study the kinematics, we extractvelocity radial profiles along the kinematic major axisto avoid any effect produced by the disk azimuthal pro-jection. These velocity radial profiles, available for fiveof the six sources (Figure 9), are derived in the samemanner as the intensity radial profiles, but we simu-late a slit aligned with the kinematic major axis insteadof the photometric major axis. We can measure rota-tion curves to (cid:38) − − ;the masked region is denoted by the shaded blue area inFigure 9. Representative values of the rotation velocity( v rot ; Table 4) are estimated by the mean value of theunmasked data points on the rotation curve, and the 1 σ uncertainty is estimated from their standard deviation. Table 4.
Kinematic ParametersObject (cid:102) PA kin k v, /k v, v rot σ v ( ◦ ) (km s − ) (km s − )(1) (2) (3) (4) (5)PG 0050+124 304 0.010 337 ±
92 36 ± ±
20 6 ± −
040 151 0.027 141 ± ± −
041 326 0.010 275 ± ± ±
15 12 ± Note — (1) Source name. (2) Median position angle of thekinematic major axis (with respect to the receding side) mea-sured from kinemetry analysis (North = 0 ◦ , East = 90 ◦ ).(3) Median value of the ratio of the fifth-order amplitudecoefficient over the first-order coefficient; the typical 1 σ un-certainty is 0.001. (4) Average rotation velocity derived fromthe rotation curve, corrected for inclination. (5) Average ve-locity dispersion derived from the velocity dispersion profile. We correct these values for inclination projection. Ananalogous treatment is applied to the velocity disper-sion ( σ v ), using the same mask applied to the rotationvelocity curve so as to avoid the central pixels affectedby beam-smearing (e.g., Davies et al. 2011; Wisnioskiet al. 2015; Stott et al. 2016).We find v rot /σ v = 9–42 for the five analyzed hostgalaxies, indicating that the gas kinematics are dom-inated by rotation. This is consistent with the rangemeasured for the EDGE-CALIFA galaxies ( v rot /σ v =10–28; Levy et al. 2018). However, different from theEDGE-CALIFA systems, σ v varies considerably fromobject to object. For example, the average dispersionof PG 0050+124 ( σ v = 36 ±
27 km s − ) is 6 times higherthan that of PG 0923+129 ( σ v = 6 ± − ). The av-erage σ v values reported for the EDGE-CALIFA galax-ies are in the range of ∼ −
19 km s − (Levy et al.2018). The average value of PG 0923+129 is consistentwith the estimate traditionally adopted for nearby star-forming galaxies ( ∼ − ; Leroy et al. 2008). Onthe other hand, PG 0050+124 presents an average σ v value larger than any of those reported for the EDGE-CALIFA systems. These values are calculated over por-tions of the dispersion velocity profile that should beimmune from beam-smearing, and thus these may re-flect different physical properties of each host galaxy ordegree of AGN activity. Detailed study of the molecu-lar gas velocity dispersion will be addressed in a futurework. 3.5. Compact Molecular Gas Distribution
In terms of global quantities, the interstellar mediumcontent of the PG quasars seems to be similar to that of4
Molina, J. et al.
Figure 10.
Distribution of molecular gas mass versus the half-light radius of the CO distribution, as measured for the six PGquasars and the star-forming galaxies in the EDGE-CALIFA survey. The EDGE-CALIFA data are presented in light blue andblue squares, where the latter symbol represents the subsample that better matches the PG quasars in terms of SFR, M (cid:63) , M H ,and f H . The dashed curves correspond to values of constant molecular gas surface density measured in [ M (cid:12) pc − ] units. Theblack bar in the lower-right corner illustrates the variation in M H induced by the choice of α CO : the range above and belowthe black dot represent the effect of varying α CO between the value for the Milky Way and ultra-luminous infrared galaxies.The filled orange symbols indicate the galaxies that belong to a multiple system. The observed PG quasars tend to have morecompact molecular gas distributions compared to the EDGE-CALIFA galaxy subsample. normal star-forming galaxies of the same stellar mass(Shangguan et al. 2018, 2020a). This is emphasizedagain in Figure 4, where we show that the majority ofthe host galaxies with CO measurements have similar M H , f H , M (cid:63) , and SFR relative to inactive galaxies ob-served by the xCOLD GASS (Saintonge et al. 2017) andEDGE-CALIFA surveys (Bolatto et al. 2017). However,one missing key information is the spatial distribution ofthe molecular gas, which now can be investigated withour new observations.Figure 10 shows the molecular gas content of the PGquasars plotted as a function of the half-light radius( R / ) measured from the intensity maps by implement-ing a tilted ring approach using the best-fit S´ersic modelparameters. We use R / instead of R e for consistencywith the methodology employed by Bolatto et al. (2017)for the EDGE-CALIFA galaxies, which serve as a di-rect comparison. The blue squares in Figure 10 high-light the subset of EDGE-CALIFA galaxies that bettermatches the quasars in terms of their global physicalproperties (Figure 4). The R / estimates agree with the R e values for the six hostgalaxies (Appendix B). The six quasar hosts tend to have compact molecu-lar gas distributions: R / = 0 . − . H ≡ M H /π (2 R / ) , the sample ischaracterized by Σ H (cid:38) M (cid:12) pc − (see dashed curvesin Figure 10). The PG quasars tend to be located atthe high-Σ H , small- R / region of the CO mass-sizeplane compared to the EDGE-CALIFA galaxies. Wenote, however, that a minority of inactive galaxies alsohave comparable Σ H and R / , and compact molecu-lar gas distribution is not a property unique to AGNs. We further distinguish between the galaxies that belongto a multiple system (filled orange symbols). This isbased on the visual detection of any system (indepen-dent of its size or brightness) within the HST imagefield-of-view for the PG sources (Figure 3) and litera-ture morphological classifications of the EDGE-CALIFAgalaxies, assembled by Bolatto et al. (2017). As Bolattoet al. (2017) have noted, EDGE-CALIFA galaxies in amultiple system tend to be more compact. The sample The five AGN hosts previously discarded from the EDGE-CALIFA sample (not shown in Figure. 10) have R / ∼ − H ∼ M (cid:12) pc − . olecular Gas in Quasar Hosts PG0050+124
PG0923+129
PG1011-040
PG1126-041
PG1244+026
PG2130+099
Figure 11.
Maps of the locations of spectra with shapesclassified as “symmetric” (dark blue), “asymmetric” (green),and “complex” (yellow), according to the multi-Gaussian fit-ting procedure described in Section 2.2. In each panel, thesynthesized beam is shown in the bottom right corner. of quasars (6 out of 40 sources at z < .
3) is too smallto draw any meaningful conclusions in this regard.Our statements on Σ H obviously depend on thechoice of α CO . This is especially worrisome if α CO dropsin environments of high surface density (Bolatto et al.2013). Shangguan et al. (2020a), in comparing COobservations of the PG quasars with independent gasmasses inferred from far-infrared dust emission, find noevidence that quasar host galaxies possess abnormal val-ues of α CO as it is found for luminous and ultra luminousinfrared galaxies ( ∼ . M (cid:12) (K km s − pc ) − ; Herrero-Illana et al. 2019). In any event, our main conclusionsdo not qualitatively change within the likely range of α CO (0.8–4.36 M (cid:12) (K km s − pc ) − , vertical bar in Fig-ure 10).3.6. Asymmetric and Complex Line Shapes -300 -150 0 150 300
Velocity [km s ] N o r m a li z e d f l u x d e n s i t y v v v v + v v = v + + v DataModel
Figure 12.
Example of the non-parametric approach im-plemented to characterize the emission line. The quantities v , v , and v correspond to the velocity values of the 5th,50th, and 95th percentiles estimated from the cumulativedistribution of the emission-line flux density. The quantities∆ v + and ∆ v − correspond to the velocity difference of v and v with respect to v , respectively. Thus, ∆ v + is posi-tive and ∆ v − is negative. A symmetric line is characterizedby ∆ v = 0. As discussed in Section 2.2, the pixelwise emissionlines come in three generic shapes that can be decom-posed into one, two, or three Gaussians, which we call“symmetric,” “asymmetric,” and “complex,” respec-tively. Figure 11 gives the spatial distribution of theline shapes, constructed from the best-fit models of theemission line in each pixel. Symmetric lines are gen-erally located at large radii, while complex profiles areconcentrated toward the central regions of the galaxies,often associated with sub-structures in the gas distribu-tion.To investigate the origin of the asymmetric and com-plex profiles more quantitatively, we implement a non-parametric scheme to quantify the line profile, moti-vated by procedures commonly used to study the sig-natures of ionized gas outflows in galaxy spectra (e.g.,Liu et al. 2013; McElroy et al. 2015; Harrison et al. 2016;Sun et al. 2017; Husemann et al. 2019). To character-ize the line asymmetry, we first measure the velocitiesat which 5% ( v ), 50% ( v ) and 95% ( v ) of the totalline flux is enclosed. We compute the velocity differences∆ v + = v − v and ∆ v − = v − v , which describethe line flux distribution on each side of v . Note that,by construction, ∆ v + is always positive and ∆ v − is al-ways negative (Figure 12). Then, the line asymmetryis given simply by ∆ v = ∆ v + + ∆ v − . A line is asym-metric when ∆ v is negative or positive, which meansthat most of the line flux is, respectively, blueshifted orredshifted relative to its centroid. For a symmetric line,∆ v = 0. We calculate ∆ v from the noiseless, best-fitemission-line model.6 Molina, J. et al. - -
200 100 0 100 200 v LOS [km s ] v [ k m s ] r = 0.2 p < 0.01 -150 150 km/s - -
100 0 100 200 v LOS [km s ] v [ k m s ] r = 0.5 p < 0.01 -150 150 km/s - - -
40 20 0 20 v LOS [km s ] v [ k m s ] r = 0.6 p < 0.01 -50 50 km/s -
200 100 0 100 200 300 v LOS [km s ] v [ k m s ] r = 0.5 p < 0.01 -250 250 km/s - -
60 40 20 0 20 40 v LOS [km s ] v [ k m s ] r = 0.8 p < 0.01 -50 50 km/s - -
200 0 200 v LOS [km s ] v [ k m s ] r = 0.5 p < 0.01 -200 200 km/s Figure 13.
Asymmetry maps of the CO(2–1) emission line for the PG quasars. The iso-velocity contours are taken from theLOS velocity maps and labeled in units of km s − , consistent with colorbar. The thick dashed line represents the kinematicmajor axis, when available. The synthesized beam is shown above the scale bar. To the right of each map, we show thecorresponding pixel-by-pixel line asymmetry (∆ v ) plotted against the LOS velocity ( v LOS ); ∆ v inversely correlates with v LOS .We give Pearson’s correlation coefficient r and the p -value in the top-right corner. Figure 13 shows that ∆ v varies smoothly with theLOS velocity, from negative to positive values across themaps and along the kinematic major axis of the system.Central ∼ kpc-scale features are seen in PG 0923+129,PG 1126 −
041 and PG 2130+099 systems. These seemto spatially correlate with the pixels that show complexemission-line profiles (Figure 2) and may be produced bythe presence of a central unresolved molecular gas sub-structure. In all cases ∆ v correlates inversely with v LOS ,with Pearson’s coefficient r = − . − . p < . DISCUSSIONAlthough our sample is admittedly small and possiblybiased (Figure 1), it nevertheless offers a first glimpseinto the properties of the molecular gas in the host olecular Gas in Quasar Hosts z < .
5) PG quasars are characteristically gas-rich(Shangguan et al. 2018, 2020a) galaxies forming starswith an efficiency equal to or even exceeding that ofstar-forming galaxies on the main sequence (Shangguanet al. 2020b; Xie et al. submitted). What is respon-sible for the apparent coeval episodes of vigorous BHaccretion and star formation activity? The external en-vironment provides no obvious clue. While the stellarmorphologies of some host galaxies suggest that theyhave experienced tidal interactions or recent merger ac-tivity, not all hosts with enhanced star formation showevidence of dynamical perturbations (Shangguan et al.2020b; Xie et al. submitted). The current ALMA sam-ple exemplifies the problem well. Of the three quasarsthat clearly lie above the 1 σ scatter of the star-forminggalaxy main sequence (see top panel of Figure 4), onlyPG 0050+124 (I Zw 1) is known to belong to a multiplesystem (Lim & Ho 1999; Scharw¨achter et al. 2003). Themorphologies of PG 1126 −
041 and PG 2130+099 resem-ble those of regular systems. Their CO surface bright-ness distribution closely follow an exponential disk dis-tribution as the best-fit S´ersic indexes suggest. Theirvelocity fields (third column of Figure 3) show little evi-dence of perturbation (Figure 7), further confirming theregular disk-like kinematics.The present high-resolution ALMA observations fur-nish some potentially useful insights. First and fore-most, we find that most of the gas resides in a com-pact, central disky (rotationally dominated) structurewith half-light radii ∼ ∼ . − H (cid:38) M (cid:12) pc − ; Downes & Solomon 1998; Ionoet al. 2009; Wilson et al. 2019), as well as in the hostsgalaxies of infrared-bright quasars (CO emission exten-sion ∼ − not infrared-selected nor doour selection criteria strongly biases our sample in termsof L IR values (or SFRs, Figure 4), and none is experi-encing a major merger event or strong tidal interactions.Apart from the compact dimensions and high molec-ular gas mass surface densities, another intriguing prop- erty of the sample is the apparent misalignment betweenthe global kinematic and photometric axes of the gas,which appears to be significant in four out of five sources,as well as the detection of significant smooth variationof the kinematic position angle with radius (kinematictwisting, Krajnovi´c et al. 2008) in three of the sources(see also Ramakrishnan et al. 2019). Given the lack ofobvious signs of recent merging events or tidal interac-tions, the origin of these kinematic features is unclear.How common are these signatures? Do they play a rolein fueling the nucleus? Much better statistics are re-quired, as well as observations of more luminous AGNsand more massive BHs to extend the relevant parameterrange, to ultimately test against numerical simulations(e.g., Weinberger et al. 2017; Terrazas et al. 2020). Com-parison with a proper control sample of inactive galaxiesis also essential. CONCLUSIONSWe present new Cycle 6 ALMA observations trac-ing the CO(2–1) emission line from six nearby ( z (cid:46) .
06) Palomar-Green quasars selected from our previ-ous Cycle 5 ACA program (Shangguan et al. 2020a,b).The host galaxies have normal molecular gas content( M H ≈ . − . M (cid:12) ) for their stellar masses( M (cid:63) ≈ . − . M (cid:12) ) and star formation rates (SFR ≈ − M (cid:12) yr − ). The ALMA observations were de-signed to resolve the molecular gas in the host galax-ies on scales of (cid:46) ≈ . (cid:48)(cid:48) − . (cid:48)(cid:48) • The quasar host galaxies tend to have disk-like(S´ersic index n = 0 . • The velocity field of the molecular gas in five ofthe quasar hosts is dominated by regular rotation,with v rot /σ v (cid:38)
9, and σ v = 6 −
36 km s − . Theremaining system (PG 1244+026) is too compactto be analyzed with confidence. Four host galaxiesexhibit a flat rotation curve out to radii (cid:38) −
10 kpc from the nucleus; however, one of thesesystems (PG 1011 − • Among the five objects amenable to tilted-ringanalysis, the kinematic position angle deviates8
Molina, J. et al. from the photometric angle on average by ∼ ◦ .Three sources show evidence of a significant kine-matic twist. • The pixelwise spectra show a diversity of lineshapes, from symmetric to asymmetric and morecomplex profiles. We argue that asymmetric pro-files predominantly arise from beam-smearing ef-fects, while complex lines are predominantly as-sociated with gas sub-structures in the central re-gions of the galaxies. ACKNOWLEDGMENTSWe thank to an anonymous referee for thoughtfulcomments and suggestions. We acknowledge supportfrom the National Science Foundation of China grant11721303 and 11991052 (LCH, RW), the National KeyR&D Program of China grant 2016YFA0400702 (LCH,RW), ANID-Chile grants Basal AFB-170002 (FEB, ET),FONDECYT Regular 1160999 (ET), 1200495 (FEB,ET) and 1190818 (ET, FEB), and Anillo de ciencia ytecnolog´ıa ACT1720033 (ET), and Millennium ScienceInitiative ICN12 009 (FEB). We thank Jing Wang foruseful insights and comments. We thank MinJin Kim forsharing his stellar photometric model for PG 2130+099source. This paper makes use of the following ALMAdata: ADS/JAO.ALMA
Software:
Astropy (Astropy Collaboration et al.2013; Price-Whelan et al. 2018),
CASA (McMullin et al.2007), emcee (Foreman-Mackey et al. 2013), kineme-try (Krajnovi´c et al. 2006), lmfit (Newville et al. 2014), matplotlib (Hunter 2007), numpy (Oliphant 2006), scipy (Virtanen et al. 2020).APPENDIX A. NOTES FOR INDIVIDUAL HOST GALAXIES
PG 0050+124:
The CO(2–1) observation showsa compact central distribution along with twoextended molecular gas outer spiral arm struc-tures that match the morphology in the restframe I -band HST/WFC3 image. The LOS velocitymap clearly shows a disk-like rotational patternwith projected peak-to-peak rotational velocity v rot sin i ≈
440 km s − . The detection of the outerspiral structures allow us to trace the rotation ve-locities at R (cid:38) R ≈ . . PG 0923+129:
The CO(2–1) line emission isdistributed along a central disk-like zone from which two molecular gas spiral arms extend to-ward the outskirts of the HST broad-band im-age surface brightness distribution. We detect aninner ring-like structure. The host galaxy dis-plays a rotational pattern, and we observe a ve-locity profile that extends up to the flat part ofthe rotation curve with projected peak-to-peak v rot sin i ≈
330 km s − . The velocity dispersionsincrease mainly smoothly toward the galactic cen-ter, and a sharp increase is seen at R ≈ . PG 1011 − The complex CO(2–1) morphol-ogy can be roughly separated into a central bulge-like component surrounded by a ring-like, clumpystructure. The LOS velocity map shows complexdynamics with a smooth variation of the veloc-ity field across the field-of-view. This suggeststhat the two main structures are spatially con-nected. Comparison with the restframe I -bandHST/WFC3 image reveals a possible spatial corre- olecular Gas in Quasar Hosts PG 1126 − A clear disk-like morphology androtational pattern can be seen. The central CO(2–1) emission seems to be distributed in a bar-likestructure. The system also shows clumpy mor-phology. The intensity map is spatially correlatedwith the HST/NICMOS image, suggesting thatthe stellar and molecular gas components followsimilar distribution. The CO emission traces theflat part of the rotation curve, from which we mea-sure a projected peak-to-peak rotational velocityof v rot sin i ≈
390 km s − . The velocity disper-sion increases smoothly toward the galactic center,reaching ∼
160 km s − . PG 1244+026:
The CO(2–1) observation showsa compact molecular gas distribution with iso-velocity contours that depart somewhat from purerotation. The velocity dispersion map shows veryhigh values across most of the disk, suggestingstrong beam smearing due to the compactness ofthe gas distribution.
PG 2130+099:
The regular gas distributionshows an outer spiral-like structure that matchesthe morphology seen in the HST/WFC3 broad-band image. Disk-like rotation can be seen fromthe LOS velocity map, which is sufficiently ex-tended to allow the detection of the flat part of therotation curve. We measure a projected peak-to-peak rotational velocity of v rot sin i ≈
470 km s − from the galaxy outskirts. The outer sub-structureseen at R (cid:38)
10 kpc presents systematically lowervelocity dispersion compared to the central, moredisky region. B. COMPARISON OF RADIAL SIZESFigure 14 compares the half-light radius obtained fromthe best-fit S´ersic model (i.e, the “effective radius” R e )and the estimates obtained using the tilted-ring ap-proach ( R / ) following Bolatto et al. (2017). We find agood agreement between both measurements for the sixhost galaxies. This suggests that the main CO(2–1) lineemission components closely resemble those expectedfrom the ideal axisymmetric two-dimensional S´ersic pro-files. R e [kpc] R / [ kp c ] PG 0050+124PG 0923+129PG 1011-040PG 1126-041PG 1244+026PG 2130+099
Figure 14.
Comparison between the half-light radius de-rived from the tilted-ring approach and the effective radiusobtained from the best-fit S´ersic model. We find a goodagreement between both estimates.C.
BEAM SMEARING EFFECT ON THEOBSERVED EMISSION LINE SHAPESWe simulate a galaxy observation datacube to exem-plify the effect of beam-smearing on the observed emis-sion line shapes. The simulated galaxy has an expo-nential surface density profile with no dark mater halocontribution. The half-light radius is set to 1.5 kpc fol-lowing our estimates reported in Table 3 for the hostgalaxies. We assume intrinsic Gaussian emission lineshapes with line widths equal to 20 km s − (Table 4)for simplicity. For sky-projection effects, we consider aredshift of 0.06, b/a = 0 . x -axis. Following ourobservational setup, we construct the datacube by set-ting the channel resolution equal to 11 km s − and thepixel size to 0 . (cid:48)(cid:48) x × FWHM y = 1 (cid:48)(cid:48) × . (cid:48)(cid:48)
1. We choose thisbeam setup to simulate the beam-smearing effect onlyalong the galaxy major axis direction and to improve thevisualization of the intrinsic emission lines contributionto the resultant asymmetric emission line shape.Figure 15 shows our simulated galaxy observation andthe emission line shape extracted from the green-coloredpixel. The luminosity-weighted convolution producesthe emission line asymmetry as the relative contribu-tion of the brighter neighbor pixels (representative ofthe systemic velocity) is greater than that of the fainterones. The convolved emission line tail is observed at0
Molina, J. et al.
ABC DEF
Mock observation 100 0 100Velocity [km s ] 0.00.51.0 N o r m a li z e d U n i t s Convolved lineSelf-contribution
A B C DEF
Figure 15.
Example of the beam-smearing effect on theobserved emission line shapes.
Left:
Map of the mock galaxyobservation. The green square indicates the pixel from whichwe extract the convolved emission line profile. The purplesquares (listed by letters) correspond to the neighbor pixelsthat mainly contribute to the asymmetric line shape.
Right:
Asymmetric emission line profile. We show the relative con-tribution of the intrinsic (Gaussian) emission lines measuredfrom the target pixel (‘self-contribution’) and the neighborpixels (colored-dashed lines) produced by beam-smearing. the opposite side of the Doppler shift direction in thespectral domain.REFERENCES