The rocky road to quiescence: compaction and quenching of quasar host galaxies at z ∼ 2

ABSTRACT

1 INTRODUCTION

In the last two decades, large-area submillimetre surveys have revealed a population of high-redshift galaxies with extreme levels of star formation, which were largely undetected in optical surveys as their ultraviolet (UV) emission is obscured by dust (see Blain et al. 2002; Casey, Narayanan & Cooray 2014; Hodge & da Cunha 2020, for reviews). These dusty star-forming galaxies (DSFGs) are expected to be precursors to locally observed massive elliptical galaxies, which are characterized by dense, old stellar populations with dispersion-dominated dynamics (Hopkins et al. 2008). A key aspect of the study of galaxy evolution is understanding how these galaxies formed such high stellar densities and grew concurrently with their central supermassive black holes (Magorrian et al. 1998). Recent near-infrared surveys have revealed that the population of compact quiescent galaxies start to appear at z ∼ 2 and rapidly increase in number density before apparently declining at z ∼ 1 (Trujillo et al. 2006; van Dokkum et al. 2008, 2015). These galaxies are around four times smaller in size than z ∼ 0 massive ellipticals and are thought to form the centre of these galaxies, which later grow in size (but little in mass) due to a series of gas-poor, minor mergers (Naab et al. 2007; Naab, Johansson & Ostriker 2009).

The characteristics and rapid formation of compact quiescent galaxies can be reproduced if a very high density of gas is concentrated within a region of ∼1 kpc to generate a brief, intense starburst (Valentino et al. 2020). However, the mechanisms that cause such rapid morphological change and quenching of star formation are currently unclear. Compaction may be due to a net loss of angular momentum that results from dynamical instabilities induced by gas-rich mergers or tidal interactions (Mihos & Hernquist 1996; Hopkins et al. 2008), or by non-axisymmetric structures caused by rapid accretion and clumpy star formation (Dekel et al. 2009; Dekel & Burkert 2014; Zolotov et al. 2015). Alternatively, these galaxies could have formed secularly, in situ at earlier epochs without the need for mergers or rapid evolution (Damjanov et al. 2011; Carollo et al. 2013; Williams et al. 2014; Wellons et al. 2015).

Star formation terminates as a result of depletion or cessation of the supply of molecular gas. In compact quiescent galaxies, this may happen as a result of compaction, where star formation is self-quenched by radiation pressure from massive stars and supernovae-driven winds (Murray, Quataert & Thompson 2005; Andrews & Thompson 2011; Diamond-Stanic et al. 2012). Eddington-limited ‘maximum’ starbursts have been discovered in DSFGs at high redshift (Riechers et al. 2013; Oteo et al. 2016; Cañameras et al. 2017; Spilker et al. 2019) as well as in nearby ultraluminous infrared galaxies (ULIRGs; Barcos-Muñoz et al. 2017). These can be sufficiently vigorous to drive large-scale outflows and quench their hosts (Cañameras et al. 2017; Spilker et al. 2018), from the inside–out, on time-scales of tens Myr (Spilker et al. 2019).

Hydrodynamical simulations and semi-analytical models of galaxy formation find that feedback from active galactic nuclei (AGNs) is necessary, in addition to stellar feedback, to reproduce observed galaxy stellar populations and luminosity functions (Di Matteo, Springel & Hernquist 2005; Sijacki et al. 2007; Somerville et al. 2008; Schaye et al. 2015). AGN feedback occurs in the form of jets or radiative winds, which can suppress star formation by mechanically coupling to the molecular gas and/or by preventing (re)accretion from the circumgalactic medium (see Fabian 2012, for review).

Due to computational limitations, current cosmological simulations involve only phenomenological implementations of feedback, which are calibrated to match observational data (see Somerville & Davé 2015). Observations find circumstantial evidence that AGNs could play a role in the evolution of their hosts, in that the compaction phase seems to coincide with the appearance of an AGN (Barro et al. 2013; Kocevski et al. 2017). However, studies have found conflicting results as to whether AGNs have any effect on their host galaxies. While some earlier studies found evidence of suppressed star formation in AGN hosts (e.g. Page et al. 2012), more recent studies of statistical samples have found no correlation between star formation and black hole accretion (Harrison et al. 2012; Rosario et al. 2013; Harris et al. 2016; Pitchford et al. 2016; Stanley et al. 2017; Kirkpatrick et al. 2019; Schulze et al. 2019). This could be because the effects of AGN feedback can only be detected later, or because black hole accretion is stochastic (Gabor & Bournaud 2014; Hickox et al. 2014). Even in simulations, star formation rates (SFRs) averaged over ∼100 Myr can show no clear correlation with black hole accretion (Harrison 2017), suggesting the effect of AGNs may not be obvious in studies of the global properties of quasar hosts.

These key tests of the evolutionary sequence require investigations of the size, structure, and dynamical properties of individual quasar hosts during the cosmic peak of galaxy growth (z ∼ 2). This demands high spatial resolution (hundreds of pc), which is most efficiently achieved by observing objects that are gravitationally lensed (e.g. Swinbank et al. 2010; ALMA Partnership 2015). By observing galaxies that are strongly gravitationally lensed, it is possible to recover the properties of the source with greater angular resolution and sensitivity. Stacey et al. (2018) conducted a survey of gravitationally lensed quasar systems with the Herschel Space Observatory to measure the level of obscured star formation in the quasar host galaxies. This paper presents observations with the Atacama Large (sub)Millimetre Array (ALMA) of a subsample of seven optically luminous quasars from the parent sample of Stacey et al. to resolve their host galaxy emission. Using a pixellated lens modelling technique applied to the interferometric data, we reconstruct dust and gas in the host galaxies and derive their intrinsic properties. In Section 2, we describe the targets, observations, and data reduction process. In Section 3, we describe our lens modelling and source reconstruction technique. Section 4 reports the results of the structure of dust and gas, gas dynamics, and star formation properties of the individual objects. We compare the morphological and star formation properties of the quasar host galaxies with a sample of DSFGs. In Section 5, we discuss the implications of our results in the context of evolutionary models and possible avenues to test our conjectures. Section 6 presents a summary of our findings and avenues for future work.

Throughout, we assume the Planck Collaboration XIII (2016) instance of a flat Lambda cold dark matter cosmology with H₀ = 67.8 km s⁻¹ Mpc⁻¹, Ω_m = 0.31, and Ω_Λ = 0.69.

2 SAMPLE AND OBSERVATIONS

In this section, we first summarize the sample, the observations with ALMA, and the data reduction processes, before giving a detailed description of the properties of each target studied here.

2.1 Summary of the sample, observations, and data reduction

The targets in this work are seven four-image gravitationally lensed quasar (Type 1 AGN) systems. Four targets, SDSS J0924+0219, PG 1115+080, WFI J2033−4723, and WFI J2026−4536, are from our own programme (PI: McKean) and were not selected on the basis of their FIR/sub-mm properties, but for the configuration of their lensed quasar images in optical data for reasons that are not relevant to this work. We also include here archival data of HS 0810+2554 (PI: Chartas), RX J0911+0551 (PI: Leung), and H1413+117 (PI: van der Werf), which are known to have bright sub-mm emission (Barvainis & Ivison 2002). Details of the observations are given in Table 1. From previous modelling of their spectral energy distributions (SEDs; Stacey et al. 2018; see Fig. A1), we expect the sub-mm emission to be entirely dominated by thermal dust. As these latter systems were primarily selected in optical imaging and subsequently selected on the basis of their far-infrared (FIR) luminosity, their selection is largely insensitive to dust temperature. We compile all accessible data with sufficiently high angular resolution to resolve the size and structure of the emission from cold dust at sub-mm/mm wavelengths. As a result, this sample is heterogeneous, and the observations probe different CO line transitions (see Table 2).

The raw data were calibrated using the ALMA pipeline in the Common Astronomy Software Applications package (casa; McMullin et al. 2007) to produce calibrated visibilities. The data were inspected to confirm the quality of the pipeline calibration and determine whether further flagging was required. RX J0911+0551, PG 1115+080, H1413+117, and WFI J2026−4536 were self-calibrated using the line-free spectral windows with solution intervals of each scan length. Self-calibration was attempted for SDSS J0924+0219 and WFI J2033−4723, but was not successful, probably because of the lower signal-to-noise ratio of the surface brightness. We did not attempt to self-calibrate the data for HS 0810+2554 as the continuum emission is extremely weak (see below).

The targets were imaged with natural weighting of the visibilities and deconvolved using clean (Högbom 1974). The deconvolved images are shown in Figs 2 and 3. The continuum emission is consistent with the level of thermal dust emission expected from SED fitting (Stacey et al. 2018). We detect the host galaxies of the quasars that, in some cases, form Einstein rings and gravitational arcs. We also find compact dust emission that, in all cases, approximately coincides with the quasar positions seen in optical data. No emission is detected from any of the lensing galaxies.

All data sets, except in the case of PG 1115+080, included observations of an emission line of CO. The spectral line data were prepared by fitting a linear model for the continuum to the line-free spectral windows and subtracting this from the visibilities. This produced a visibility data set that included only the line emission. The resulting line profiles, fit with single Gaussians, are shown in Fig. 1. In some cases, the line peaks are offset from the systemic velocity (the rest frame of the galaxy), which is possibly due to an uncertainty in the inferred redshift from optical spectroscopy: redshifts of quasars determined from optical spectroscopy can trace outflows of ionized gas, leading to incorrect inference on the systemic velocity. A summary of the spectral line observations is presented in Table 2.

Figure 1.

Line profiles for the six objects with CO observations. The red line shows the data; the blue dotted line is a Gaussian fit to the data. The systemic velocity is relative to the best redshift from the literature, using the radio definition of velocity.

Open in new tabDownload slide

The spectral line moment maps (velocity-integrated flux density, velocity field, velocity dispersion) are shown in Figs 2(a)–(f). The velocity field and velocity dispersion are made by masking channel pixels below a signal-to-noise ratio threshold (stated within the figure captions).

Figure 2.

Deconvolved ALMA images for the six systems with spectral line data, using natural weighting of the visibilities. Left to right: the continuum; moment 0 (velocity-integrated line intensity); moment 1 (line velocity); moment 2 (velocity dispersion). The beam FWHM is shown in the lower left corner.

Open in new tabDownload slide

2.2 HS 0810+2554

HS 0810+2554 is a quasar at z_s = 1.51 that is lensed into a characteristic fold configuration (Reimers et al. 2002; Hewett & Wild 2010). The redshift of the foreground galaxy is unknown, but is estimated to be 0.9 based on the lens population distribution (Mosquera & Kochanek 2011). This quasar is faint at radio wavelengths, but VLBI investigations find evidence of compact, low-luminosity radio jet structure that dominates the radio emission (Hartley et al. 2019). HS 0810+2554 was observed with ALMA at 145 GHz, which also targeted the CO (3–2) emission line.

We resolve thermal dust emission from this system with ALMA at 145 GHz that has a peak surface brightness of 6σ (where σ is the rms noise per beam). The achieved rms noise level is 13 μJy beam⁻¹ with natural weighting. Emission is detected around the close images, where the magnification is high and there is a hint of emission from the counter images. The total flux density of 0.3 mJy is consistent with the level of thermal dust emission expected at the observing frequency, based on SED fitting (see Fig. A1). The weak synchrotron emission seen at cm wavelengths is not expected to be detectable at the higher frequencies investigated here.

The CO (3–2) line profile for HS 0810+2554 is shown in Fig. 1, which has a full width at half-maximum (FWHM) of about 380 ± 10  km s⁻¹ based on a single Gaussian fit. We find that the peak of the line emission is shifted −465 ± 6  km s⁻¹ from the assumed optical systemic velocity, likely due to an uncertain redshift estimate. ALMA imaging of the line emission shows that the molecular gas is extended and lensed into an Einstein ring (see Fig. 2a).

2.3 RX J0911+0551

RX J0911+0551 (RX J0911.4+0551) is a quasar at z_s = 2.79 that is lensed by a foreground galaxy at z_l = 0.77 (Burud et al. 1998; Kneib, Cohen & Hjorth 2000). The lens system has a characteristic cusp configuration with three close images. The environment of this lens system is quite complex: the primary lens has a satellite galaxy within the Einstein radius, and the lens galaxy is part of a cluster that contributes a high level of tidal shear (Kneib et al. 2000). Investigations at radio wavelengths have not been able to determine whether there is radio jet emission (Jackson et al. 2015). However, its radio luminosity is consistent with the expectations for star formation based on the radio–infrared correlation (Stacey et al. 2018). For our analysis, we use ALMA imaging at 145 GHz that also targeted the CO (5–4) emission line from RX J0911+0551.

We detect thermal dust continuum for this system with a total flux density at 145 GHz of 1.4 mJy, consistent with the expectations from SED fitting (see Fig. A1). We achieve an rms noise of 7 μJy beam⁻¹ with natural weighting of the visibilities. Imaging of the dust shows compact emission that is resolved around the triplet images (see Fig. 2b). From the lens configuration, it is clear that the galaxy is crossing the cusp of the lens caustic.

The CO (5–4) line profile for RX J0911+0551 is presented in Fig. 1. The spatially integrated line profile has an FWHM of 133 ± 3  km s⁻¹  based on a single Gaussian fit, which is similar to the 120 ± 14  km s⁻¹  reported for the CO (7–6) line by Tuan-Anh et al. (2017). Imaging of the CO (5–4) line emission shows compact structure that is similar to the dust continuum (see Fig. 2b). The imaging is visually similar to the CO (11–10), CO (10–9), CO (7–6), and CO (1–0) emission reported by Tuan-Anh et al. (2017) and Sharon et al. (2016): none of which are extended into rings or arcs, suggesting the emission is similarly compact.

2.4 SDSS J0924+0219

SDSS J0924+0219 is a radio-quiet quasar at z_s = 1.525 that is lensed by a galaxy at z_l = 0.39 (Inada et al. 2003; Hewett & Wild 2010). This lens system was discovered at optical wavelengths and has a characteristic fold configuration. Observations with ALMA at 358 GHz, which also targeted the CO (8–7) line, were taken as part of our own observing programme.

The continuum emission from SDSS J0924+0219 is detected with a total flux density of 9.5 mJy at 358 GHz, consistent with the expectations for thermal dust emission from SED fitting (see Fig. A1). We achieve an rms noise of 28 μJy beam⁻¹ with natural weighting of the visibilities. The continuum imaging shows an Einstein ring of clumpy dust with compact emission at the locations of the quasar images. The extent of the dust seems consistent with a ‘red’ ring seen in optical/infrared imaging (Eigenbrod et al. 2006). The same study also finds lensed arcs from a ‘blue’ component, of which there is a suggestion in the ALMA imaging where the Southern arc appears to split. These different optical components could be evidence of an ongoing merger or interaction (e.g. Rybak et al. 2015b).

The CO (8–7) line profile is presented in Fig. 1. The line profile is found from a Gaussian fit to have an FWHM of 176 ± 11  km s⁻¹  that is shifted by 36 ± 5  km s⁻¹  from the optical systemic velocity. Imaging of the CO (see Fig. 2c) shows resolved emission at the position of the quasar images and evidence for a velocity gradient in the two close images. The velocity gradient is similar to that reported by Badole et al. (2020) for CO (5–4), where the CO emission is more extended. The line profile has a similar FWHM, but the shape of the CO (8–7) is more peaked than for the CO (5–4) and the surface brightness sensitivity is lower.

2.5 PG 1115+080

PG 1115+080 is a quasar at z_s = 1.74 that is lensed by a galaxy at z_l = 0.31 (Weymann et al. 1980; Hewett & Wild 2010). Investigations at radio wavelengths have not been able to determine whether the quasar has radio jet emission (Jackson et al. 2015), however, its radio luminosity is consistent with the expectations for star formation based on the radio–infrared correlation (Stacey et al. 2018). PG 1115+080 was targeted as part of our own programme with ALMA at 346 GHz.

We detect continuum emission with a total flux density of 2.6 mJy, consistent with the expectations for thermal dust emission from SED fitting (see Fig. A1). The imaging of the dust emission for PG 1115+080 is presented in Fig. 3 and shows compact structure at the position of the quasar images, without any evidence of an extended Einstein ring. Given that the total flux density is close to what is expected from this system, we do not believe that significant emission from an extended component has been resolved out with these data. However, the emission is marginally radially resolved around the two close images. Unlike in the cases of the other quasar hosts in this sample, the unobscured star formation seen in the optical and near-infrared appears more extended than in the obscured star formation we observe here (Sluse et al. 2012; Chen et al. 2019).

Figure 3.

ALMA 346 GHz continuum images of PG 1115+080 with natural weighting of the visibility data. The synthesized beam FWHM is shown in the bottom left-hand corner.

Open in new tabDownload slide

2.6 H1413+117

H1413+117 (the Cloverleaf) is a quasar at z_s = 2.56 that is lensed by a galaxy at unknown redshift (Magain et al. 1988; Riechers et al. 2011b). Due to its extreme luminosity at FIR wavelengths, the lens system is one of the best-studied quasar starbursts. While it is not extremely radio luminous, this system does have a radio excess relative to its infrared luminosity (Stacey et al. 2018) and has evidence of structure at radio wavelengths that may be from a low-luminosity radio jet (Stacey et al., in preparation). We use for our analysis archival data for H1413+117 that was taken at 285 GHz with ALMA, which also targeted the CO (9–8) emission line.

Thermal dust emission is detected in the ALMA imaging with a total flux density of 26 mJy and an rms noise level of 93 μJy beam⁻¹, using a natural weighting of the visibilities. The CO (9–8) line emission is detected with a spatially integrated peak flux density of 83 mJy and an FWHM of 394 ± 8  km s⁻¹, from a Gaussian profile fit (see Fig. 1); this is the brightest line emission detected from our sample.

We resolve the four images of the quasar host galaxy in continuum and the CO (9–8) molecular line emission, with high surface brightness peaks at the position of the quasar images and an extended, almost complete Einstein ring that connects them (see Fig. 2d). The structure of the emission appears similar to previous observations of dust and gas in this system (Alloin et al. 1997; Ferkinhoff et al. 2015). We also identify an additional source with a flux density of ≃ 3.5 mJy within the field at a distance ∼6 arcsec North of the lens. This is a potential source of confusion for the SED fitting, where the photometric measurements are typically made at a low angular resolution (e.g. IRAS or Herschel/SPIRE). This additional source may also explain the scatter in the photometric measurements from H1413+117 (e.g. see Stacey et al. 2018; Fig. A1). From a cursory inspection of the available ALMA archival data at a frequency of ∼100 GHz (project code 2015.1.01309.S), we find that the additional source is detected with a similar flux density ratio as at 290 GHz with respect to H1413+117, suggesting that these galaxies have similar redshifts.

2.7 WFI J2026−4536

WFI J2026−4536 is a quasar at z = 2.22 (Morgan et al. 2004) that is lensed by a galaxy at unknown redshift into four images in a typical fold configuration. There are no radio observations of this system available to characterize its radio properties at 1.4 GHz, but the lack of detection in the Sydney University Molonglo Sky Survey (SUMSS) Source Catalog (Mauch et al. 2003) suggests a radio flux density of <12 mJy at 1.4 GHz, assuming a typical synchrotron spectral index of −0.7. WFI J2026−4536 was observed with ALMA at 350 GHz as part of our ALMA programme, which also observed the CO (10–9) spectral line transition.

The ALMA observations detect continuum emission with a flux density of 28.3 mJy and achieve an rms noise level of 30 μJy beam⁻¹ with natural weighting. This flux density is in good agreement with thermal dust emission, based on SED fitting (Fig. A1). The clean ed image shows bright emission from the location of the quasar images and extended emission at lower surface brightness.

The CO (10–9) line emission is also spatially resolved with a similar morphology to the dust emission. The spatially integrated CO line profile has a peak flux density of 49 mJy and shows a hint of a double-horned feature characteristic of disc kinematics. The FWHM of the line is 252 ± 7  km s⁻¹  and is offset by ∼1000  km s⁻¹  from the systemic velocity assumed from optical spectroscopy.

2.8 WFI J2033−4723

WFI J2033−4723 is a quasar at z_s = 1.66 that is lensed by a galaxy at z_l = 0.66 (Morgan et al. 2004) into four images in a typical fold configuration. The lack of detection in the SUMSS suggests a radio flux density of <12 mJy at 1.4 GHz. No other radio observations are available for this lens system. The target was observed with ALMA at 341 GHz as part of our programme, which also observed the CO (8–7) emission line from the lensed quasar.

The ALMA observations detect thermal dust emission with a total flux density of 10.4 mJy and achieve an rms noise of 30 μJy beam⁻¹ with natural weighting of the visibilities. Imaging of WFI J2033−4723 shows a ring of extended, clumpy dust emission and compact emission from the approximate location of the quasar images (see Fig. 2f). Arcs of faint emission from the host galaxy have also been observed at optical wavelengths, similar to the extended dust emission observed here (Rusu et al. 2020). The observed flux density at 341 GHz is in good agreement with the expectations of SED fitting to a thermal dust model (see Fig. A1).

We detect the CO (8–7) line emission with a spatially integrated peak flux density of 11 mJy. The FWHM of the line profile is quite narrow, at 83 ± 10  km s⁻¹, and is offset by 258 ± 4  km s⁻¹  from the optical systemic velocity (see Fig. 1). The line profile shows a broad blueshifted component to the main peak, but it is not clear whether this is a real feature or due to the low signal-to-noise ratio of the data. Imaging of the CO line shows resolved, compact emission that is coincident with the quasar positions (see Fig. 2f).

3 LENS MODELLING

Our lens modelling analysis is adapted from the grid-based modelling technique of Vegetti & Koopmans (2009). The technique simultaneously optimizes for the parameters of the lens galaxy mass distribution and the source surface brightness to produce a pixellated reconstruction of the source. This method is appropriate for optical or infrared images. However, interferometers do not measure the sky surface brightness distribution directly, but measure complex ‘visibilities’ in the Fourier (uv) plane. Imaging interferometric data is a non-linear process, and the resulting image fidelity can depend heavily on the sparsity of the uv-plane coverage and choice of deconvolution method (e.g. clean; Högbom 1974).

These systematic uncertainties can be overcome by performing model fitting on the visibility data directly. This approach has been extended to lens modelling, with parametric source models (e.g. Bussmann et al. 2013, 2015; Hezaveh et al. 2013) and also pixellated sources (e.g. Rybak et al. 2015a,b; Hezaveh et al. 2016; Dye et al. 2018). However, these frameworks are limited by the size of the data and require averaging in time and/or frequency in order to be computationally tractable. Recently, Powell et al. (2020) introduced several algorithmic improvements to overcome this limitation and model large interferometric data sets directly. We employ this method for optimizing the visibility plane posterior probability and source plane regularization.

Under the assumption of a smooth mass distribution, we parametrize the primary lens potential as a singular power-law ellipsoid with external shear, i.e. ρ(r)∝r^−ζ (Kormann, Schneider & Bartelmann 1994). For all objects, with the exception of HS 0810+2554, we perform the model optimization on the continuum data. For HS 0810+2554, we perform the lens modelling with the integrated CO (3–2) spectral line, which provides more extended emission and a higher signal-to-noise ratio than the continuum.

The lens redshifts are not known for HS 0810+2554, H1413+117, or WFI J2026−4536. However, the lens redshift only scales the angular size and mass of the perturber, so the inferred properties of the source are not affected.

Robust lens models have been produced for PG 1115+080 and WFI J2033−4723, based on higher quality optical imaging of the quasars and extended arcs as part of the H₀ Lenses in COSMOGRAIL’s Wellspring (H0LiCOW) project (Chen et al. 2019; Rusu et al. 2020). Therefore, in these cases we keep the lens slope, ellipticity, and position angle fixed to literature values, and reoptimize only the lens galaxy position (which is not observed in our data) and shear parameters (for which we do not consider second-order corrections). In all other cases, we assume isothermal mass distributions (i.e. ζ = 2) and optimize for the lens galaxy position, ellipticity parameters, and external shear parameters.

For WFI J2033−4723, our lens model also includes a luminous satellite galaxy close to the Einstein radius of the primary lens and a massive galaxy 4 arcsec from the primary lens, both of which we parametrize as a single isothermal sphere (SIS). We find the source size changes by ≥10 per cent with the addition of these galaxies, so it is necessary to explicitly include them in our model. The position and mass of these galaxies are fixed to the values found by Rusu et al. (2020). The lensing galaxy of RX J0911+0551 also has a small satellite, however, this is well within the Einstein radius and we find it does not have a significant effect on our lens modelling.

The uncertainties on the lens model parameters were obtained from the posterior distributions derived using multinest (Feroz et al. 2013), with the source regularization constant fixed to the maximum a posteriori value. Flat prior ranges were assumed for each free parameter, typically ±20 per cent of the optimized values. For WFI J2026−4536, it was necessary to keep the lens position fixed to the maximum a posteriori (consistent within 20 mas of the optical position) for stability during the nested sampling. The likelihood-weighted posterior probability distributions for the lens parameters are given in Table 3.

Reconstructed velocity cubes were generated from the spectral line data by optimizing for the source regularization, using the best-fitting lens model. For WFI J2033−4723 and SDSS J0924+0219, we use only three velocity channels and restrict the uv-data to within a maximum baseline of 300 m to improve the surface brightness sensitivity.

4 RESULTS

In this section, we present the results of the source reconstruction and an analysis of their physical properties.

4.1 Source reconstructions

The maximum a posteriori lens models are shown in Figs A2 and A3 of the Appendix. The residual maps are a dirty image produced after the model has been subtracted from the visibility data, normalized to the rms noise in the data. We find our lens models fit the data well and the residual surface brightness features are below the 4σ level.

The reconstructed sources and moment maps are shown in Figs.   4 and 5. For SDSS J0924+0219 and WFI J2033−4723, the dust emission is extended and resolved into clumpy features on kpc-scales. For HS 0810+2554, RX J0911+0551, PG 1115+080, H1413+117, and WFI J2026−4536, the reconstructed dust emission appears to be very compact and does not show evidence for clumpy features on the same scales. We caution that, for HS 0810+2554, there may be extended dust emission that is not detected here due to the low signal-to-noise ratio of the data. The apparent smoothness of the compact systems may be due to angular resolution limitations, however, we note that other unlensed dusty galaxies at z ∼ 2 are found to have smooth dust distributions (e.g. Falgarone et al. 2017). The radially averaged surface brightness distributions of the reconstructed continuum sources are shown in Fig. 6.

Figure 4.

Reconstructed dust emission and moment maps for the six lensed quasar systems with spectral line data. The dust emission is shown in units of Jy arcsec⁻² with the lensing caustics shown in grey. Dashed contours show signal-to-noise ratio contours, starting at 3. Continuum contours are shown in grey over the moment 0 image. The bar shows 1 kpc at the redshift of the source. The line reconstructions are generated by masking pixels in each channel below 3σ. The CO velocity is corrected to the redshift derived from a Gaussian fit to the CO line profile.

Open in new tabDownload slide

Figure 5.

Reconstructed dust emission of PG 1115+080, shown in units of Jy arcsec⁻² with the lensing caustics shown in grey. The white bar shows 1 kpc at the redshift of the source.

Open in new tabDownload slide

Figure 6.

Normalized, azimuthally averaged surface brightness profiles of the reconstructed dust emission from the quasar hosts, accounting for ellipticity. The shaded regions show the respective standard deviation of the dust emission. The radius of HS 0810+2554 and PG 1115+080 is scaled to the size of the minor axis, based on the axial ratio from the Sérsic model fits (see Section 4.1). The black dotted line shows the mean profile of (unlensed) ALESS DSFGs from Sérsic fits (Hodge et al. 2016, 2019).

Open in new tabDownload slide

The robustness of the reconstructed surface brightness distribution of the lensed source depends on the quality of the data. As our source is non-parametric, our fitting down to the level of the noise can allow correlated noise features in the dirty image to be absorbed into the source reconstruction. Our investigations find that the uncertainties in the source surface brightness resulting from noise artefacts are significantly larger than the uncertainties due to the lens model parameters (Rizzo et al., in preparation). To estimate these uncertainties, we assume our maximum a posteriori source model and create mock data with 100 different realizations of the noise at the level measured in the real data.¹ We create reconstructions of these mock data sets and measure the mean and standard deviation of the surface brightness in each pixel (in each channel, for spectral line data). This allows us to discriminate features in the reconstructed source from noise artefacts.

To estimate the uncertainty on the source size, we fit a Sérsic profile (Sérsic 1963) to the source from each noise realization of continuum and integrated line emission using a basin-hopping optimization algorithm within the scipy package (Wales & Doye 1998; Jones et al. 2001). For the integrated line emission of SDSS J0924+0219 and WFI J2033−4723, we fix the index of the Sérsic profile to that found for the continuum as the fitting was found to be unstable due to the low signal-to-noise ratio of the data. We assume the standard deviation of the inferred Sérsic parameters as their uncertainty. We estimate the magnification of each lens system by generating a model for the lensed emission for each noise realization, where source pixels are masked below a signal-to-noise ratio of 4. This value is chosen so as not to include any noise features that may be present in the individual realizations that may bias the magnification to lower values. We take the mean and standard deviation of the magnifications as the inferred value and its uncertainty (given in Tables 4 and 5).

In order to test the robustness of our inferred source sizes, we created mock images of Gaussian sources lensed by the maximum a priori lens model for each object. The Gaussian sources were located at the position of peak surface brightness of the maximum a priori source and given a range of effective radii. Mock visibilities were created with the same uv-coverage, baseline-based rms noise, and image peak surface brightness as the real observations. We modelled and analysed these mock data sets with the same methodology as the real data, optimizing for the source regularization parameter. We found that the inferred source sizes were inflated by a constant value that can be related to the strength of the source regularization. This inflation has little effect where the source is extended on kpc-scales, but becomes significant where the source is very compact. We use these simulations to infer a correction to the sizes derived from the Sérsic fits.

From our simulations, we confirmed that we can recover the source size from the real observations. However, for PG 1115+080 and HS 0810+2554 continuum, we determine that the lensed emission is only marginally resolved in the tangential direction. This results in a source that is constrained only by the angular resolution of the data in one dimension. Therefore, for these two cases, we assume the minor axis of the fitted Sérsic profile as the source size. The inferred magnifications from the real data are generally consistent with those inferred from the mock data sets, allowing for minor differences due to source structure. However, the inferred magnifications of PG 1115+080 and HS 0810+2554 (continuum) were underestimated due to the effect of source inflation, so we infer a correction using the mock data sets.

A consequence of the lack of radial resolution is that there is a degeneracy between the normalized convergence (κ₀) and the slope of the lens mass density profile (ζ). The effect of this degeneracy is a geometric scaling of the reconstructed source. In most cases, we assume an isothermal profile for our lens models (i.e. ζ ≡ 2), however, the empirical scatter in the total mass slope of early-type galaxies (0.16; Auger et al. 2010) implies an additional uncertainty in the source size. We tested this by modelling PG 1115+080 with a profile fixed to a value ±1σ from the median found by Auger et al. (2010). We find that the size and normalization of a Sérsic fit to the reconstructed source changes by ≤10 per cent. Therefore, where we assume only an isothermal profile is fit to the data, we propagate an additional ±10 per cent error into the uncertainty in our estimated sizes to account for the empirical scatter in ζ.

All the derived parameters of the continuum and line emission for each source can be found in Tables 4 and 5, respectively. Note that the given magnifications are mean magnifications, as the magnification changes across the source and between velocity channels.

4.2 Dust properties

We apply a Markov chain Monte Carlo (MCMC) analysis using the python implementation emcee (Foreman-Mackey et al. 2013) to infer the marginalized posterior distributions for each free parameter and the integrated FIR luminosity. The SED models are shown in Fig. A1 of the Appendix. The parameters from the SED modelling, corrected for the lensing magnification, are listed in Table 4.

The dust temperature we infer from our SED fitting is an effective dust temperature, which should not be considered a true measure of dust temperature, but weighted by multiple components of dust emission. This may include a cold dust component from cirrus dust of the interstellar medium (ISM; ∼20 K) and a warmer dust component associated with star-forming regions (∼50 K). Indeed, Swinbank et al. (2014) find DSFGs require at least three dust components to account for the full mid-infrared (MIR) to FIR SED. AGN host galaxies may also have contributions from AGN-heated dust emission. For HS 0810+2554 and H1413+117, we account for a warm dust component, which we assume is associated with the AGN. For HS 0810+2554, β is kept fixed in the SED modelling as we do not have sufficient data points to fit both β and two dust components. We fix β to a value of 2, as 1.5 cannot produce a satisfactory fit to the data.² With the exception of H1413+117 and HS 0810+2554, we do not have data at rest-frame wavelengths ∼10 to 100 μm to consider additional dust emission components.

Note that, while we leave β as a free parameter, its value is not necessarily physically significant as, for poorly sampled SEDs, the strong temperature–β degeneracy absorbs the observational noise in the data (Juvela & Ysard 2012). We allow this as a free parameter to better estimate the true uncertainty on the fitted temperature and luminosity.

It has often been found that a dust mass estimated from single-temperature SED model may underestimate the total dust content, as most of the mass is in cold cirrus dust and not the dust that contributes most to the effective (luminosity-weighted) temperature (e.g. Scoville et al. 2014). However, there is also evidence that multicomponent dust models could lead to overly large dust masses (Cortzen et al. 2020) and it has not been established to what extent the approach is appropriate for quasar host galaxies. Note that, if the temperature of dust in the diffuse ISM is much lower than the effective dust temperature for quasar hosts, the dust masses derived for our sample could be underestimated. We attempted fitting two-component models to our SEDs but, as most are sparsely sampled, the uncertainties on the mass of the cold component was very large and encompasses the masses derived with a single-temperature fit. Thus, we report the dust masses from a single-temperature fit, but with this caveat.

The dust temperatures, emissivities, intrinsic luminosities, SFRs, SFR surface densities, and dust masses for the sample are given in Table 4. For H1413+117, the nearby field source we identify in our ALMA imaging (Section 2.6) may cause modest confusion in SED fitting, which leads to overestimates of the inferred star formation properties of the source of interest. Photometric measurements at multiple frequencies will be required to accurately deblend this emission from the total infrared luminosity. Such an analysis is beyond the scope of this work, however, the consistent relative flux density of the objects at 100 and 290 GHz suggests that the neighbour contributes ∼20 per cent to the photometry of H1413+117. Therefore, our measurements based on infrared luminosity are likely to be overestimated by a similar percentage. We estimate the effect of this by scaling all the photometric measurements above 300 GHz by a factor of 0.8. The physical properties based on the rescaled data are also shown in Table 4.

4.3 Molecular gas properties

Fig. 7 shows the intrinsic effective radius of the dust and CO emission against effective dust temperature and total infrared luminosity (8–1000 μm). The sample in this work is quite heterogeneous and different rotational transitions of CO were observed. These CO lines correspond to different physical conditions in the ISM. At one extreme, CO (10–9), CO (9–8), and CO (8–7) trace warm gas (150–200 K) in regions of intense star formation or that are associated with AGNs. For the four objects where we observe one of these lines (SDSS J0924+0219, H1413+117, WFI J2026−4536, and WFI J2033−4723), the molecular gas is a similar size or more compact than the dust continuum. For HS 0810+2554 and RX J0911+0551, in which we observe CO (5–4) and CO (3–2), which are usually associated with star-forming regions, we find the gas to be more extended than the dust. Different relative sizes of the emission regions of CO line transitions have been found for other DSFGs (e.g. Apostolovski et al. 2019). While little can be inferred about the relative contributions to the energy budget from observations of a single line transition, the observed sizes are consistent with high levels of star formation and with the expectations for a radially decreasing gas column density and temperature (Weiß et al. 2007).

Figure 7.

The effective radius (R_eff) of the dust continuum (solid circles) and CO line emission (open triangles) against effective dust temperature (T_d) for the seven lensed quasar hosts investigated here.

Open in new tabDownload slide

We convert our line flux densities to luminosities (L_CO and |$L^{\prime }_{\rm CO}$|⁠) with relations as given by Solomon & Vanden Bout (2005). We use the source plane reconstructions to derive the intrinsic line properties to account for differential magnification across the velocity channels. As some line velocity components may be more strongly magnified than others, the integrated line emission measured in the lens plane is not necessarily the same as the intrinsic line emission multiplied by the mean magnification (i.e. |$I^{\rm lensed}_{\rm CO} \ne \bar{\mu }\ \times I_{\rm CO}$|⁠). The intrinsic line luminosities (corrected for lensing magnification) are given in Table 5. The reconstructed line emission in units of luminosity surface density, as well as the velocity field and velocity dispersion, are shown for each source in Figs 4(a)–(f).

In addition to the availability of gas, the level of star formation in a galaxy depends on the local balance between turbulent pressure and self-gravity. In dynamically unstable systems, self-gravity dominates over turbulent pressure allowing gas to fragment and efficiently form stars. As the systems in this work are known to have high SFRs and, in some cases, have clumpy non-axisymmetric dust structures, we expect the gas discs to be globally dynamically unstable. As shown in Figs 4 and 5, the reconstructed velocity maps appear to show radially increasing rotational velocity and high central dispersions. However, these two-dimensional maps may be strongly affected by beam smearing, which can misguide interpretation of the dynamics of the systems (velocity dispersion in particular, e.g. Lelli, Fraternali & Sancisi 2010). Also note that for SDSS J0924+0219 and WFI J2033−4723, we use only three velocity channels for the reconstructions. We defer a 3D analysis to future work, combining lens and source kinematic modelling as introduced by Rizzo et al. (2018; see also Rizzo et al. 2020).

4.4 Comparison with dusty star-forming galaxies

We make a comparison between our lensed quasar hosts and high angular resolution ALMA observations of DSFGs. These samples include galaxies from the LABOCA ECDFS Submm Survey (ALESS; Hodge et al. 2016), strongly lensed DSFGs selected with Herschel surveys, and several DSFGs at z > 4. The choice of these samples is motivated by the different selection biases that dominate in FIR/sub-mm surveys: we choose these samples to represent a diverse range of properties observed for DSFGs and combat the effects of selection bias.

DSFGs in the ALESS sample cover a similar redshift range (z ∼ 2, including one at z = 3.4) and infrared luminosities to our lensed quasar sample. We perform SED fitting for the ALESS sample with photometry from Swinbank et al. (2014) and Hodge et al. (2019). For consistency, we use the same methodology as for the quasar hosts to obtain effective dust temperatures and infrared luminosities with an optically thick dust model. Following this methodology, we only include the 11 ALESS galaxies where there are sufficient detections at FIR–sub-mm wavelengths to constrain the dust temperature.

We note that sub-mm source selection may preferentially identify sources with lower dust temperatures at a given redshift, due to the effect of negative k-correction (Chapin et al. 2011). This does not strongly influence the selection of the quasar hosts (which were mostly selected on the basis of their optical properties and lens configuration, as detailed in Section 2) but may influence the selection of the coeval DSFGs. To account for this selection bias, we include a sample of DSFGs at z > 4 and several lensed DSFGs. The z > 4 sample includes one of the ALESS DSFGs (Hodge et al. 2019) and four other objects that are spatially resolved, and have spectroscopic redshifts and multiple detections at FIR–sub-mm wavelengths. These include GN20 (Hodge et al. 2015), HFLS3 (Riechers et al. 2013), AzTEC-3 (Riechers et al. 2014), the two components of SGP38326 (Oteo et al. 2016), and the two components of ADFS-27 (Riechers et al. 2017). A summary of the DSFGs is given in Table A1 of the Appendix.

The lensed DSFGs were selected in FIR surveys with Herschel/SPIRE (Negrello et al. 2017). We perform the same analysis for the lensed DSFGs as for the lensed quasar sample; further details of the lens modelling and source reconstructions will be presented in a follow-up paper (Stacey et al., in preparation). The selection of lensed DSFGs is expected to be strongly biased towards compact systems, where the flux density is boosted by high magnification (Hezaveh, Marrone & Holder 2012; Serjeant 2012). We therefore expect these lensed systems to be examples of more compact DSFGs.

Fig. 8 (left) shows the continuum source size (effective radius) of the quasar hosts in this work and the z ∼ 2 DSFGs. This shows that more compact dust emission is associated with higher effective dust temperatures. A Kendall rank test (Kendall 1945) for correlation between dust size and temperature yields a coefficient of τ = −0.37 with a significance of p = 0.04 (where p = 0.05 is often taken as significant). Such a relationship is expected as a natural consequence of the Stefan–Boltzmann law, which relates the luminosity of a blackbody to its temperature and physical size. If the luminosity and temperature of the thermal dust emission we observe from the quasar hosts is related only to star formation, they should follow this relation.

Figure 8.

Left: Effective radius (R_eff) against effective dust temperature (T_d), coloured by SFR. Right: Effective radius against mean SFR surface density (Σ_SFR), where the grey region shows our estimated Eddington limit and its uncertainty. Σ_SFR and the dashed curves show the tracks for log L_IR values of 11.5, 12, 12.5, 13, and 13.5. Circles show the quasar hosts in this work (with H1413+117 corrected for SED blending); stars show coeval DSFGs at z ∼ 2 (Hodge et al. 2016, 2019); triangles show DSFGs at z > 4 (Riechers et al. 2013, 2014, 2017; Hodge et al. 2015, 2016). Four of the seven quasar hosts have SFR surface densities close to the Eddington limit, similar to DSFGs selected at z > 4.

Open in new tabDownload slide

Fig. 9 shows the intrinsic infrared luminosity (and equivalent SFR) against effective dust temperature for our sample and the DSFGs. The relationship expected for the modified Stefan–Boltzmann law from Yan & Ma (2016) is also shown for different source sizes, assuming a dust model with β = 1.5 and ν₀ = 3 THz. The DSFGs and quasar hosts can be seen to broadly follow this relationship with some scatter, as expected due to source structure as well as variations in dust emissivity and opacity. Fig. 9 shows that two lensed quasars have sizes in dust emission larger than implied by their dust temperature, which could indicate the sub-mm emission has a significant contribution from AGN heating. However, we note that for these systems (SDSS J0924+0219 and WFI J2033−4723), the dust distribution is resolved into multiple clumps and a single Sérsic component is likely a poor descriptor of their surface brightness distribution (Figs 4c and f). This is also the case for some lensed DSFGs. This blending effect was also suggested by Yan & Ma (2016), who found some DSFGs have sizes in dust emission larger than expected from the Stefan–Boltzmann relation.

Figure 9.

Effective dust temperature against infrared luminosity, coloured by effective radius. The dashed lines show the modified Stefan–Boltzmann relation from Yan & Ma (2016) for sources of indicated physical sizes, assuming a dust model with β = 1.5 and ν₀ = 3 THz. Sources from this work are shown with circles. DSFGs at z ∼ 2 are shown with stars; DSFGs at z > 4 are shown with diamonds; lensed DSFGs are shown with triangles. The quasar hosts in this sample have higher effective dust temperature and are generally more compact than dusty star-forming galaxies of similar infrared luminosity.

Open in new tabDownload slide

Fig. 9 suggests that the dust emission of quasar hosts generally has smaller sizes and higher temperatures than DSFGs of similar infrared luminosity, and smaller sizes and higher luminosities than DSFGs of similar temperature. The z ∼ 2 sub-mm-selected DSFGs are a factor of ∼3 larger than quasar hosts, the z > 4 DSFGs have intermediate sizes between the z ∼ 2 DSFGs and quasar hosts, and lensed DSFGs span a broad range of sizes in dust emission (0.5 < R_eff < 3.8) the smaller of which are similar to lensed quasars.

4.5 Intensity of star formation

HS 0810+2554 is the only object where we observe CO (3–2), which is a common proxy for molecular hydrogen (H₂; Greve et al. 2014). We assume α_CO = 0.8 M_⊙ (K km s⁻¹ pc²)⁻¹, as is often assumed for dusty starbursts, to convert |$L^{\prime }_{\rm CO}$| to total gas mass (M_gas). For this object, we find M_gas = (2.1 ± 0.1) × 10⁹  M_⊙  and Σ_gas = 4.5 ± 0.7 × 10⁹  M_⊙  kpc⁻² (derived in the same manner as equation 7). This implies a very low gas fraction of just 0.06 ± 0.04 compared to its dynamical mass, and a dust-to-gas ratio is 0.006 ± 0.002 (slightly below solar metallicity). The implied gas depletion time-scale (i.e. t_dep ≡ M_gas/SFR) is 50 ± 40 Myr.

5 DISCUSSION

5.1 Evidence for extreme star formation in quasar host galaxies

The seven objects investigated here were selected from a larger survey of 104 lensed quasars with Herschel/SPIRE that was reported by Stacey et al. (2018), who performed SED fitting to constrain the level of dust-obscured star formation in the sample. In all but three cases (including H1413+117), a single dust component was used for the fitting due to the lack of data in the MIR to FIR. From this analysis, Stacey et al. (2018) found a median dust temperature of |$38^{+12}_{-5}$| K (based on optically thin dust models) and a median SFR of |$120^{+160}_{-80}$| M_⊙ yr⁻¹, both of which are typical of DSFGs. A concern from such a simple SED model is that AGN-heated dust emission could contribute to the FIR and result in a higher dust temperature and, hence, a higher effective dust temperature (e.g. Kirkpatrick et al. 2012). Indeed, in some cases, effective dust temperatures were found that are higher than typical for star formation. However, we find that the quasar hosts are generally consistent with the Stefan–Boltzmann relation between temperature and luminosity when considering an optically thick dust model (see Fig. 7), suggesting that the dust is uniformly heated to higher temperatures by an obscured starburst, as opposed to the central AGN. This does not exclude the possibility that the AGN contributes to the cold dust emission, but in such a scenario it is unlikely be a significant fraction without a conspiracy between the effective radius, temperature, and luminosity. The higher dust temperatures found for some quasar hosts relative to DSFGs can be explained by the more compact physical size of their dust emission.

We estimate that five compact quasar hosts in our sample have high SFR surface densities which may be Eddington-limited. This finding is consistent with previous investigations of H1413+117, which have independently verified an Eddington-limited starburst with radiative transfer modelling of spectral line ratios (Bradford et al. 2009; Riechers et al. 2011a; Uzgil et al. 2016). In general, while an X-ray radiation field may penetrate large column densities of molecular gas, models suggest it will not efficiently heat dust grains over large (∼kpc) distances (Meijerink, Spaans & Israel 2007). These factors further suggest that the AGNs are unlikely to be responsible for a significant fraction of the sub-mm dust emission, but is more ambiguous for more compact sources with an effective radius of ∼100 pc and super-Eddington starbursts. We cannot determine the origin of the molecular gas heating from single CO lines, but a comparison between multiple transitions will help constrain the contribution from the AGN to the emission reported here. Radiative transfer models of quasar hosts find that these are very likely to be enhanced by X-rays from the AGN (Vallini et al. 2019), so could probe the effect of radiative feedback on the host galaxy ISM.

Another potential source of uncertainty in the analysis by Stacey et al. (2018) was the assumption on the unknown lensing magnification of the FIR emission, which was set to be |$10_{-5}^{+10}$| where the FIR magnification is unknown. If the lensing magnifications were significantly higher than was conservatively assumed, then the level of the inferred star formation would be much lower. From our analysis of the seven lensed quasars investigated here, we find dust magnifications factors in the range 11–25, roughly consistent with the general estimate of 5–20 by Stacey et al. (2018). While on the higher side of the estimate, these are four-image lens systems so high magnifications are typical. Accounting for the inferred magnifications, we find that the SFRs for the sample here range from 40 to 1600 M_⊙ yr⁻¹.

5.2 Evidence for compact quasar hosts

The quasar hosts in our study can be divided into two types: DSFGs characterized by clumpy dust distributions, with sizes and SFR densities similar to typical sub-mm-selected DSFGs (i.e. WFI J2033−4723 and SDSS J0924+0219), and DSFGs characterized by compact (R_eff < 1 kpc) sizes and generally high star formation densities (see Fig. 8) with no evidence of clumpy features on kpc-scales.⁴

The sizes of the compact quasar hosts in this work are similar to the sub-mm sizes of compact star-forming galaxies (Ikarashi et al. 2015; Barro et al. 2016, 2017). The lensed DSFGs (which we expect to represent more compact DSFGs due to their selection bias, e.g. Serjeant 2012) show a large range of sizes, yet the quasar hosts of similar luminosity are generally more compact and have higher dust temperatures. This may indicate that quasars are preferentially hosted in more compact systems. Even considering that two of these quasar hosts are likely to be selected on the basis of their FIR properties, the remaining optically selected systems are distinct from the remaining sample.

In this work, by selecting hosts of Type 1 quasars, we have explicitly chosen galaxies whose AGNs are rapidly accreting and generating relativistic winds that have exposed their accretion discs. It is predicted that a significant fraction of DSFGs host AGN (e.g. Hickox et al. 2014). While it has been found that at least one in five DSFGs have an X-ray luminous AGN (Wang et al. 2013), this is true for at least half of compact star-forming galaxies (Barro et al. 2014; Kocevski et al. 2017).⁵ The relative prevalence of X-ray luminous AGNs may point towards elevated black hole growth in compact DSFGs. The small size of compact galaxies requires that there has been a rapid infall of gas such that the bulk of the star formation occurs within the central region of the galaxy: this could allow the AGN to accrete more efficiently from the dense ISM.

The sub-mm sizes of both the compact DSFGs and quasar hosts are similar to the optical/infrared sizes of compact quiescent galaxies at similar redshifts (⁠|$\tilde{R}_{\rm eff}=0.9$| kpc; van Dokkum et al. 2008). Assuming the sub-mm dust emission traces the bulk of the stellar component, this is consistent with the hypothesis that compact starbursts are progenitors of compact quiescent galaxies (and ultimately massive ellipticals). The still-high SFRs and molecular gas densities of these hosts could mean that we catch them at the moment of compaction, at the onset of gas depletion and quenching. This would require that at least some compact quiescent galaxies form rapidly, rather than through a slow evolution (Wellons et al. 2015).

Overall, the sizes of < 1–3 kpc we find for our sample are similar to those of Silverman et al. (2019), who found the optical/infrared sizes of quasar hosts at z ∼ 1.5 are in the range <1–6 kpc (⁠|$\bar{R}_{\rm eff}=2.2$| kpc), intermediate between main sequence and quiescent galaxies. Additionally, Ikarashi et al. (2017) also found that z ∼ 1–3 quasar–starburst composites have smaller sub-mm dust sizes than both starburst-dominated DSFGs and quasars whose SEDs have a more dominant AGN fraction in the MIR. D’Amato et al. (2020) also found similar sizes for the dust components of X-ray-selected quasars at 2 < z < 5. These results and ours are inconsistent with models that predict the size of quasar hosts to be larger than those of normal star-forming galaxies, due to adiabatic expansion that results from negative AGN feedback (Fan et al. 2008) or positive AGN feedback on kpc-scales (Ishibashi, Fabian & Canning 2013).

5.3 Mechanism of formation

In the context of the ongoing debate about how compact quiescent galaxies form their very high stellar densities, our finding that the host galaxies of quasars are compact is in agreement with a model whereby they form rapidly in a period of dissipative contraction. Simulations suggest that star-forming galaxies are able to maintain high SFRs in quasi-stable discs, fed by smooth accretion from the cosmic web (e.g. Kereš et al. 2005). For a galaxy to compactify, dynamical instabilities must be induced by mergers or intense inflows that can drive gas into the centre of the galaxy on time-scales shorter than the SFR (Dekel & Burkert 2014; Zolotov et al. 2015). The high SFRs and the existence of clumpy, non-axisymmetric dust features that we detect strongly suggest that quasar hosts are globally unstable due to gravitational fragmentation, as has been frequently observed for DSFGs (Iono et al. 2016; Oteo et al. 2017; Hodge et al. 2019).

Optical imaging of DSFGs has often found that their molecular gas and obscured star formation is significantly more compact than pre-existing stellar distributions (Ikarashi et al. 2015; Simpson et al. 2015; Kaasinen et al. 2020). The optical/UV-luminous stellar emission observed for PG 1115+080 (see Peng et al. 2006) must also be more extended or offset from the dust, as it is lensed into an Einstein ring, whereas the dust is not. These unobscured stellar features may exist for the other quasar systems we observe here, but are too faint to be seen in optical imaging or cannot be distinguished from the bright quasar emission. This UV-luminous emission likely contributes only a small fraction of the star formation, which is largely obscured, but these features may hint at the formation histories of the galaxies. For example, the unobscured stellar emission may be quenched as a result of gas compaction and stellar feedback (Maiolino et al. 2015). However, this may be difficult to interpret due to the much higher sensitivity to star formation achieved by optical/UV imaging compared to sub-mm imaging.

5.4 Mechanism of quenching

Simulations find that compaction can naturally lead to quenching, through rapid gas consumption and stellar feedback, coupled with increasing dynamical stability (e.g. Zolotov et al. 2015). Observations of dense starbursts suggest that radiation pressure from stars could be a feasible mechanism to suppress further star formation (Murray et al. 2005; Andrews & Thompson 2011). If all the dust emission we observe here is associated with star formation, the five most compact sources in our work are in the regime of Eddington-limited maximum starbursts. This is consistent with observations of compact quiescent galaxies, whose star formation histories suggest a short starburst of ∼50 Myr before quiescence (Valentino et al. 2020). In the absence of reliable stellar mass estimates, we cannot determine whether star formation in the compact quasar hosts has begun to quench by direct comparison with DSFGs. However, the higher implied SFR surface densities imply that stellar feedback could play the primary role.

For the one system where we can trace the extent of the gas reservoir, we find a very low gas fraction in comparison to its dynamical mass. This suggests that this galaxy may be transitioning into quiescence following the depletion and removal of gas, as observed for compact star-forming galaxies at z ∼ 2 (Spilker et al. 2016, 2019). The implied depletion time-scale of 50 ± 40 Myr is comparable to the finding of ∼ 100 Myr for compact star-forming galaxies by Spilker et al. (2016) and consistent with ∼50 Myr found by Valentino et al. (2020).

The compact starbursts and implied short depletion time-scales could suggest that AGN feedback does not play an important role in the immediate quenching of star formation, contrary to simulations. However, long-term maintenance of quenching requires that the supply of fresh cold gas into the galaxy is halted. As the quasars in this work are likely hosted in massive haloes (Kormendy & Ho 2013), virial shock heating alone could help suppress accretion of cold gas from the cosmic web (Dekel & Birnboim 2006). It is likely that AGN feedback becomes important at some point to prevent reaccretion of ejected gas (Croton et al. 2006), as the energy input from star formation falls far short of the level required to completely unbound the gas from the galaxy and maintain a quenched state. Indeed, we expect that compact galaxies quench abruptly (Valentino et al. 2020) and that ‘maintenance mode’ (or ‘jet mode’) AGN feedback tends to be observed after galaxies have already begun to quench (Hardcastle, Evans & Croston 2007).

Observations of CO (1–0) in high-redshift lensed AGN starbursts have revealed disturbed morphologies, suggestive of feedback (AGN or stellar) or ongoing major mergers (Thomson et al. 2012; Spingola et al. 2020), while others do not appear to show such features (Riechers et al. 2011b; Sharon et al. 2016). Such observations for this sample will be useful to understand to what extent the CO lines investigated here probe the mass and kinematics of the gas reservoir, particularly in comparison to the population of DSFGs that do not appear to host rapidly accreting AGNs.

5.5 Selection effects and confusion

A source of caution in the interpretation of our findings is the combination of systematic biases that stem from the comparison of samples with different selection effects. Ideally, we would like to compare coeval samples of DSFGs and quasar hosts. However, the detection of DSFGs is strongly influenced by selection effects, such that systems with lower dust temperatures are preferentially detected at lower redshifts due to the shape of the modified blackbody spectrum. The selection of strongly lensed DSFGs at FIR wavelengths is much less dependent on dust temperature and expected to be more strongly dependent on source size (Hezaveh et al. 2012; Serjeant 2012). Together, these represent a diverse selection of DSFGs. The quasar systems, on the other hand, were optically selected, based on the emission from their unobscured accretion discs, and their selection was in most cases based on their lens configuration (RX J0911+0551 and H1413+117 being exceptions). Therefore, we do not implicitly select quasar hosts that have high temperatures or are more compact. However, analyses with larger samples, including less strongly magnified systems, will provide a better statistical comparison between coeval galaxy populations selected with different methods.

Notably, we identify a nearby source of one object (H1413+117) that likely causes photometric confusion at shorter wavelengths. This has the important implication that source blending may contribute to the extreme infrared luminosities found for some quasar starbursts. Observations with ALMA of (unlensed) optically selected quasars have revealed that 30 per cent of FIR identifications can be resolved into secondary counterparts that contribute at least 25 per cent to the measured infrared luminosity (Hatziminaoglou et al. 2018). This may be unsurprising as massive galaxies seem to be formed in overdense regions (Zeballos et al. 2018). The sources investigated here are gravitationally lensed, so the relative contribution from companions is expected to be lower on average than for field sources. Nevertheless, our findings suggest care should be taken to account for field sources when compiling an SED with photometric measurements obtained with low spatial resolution (e.g. Herschel/SPIRE).

6 CONCLUSIONS

The properties of cold dust and molecular gas in quasar hosts at cosmologically important redshifts can give insights into the mechanism responsible for the transformation of these galaxies from gas-rich, dusty starbursts into passive, quiescent systems. We have presented high angular resolution imaging with ALMA of the host galaxies of seven optically selected quasars at redshifts between 1.5 and 3. Using pixellated lens modelling, we reconstructed the sources and their intrinsic properties.

In comparison to unlensed DSFGs with similar redshifts and infrared luminosities (i.e. ALESS DSFGs and FIR-selected lensed DSFGs) and DSFGs at higher redshifts with similar dust temperatures, the quasar hosts in this work are generally more compact. The observed luminosities are broadly consistent with the expectations for more compact star formation, disputing the case for a significant contribution from black hole accretion to the global dust emission at sub-mm wavelengths.

We find that two of the quasar hosts are characterized by extended, clumpy dust distributions, but the remainder are compact starbursts. These differences may represent quasar hosts at different stages of morphological change into compact galaxies. The more compact systems have sizes that are already similar to quiescent galaxies at z ∼ 2, with extreme SFR densities that imply a rapid consumption of gas. This is consistent with a picture in which the inflow of gas resulting from mergers and/or dynamical instabilities triggers both the intense starburst and a period of efficient black hole accretion (Kocevski et al. 2017). These intense starbursts could be responsible for both the formation of the high stellar densities and rapid quenching of compact galaxies.

If quasar hosts are in the process of forming stellar bulges, we should find that these systems are globally characterized by unstable gas discs and rapidly depleting gas reservoirs (Zolotov et al. 2015; Tacchella et al. 2016). Matched, high angular resolution observations of low J-level CO emission would probe the size of the cold gas reservoir and determine whether the SFR efficiency and gas depletion time-scales of quasar hosts are significantly different relative to those of DSFGs. Furthermore, cold gas diagnostics will be useful to determine whether outflows are prevalent and whether AGN or stellar feedback is ongoing, as expected from simulations and semi-analytical models (e.g. Hopkins et al. 2008).

Our conclusions here are based on the data for just seven objects from the sample of lensed quasars. In the current observing cycle of ALMA, we will obtain data for at least 27 lensed quasars in Bands 6 and 7 which will further probe the size and structure of the heated dust emission from this class of objects, and probe the structure and kinematics of the molecular gas that is feeding both the star formation and AGN activity.

ACKNOWLEDGEMENTS

JPM acknowledges support from the Netherlands Organization for Scientific Research (NWO, project number 629.001.023) and the Chinese Academy of Sciences (CAS, project number 114A11KYSB20170054). SV has received funding from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation programme (grant agreement No 758853). CS is grateful for support from the National Research Council of Science and Technology, Korea (EU-16-001).

This research made use of astropy and matplotlib packages for python (Hunter 2007; The Astropy Collaboration 2018).

This paper makes use of the following ALMA data:

ADS/JAO.ALMA#2012.1.00175.S, ADS/JAO.ALMA#2015.1.01309.S, ADS/JAO.ALMA#2017.1.01368.S, ADS/JAO.ALMA#2017.1.01081.S, ADS/JAO.ALMA#2018.1.01519.S. ALMA is a partnership of ESO (representing its member states), NSF (USA), and NINS (Japan), together with NRC (Canada), MOST and ASIAA (Taiwan), and KASI (Republic of Korea), in cooperation with the Republic of Chile. The Joint ALMA Observatory is operated by ESO, AUI/NRAO, and NAOJ.

DATA AVAILABILITY

The data underlying this article will be shared on reasonable request to the corresponding author.

Footnotes

REFERENCES