Connection between steep radio spectral slopes and dust extinction in QSOs: evidence for outflow-driven shocks in dusty QSOs Open Access

ABSTRACT

Recent studies have found a striking positive correlation between the amount of dust obscuration and enhanced radio emission in quasi-stellar objects (QSOs). However, what causes this connection remains unclear. In this paper we analyse uGMRT Band-3 (400 MHz) and Band-4 (650 MHz) data of a sample of 38 1.0 < z < 1.5 QSOs with existing high-resolution |$0{_{.}^{\prime\prime}} 2$| e-MERLIN 1.4 GHz imaging. In combination with archival radio data, we have constructed sensitive 4–5 band radio spectral energy distributions (SEDs) across 0.144–3 GHz to further characterize the radio emission in dusty QSOs. We find that the dusty QSOs [those with |$E(B-V)$| > 0.1 mag] are more likely to exhibit steep spectral slopes (⁠|$\alpha$| < |$-0.5$|⁠; S |$\propto$| |$\nu ^{\alpha }$|⁠) than the non-dusty QSOs [|$E(B-V)$| < 0.1 mag], with fractions of 46|$\pm$|12 and 12|$\pm$|4 per cent, respectively. A higher fraction of the non-dusty QSOs have peaked radio SEDs (48|$\pm$|9 per cent) compared to the dusty QSOs (23|$\pm$|8 per cent). We discuss the origin of the radio emission, finding that the majority of the peaked, predominantly non-dusty, QSOs have consistent sizes and luminosities with compact jetted radio galaxies. However, the connection between steepness and dust obscuration implies an outflow-driven shock origin for the enhanced radio more commonly found in dusty QSOs. These results add to the emerging picture whereby dusty QSOs are in an earlier blow-out phase, with shocks that heat and destroy the surrounding dust, eventually revealing a typical non-dusty QSO.

1 INTRODUCTION

Quasi-stellar objects (QSOs) are the most powerful class of active galactic nuclei (AGNs), and are often classified by their extremely high bolometric luminosities (⁠|$10^{45-48}$| erg s|$^{-1}$|⁠). The majority of QSOs are optically very blue due to an unobscured view of the accretion disc which has emission that peaks in the optical/ultraviolet (UV). However, a subset have been found to display much redder colours (‘red QSOs’; Webster et al. 1995). There are several different explanations for the red colours in these QSOs, such as a moderate viewing angle through the dusty torus (Wilkes et al. 2002; Rose et al. 2013), as suggested by the unification model of AGN (Urry & Padovani 1995), an intrinsically red continuum (Whiting, Webster & Francis 2001; Young, Elvis & Risaliti 2008), an unusual covering factor of hot dust (Rose 2014), a strong synchrotron component (Whiting et al. 2001), or dust within the host-galaxy extinguishing the emission at shorter wavelengths (Glikman et al. 2012; Kim & Im 2018; Fawcett et al. 2022; Kim et al. 2024b). Interestingly, simulations of galaxy formation predict that these red QSOs could be in a short-lived ‘blow-out’ phase, whereby powerful outflows interact and shock the surrounding dust and gas, eventually clearing it and revealing a typical unobscured (blue) QSO (Hopkins et al. 2006, 2008). If red QSOs do indeed represent a transitional phase in galaxy evolution then they would provide excellent laboratories to study AGN ‘feedback’ in action (see review by Harrison & Ramos Almeida 2024).

Radio observations provide a powerful tool for distinguishing between different red QSO scenarios, since radio emission is not affected by dust. Therefore, if red and blue QSOs are intrinsically the same objects (as suggested by the orientation scenario) then we would not expect to observe any differences in the radio emission (apart from potentially enhanced radio emission in blue QSOs due to Doppler boosting of a more face-on radio jet; e.g. Lähteenmäki & Valtaoja 1999). Indeed, evidence supporting the red QSO blow-out phase scenario over the orientation scenario has largely come from differences found in the radio properties compared to blue QSOs. For example, studies have found that red QSOs are more likely to be detected in the radio compared to redshift and luminosity-matched blue QSOs, which has been observed with radio data across a large range of frequencies, spatial resolutions, and sensitivities (Richards et al. 2003; Georgakakis et al. 2009; Banerji et al. 2012; Glikman et al. 2012; Klindt et al. 2019; Fawcett et al. 2020; Rosario et al. 2020; Fawcett et al. 2021; Glikman et al. 2022). This enhanced radio emission in red QSOs has been found to be predominantly compact, on host-galaxy scales (i.e. |$\lesssim$| 10 kpc) and tends to be radio-quiet/intermediate (i.e. they do not host large-scale powerful radio jets; Klindt et al. 2019; Fawcett et al. 2020; Rosario et al. 2021). Previous work has also found that star formation is unlikely the dominant mechanism, especially for more luminous QSOs (Fawcett et al. 2020; Rosario et al. 2020; Calistro Rivera et al. 2021; Andonie et al. 2022; Yue et al. 2024). These results point towards either low-powered radio jets (e.g. Girdhar et al. 2022) or outflow-driven shocks (e.g. Haidar et al. 2024; Stepney et al. 2024) as the origin of the radio emission in red QSOs [see Panessa et al. (2019) for a review on different radio emission mechanisms].

Building on this previous work, Fawcett et al. (2023) explored the radio detection rate of QSOs from the Dark Energy Spectroscopic Instrument (DESI; DESI Collaboration 2016a, b) as a function of dust extinction and found a striking positive correlation (confirmed by Calistro Rivera et al. 2024; Petley et al. 2024). This suggests a causal link between opacity and the production of radio emission in QSOs, which is likely due to shocks from outflows, which may be driven by radiation pressure, accretion disc winds, and/or jets on the surrounding interstellar medium (ISM); see review by Harrison & Ramos Almeida (2024). These shocks will heat and eventually destroy the dust, reducing the obscuration and therefore make the QSO appear bluer. This scenario is further supported by Calistro Rivera et al. (2021), who found a hot dust excess and stronger [O iii]|$\lambda$|5007 outflows in the red QSOs, suggesting dust is getting heated by outflow-driven shocks. Indeed, this may be a more extreme example of the spatially resolved connection observed between outflows, dust, and shocks in nearby AGN (e.g. Haidar et al. 2024). If synchrotron radiation from a shocked dusty environment was found to be the driving mechanism behind the radio emission in red QSOs, then this would present strong evidence for the dusty red QSO blow-out phase.

One of the most effective ways to determine the radio emission mechanism is high spatial resolution radio imaging (e.g. very long baseline interferometry; VLBI). For example, a highly collimated structure and/or hotspots would likely be associated with a radio jet, whereas a less collimated, more diffuse structure is likely to be an outflow or star formation (Panessa et al. 2019; Kharb et al. 2021; Chen et al. 2024; Njeri et al. 2024). However, this technique has its limitations. For example, if the jet is inclined towards the accretion disc then this can drastically alter the morphology of the resulting radio emission, making the origin more unclear (Mukherjee et al. 2018; Meenakshi et al. 2024). To compare the |$\sim$| kpc scale radio morphology of red and blue QSOs, Rosario et al. (2021) obtained |$0{_{.}^{\prime\prime}} 2$| enhanced Multi Element Remotely Linked Interferometer Network (e-MERLIN) 1.4 GHz images of 20 red and 20 blue (selected based on |$g-i$| colours) luminous (⁠|$L_{\rm bol}\approx 10^{46-47}$| erg s|$^{-1}$|⁠) QSOs at intermediate redshifts (1.0 < z < 1.55) with unresolved, detected radio emission in the Faint Images of the Radio Sky at Twenty-centimeters (FIRST; Becker, White & Helfand 1995; 5 arcsec resolution; |$L_{\rm 1.4\, GHz}\approx 10^{25-26}$| W Hz|$^{-1}$|⁠). At these redshifts the host-galaxy scale radio emission could be probed; they found a statistically significant difference in the incidence of |$\sim$| 2–10 kpc-scale extended radio emission in the red QSOs compared to the blue QSOs which indicated that the enhanced radio emission in red QSOs was likely due to outflows or low-powered jets that are confined within the host galaxy (Rosario et al. 2021). However, despite the impressive angular resolution of e-MERLIN (25|$\times$| better than FIRST) the majority of the sample (⁠|$\sim$| 85 per cent) remained unresolved.

To understand what mechanism drives these |$\sim$|kpc-scale unresolved radio structures (and to gain additional insight into the nature of the resolved sources), analysing radio spectral energy distributions (SEDs) can be an effective method. For example, older radio emission produced by shocks from jets or winds (optically thin synchrotron emission) would result in a steep radio spectral slope (Faucher-Giguère & Quataert 2012; Nims, Quataert & Faucher-Giguère 2015; Laor, Baldi & Behar 2019) and a young, self-absorbed radio nucleus (optically thick synchrotron emission) would result in an inverted or flat spectrum (Blandford & Königl 1979). In particular, a well-studied class of objects known as compact steep spectrum (CSS) and gigahertz-peaked spectrum (GPS) sources can be identified by a spectral turnover in their SED around |$\sim$| 100 MHz and |$\sim$| 1 GHz, respectively (see reviews by O’Dea & Baum 1997; O’Dea & Saikia 2021). These systems likely host compact jets that are thought to be either ‘frustrated’ (i.e. confined by the galaxy ISM) on scales of |$\sim$|0.1–2 kpc (van Breugel, Miley & Heckman 1984; O’Dea, Baum & Stanghellini 1991; Orienti 2016), in a young evolutionary phase (Phillips & Mutel 1982; Carvalho 1985; Bicknell et al. 2018), or both (Patil et al. 2020). The frequency of the spectral turnover has been found to strongly anticorrelate with the projected linear size of the radio emission (⁠|$\nu$| |$\propto$| |$l^{-0.65}$|⁠; e.g. Fanti et al. 1995; O’Dea & Baum 1997; Orienti & Dallacasa 2013) which is expected for synchrotron self-absorption (Snellen et al. 2000).¹ Therefore, a peak in the radio SED of a QSO would likely suggest the presence of a compact radio jet.

Many studies have explored the radio spectral slopes of different classes of QSOs (e.g. Callingham et al. 2017; Laor et al. 2019; Zajaček et al. 2019; Shao et al. 2022; Kukreti et al. 2023; Hayashi, Doi & Nagai 2024), including between red and blue QSOs (Georgakakis et al. 2012; Rosario et al. 2020; Glikman et al. 2022; Sargent et al., in preparation). For example, utilizing data from FIRST and the VLA Sky Survey (VLASS; Gordon et al. 2020), Glikman et al. (2022) explored the 1.4–3 GHz spectral slopes and found that infrared (IR)-selected red QSOs were on average steeper compared to blue QSOs. On the other hand, utilizing data from the LOw-Frequency ARray (LOFAR; van Haarlem et al. 2013), Two-metre Sky Survey (LoTSS; Shimwell et al. 2017), and FIRST, Rosario et al. (2020) found no differences between the lower frequency 0.144–1.4 GHz radio spectral slopes of optically selected red and blue QSOs. This could suggest that the radio spectral slope of red QSOs only deviates at higher frequencies or that differences only arise in the more heavily reddened QSOs studied in Glikman et al. (2022). However, both of these studies are limited to two data points in frequency. Without additional data points information can be lost; e.g. whether or not there is a spectral turnover and, if so, the exact location of this turnover (Patil et al. 2022). Radio observations at additional frequencies can be acquired, but usually only for small samples. Therefore, a combination of both approaches is needed to fully understand the nature of the radio emission in different QSO populations.

In this paper we construct sensitive radio SEDs from 144 MHz–3 GHz of the QSOs that have been observed at |$0{_{.}^{\prime\prime}} 2$| with e-MERLIN (Rosario et al. 2021), by combining dedicated observations at 400 and 650 MHz from the upgraded Giant Metrewave Radio Telescope (uGMRT; Ananthakrishnan 1995; Gupta et al. 2017) and archival radio data. Using these data we can carefully explore the differences, if any, between the radio SEDs of dusty red and typical blue QSOs. If found, this would indicate differences in either the driving mechanism (i.e. outflow-driven shocks, frustrated jets, etc.) or evolution (i.e. young/old) of the radio emission. The presence of a spectral turnover can also determine whether the QSOs are similar to GPS/CSS-like sources, by comparing their peak frequency and measured e-MERLIN sizes to the known anticorrelation for GPS/CSS sources. Finally, finding steep radio spectral slopes in the dusty QSOs would be consistent with the shocked dust scenario as the origin of the enhanced radio emission found in red QSOs. In Section 2 we describe the sample selection, data utilized in this paper, and SED fitting procedure. In Section 3 we present our results and in Section 4 we discuss the origin of the radio emission. In this paper we define spectral index |$\alpha$| as |$S_{\nu }$| |$\propto$| |$\nu ^{\alpha }$|⁠, where |$S_{\nu }$| is the flux density at frequency |$\nu$|⁠. Throughout our work we adopt a standard flat |$\lambda$|-cosmology with |$H_0$| |$=$| 70 km s|$^{-1}$|Mpc|$^{-1}$|⁠, |$\Omega$|_M |$=$| 0.3, and |$\Omega _{\lambda }$| |$=$| 0.7 (e.g. Planck Collaboration VI 2020).

2 METHODS

In this paper we explore the radio SEDs of Sloan Digital Sky Survey (SDSS)-selected QSOs with complementary e-MERLIN 1.4 GHz imaging (Rosario et al. 2021; hereafter, R21). To construct the radio SEDs, we combined observations from the uGMRT in addition to archival radio data sets. In the following sections we describe the sample selection (Section 2.1), the uGMRT observations and data reduction process (Section 2.2), the various archival radio data sets utilized (Sections 2.3 and 2.4), how radio-loudness and radio-luminosity are calculated (Section 2.5), the radio SED fitting method (Section 2.6), and the dust extinction fitting method (Section 2.7).

2.1 Sample selection

The sample of red QSOs (hereafter, ‘rQSOs’) and control blue QSOs (hereafter, ‘cQSOs’) observed with uGMRT were chosen to have existing |$0{_{.}^{\prime\prime}} 2$| e-MERLIN observations at 1.4 GHz (R21). The sample selection is presented in detail in R21, which we describe briefly below. The sample properties can be found in Table 1.

Starting with the SDSS DR7 (Schneider et al. 2010) QSO catalogue, a sample of rQSOs and cQSOs were chosen as the top 10 percentile and middle 50 percentile of the |$g-i$| distribution, respectively, in consecutive redshift bins of 1000 sources (see top panel of Fig. 1), following the same method as Klindt et al. (2019); hereafter, K19. The QSOs were also selected to be detected in FIRST, with no visibly resolved radio emission beyond the 5 arcsec resolution. In order to robustly determine the radio morphology via visual inspection, additional cuts were applied: a 1.4 GHz flux density of > 3 mJy, a peak flux S/N > 15, and a flux ratio between FIRST and the 20 cm NRAO VLA Sky Survey (NVSS; 1.4 GHz at 45 arcsec resolution; Condon et al. 1998) of less than two. These criteria ensured a robust detection in FIRST and the removal of any sources with radio lobes detected in NVSS but resolved out by the higher resolution FIRST imaging K19.²

A narrow redshift range (1.0 < z < 1.55) was chosen to minimize Malmquist bias. A radio luminosity range of 25.5 W Hz⁻¹ < log|$L_{\rm 1.4\, GHz}$| < 26.5 W Hz|$^{-1}$| was chosen to probe the fainter end of FIRST-detected QSOs at the chosen redshift. Finally, R21 selected 20 rQSOs and 20 cQSOs from this final QSO sample, matched in redshift and |$L_{\rm 6\mu m}$|⁠, using a matching tolerance of 0.05 and 0.2 dex, respectively (following the method from K19; see Section 2.3 therein). One cQSO (1511+3428) was later found to have FIRST radio emission slightly offset from the optical centre, which therefore should have been classified as ‘extended’. Furthermore, a re-analysis of the Galactic extinction of one rQSO (0842+5804) placed its intrinsic colour outside of the red colour selections. In this paper we remove both of these QSOs from our analyses, resulting in a final sample of 19 rQSOs and 19 cQSOs. The redshift, |$L_{\rm 1.4\, GHz}$|⁠, and |$L_{\rm 6\mu m}$| distributions of the final samples are displayed in Fig. 1.

2.2 uGMRT observations and data reduction

In this paper, we present the upgraded GMRT (uGMRT; Swarup et al. 1991; Gupta et al. 2017) Band-3 (250–500 MHz; central frequency 400 MHz) and Band-4 (550–850 MHz; central frequency 650 MHz) observations taken in February 2020 (Proposal ID: 37_064, PI: V. Fawcett).

Due to two and three antennas not working for the Band-4 and Band-3 observations, respectively, 28 and 27 antennas were used for the observations with the full available bandwidth (200 MHz for Band-3 and 400 MHz for Band-4) and 4096 channels. The flux density calibrators 3C147 and 3C286 were observed for 10 min at the beginning and end of each observation. The phase calibrators (listed in the online Supplementary material) were observed for 5 min every 2–8 targets, depending on their on-sky position. The QSO targets were observed for 5 min in six observing blocks, using the same grouping as that used for the e-MERLIN observations (see section 2.4 of R21). For one rQSO (1315+2017) there is no Band-3 observation.

The data were reduced in CASA-6³ using CAPTURE⁴ (Kale & Ishwara-Chandra 2021), an automated pipeline to produce images from interferometric data obtained with the uGMRT. The LTA data files produced by the uGMRT were first converted into FITS format by the LISTSCAN and GVFITS functions. A measurement set was then produced using the IMPORTUVFITS task. The standard steps of flagging, calibration, and imaging were then followed in CAPTURE. Self-calibration with both phase-only and amplitude and phase solutions was carried out in eight iterations. The flux scale of Perley & Butler (2017) was used for the absolute flux calibration. For all sources, imaging was carried out with robust = 0, a pixel scale of 1.5 and |$1\, \rm arcsec$| for Band-3 and 4, respectively, and an image size of 8000 pixels (88.9 arcmin and 133.3 arcmin for Band-3 and 4, respectively). Due to the low flux density of some of the targets, the mean flux cut-off for flagging bad antennas in Band-3 was reduced from the default value of 0.3 to 0.08 Jy in the CAPTUREugfunctions script. This was required since for some sources > 70 per cent of the antennas were getting flagged. We note that although including noisier data in the imaging might result in lower signal-to-noise (SNR) images, all our sources have a final SNR > 6 in both Band-3 and Band-4. The median flagging fraction for all the QSOs in Band-3 and Band-4 is 24 per cent and 9 per cent, respectively.

Finally, flux density measurements were obtained using the IMFIT task in CASA, which fits 2D elliptical Gaussians to the image. For sources with extended radio lobes only the flux in the radio core was extracted. A typical root mean square (RMS) of |$\sim$| 0.17 and |$\sim$| 0.06 mJy was achieved for Band-3 and 4, respectively. The online Supplementary material contains a table with the observation details and extracted flux values for the Band-3 and 4 data for our sample.

2.3 Archival radio data used in SED fitting

In order to robustly measure the radio SEDs of the QSOs, we combined our 2-band uGMRT data with various archival radio data. At the redshifts of our sample, all the radio data used in the SED fitting are for galaxy-wide emission (⁠|$\sim$| 20–50 kpc) and sources with extended low-frequency emission are treated carefully (see Section 3.2).⁵ The archival radio data used in this paper is displayed in Table 2 and summarized below.

For the final radio SEDs, we utilized the integrated flux density, rather than the peak flux density, in order to better capture the flux density for sources that might be extended at lower radio frequencies (a comparison between the FIRST peak and integrated flux densities is shown in Fig. B1). We carefully treat sources that are only extended in LoTSS separately (see Section 3.2).

2.3.1 LoTSS DR2

The LOFAR Two Metre Sky Survey DR2 (LoTSS; Shimwell et al. 2022) is a 120–168 MHz LOFAR radio survey that aims to observe the whole northern sky at a resolution of |$6\, \rm arcsec$|⁠. The second data release covers 5740 deg|$^{2}$| down to a sensitivity limit of |$\sim$| 83 |$\mu$|Jy beam|$^{-1}$|⁠. In this paper we utilize the catalogue with associated optical and/or near-infrared counterparts from Hardcastle et al. (2023), which contains 4116 934 sources. 14 rQSOs and 13 cQSOs lie in the LoTSS DR2 sky coverage, and 12/14 rQSOs and 13/13 cQSOs are detected (⁠|$0{_{.}^{\prime\prime}} 2$| matching radius with optical positions, consistent with Rosario et al. 2020).

2.3.2 FIRST

The FIRST (Becker et al. 1995) is a 1.4 GHz VLA radio survey, covering 10 575 deg|$^{2}$| of the sky in the SDSS region at a resolution of |$5\, \rm arcsec$|⁠. The final catalogue (Helfand, White & Becker 2015) contains 946 432 sources down to a sensitivity limit of |$\sim$| 0.65 mJy beam|$^{-1}$|⁠. As our sample selection requires a robust FIRST detection, all 38 QSOs have FIRST flux densities (⁠|$10\, \rm arcsec$| matching radius, consistent with K19).

2.3.3 VLASS

The VLA Sky Survey (VLASS; Gordon et al. 2020) is a 2–4 GHz VLA radio survey, covering 33 885 deg|$^{2}$| of the sky at a resolution of |$2{_{.}^{\prime\prime}} 5$|⁠. In this paper, we utilize the Quick Look Epoch 2 catalogue, which consists of 2995 025 sources down to a sensitivity of 1 mJy beam|$^{-1}$|⁠. We included an additional error of 10 per cent of the integrated flux density to account for the known flux underestimation (Lacy et al. 2020; Gordon et al. 2021). All 38 QSOs are detected (⁠|$1\, \rm arcsec$| matching radius). We note that only five QSOs have a |$F_{\rm 3\, GHz}$| < 3 mJy, which are known to have more unreliable flux densities (Gordon et al. 2021); for these sources, we add instead an additional systematic error of 20 per cent of the integrated flux density and flag these sources in Table 3. We also note that in our SED fitting we utilize the integrated flux densities from the Epoch 2 catalogue, which are known to be more reliable than those from Epoch 1 and more reliable than the peak flux densities.⁶

In order to test the impact of radio variability on our study, we compared the integrated flux densities from both the Epoch 1 and 2 Quick Look tables; we found no significant variation in the VLASS integrated fluxes for any of our sample (Fig. A1). For more discussion on radio variability, see Appendix  A.

2.3.4 TGSS ADR1

The TIFR Giant Meterwave Radio Telescope Sky Survey Alternative Data Release (TGSS ADR1; Intema et al. 2017) is a 140–156 MHz radio survey covering 36 900 deg|$^{2}$| of the sky at a resolution of |$25\, \rm arcsec$|⁠. The |$7\sigma$| catalogue contains 623 604 sources down to a sensitivity limit of |$\sim$| 5 mJy beam|$^{-1}$|⁠. Using a |$10\, \rm arcsec$| matching radius, we found that four rQSOs and two cQSOs are detected. Despite the low angular resolution of this survey, there are two sources that lie outside of the LoTSS coverage but are detected by TGSS; including these additional data points greatly improves the fitting constraints (we indicate these sources in Table 3). There are a further nine sources that lie outside of the LoTSS coverage, but have upper limits in TGSS. The majority of these upper limits do not affect the fitting due to the low sensitivity limit of TGSS; however, there is one source for which the best-fitting model without the TGSS upper limit is inconsistent with the upper limit (1003+2727).

2.3.5 LoLSS DR1

The LOFAR⁷ LBA Sky Survey DR1 (LoLSS; de Gasperin et al. 2023) is a 41–66 MHz LOFAR radio survey that aims to observe the whole northern sky above declination 24|$^{\circ }$| at a resolution of |$15\, \rm arcsec$|⁠. The first data release covers 650 deg|$^{2}$| in the Hobby-Eberly Telescope Dark Energy Experiment (HETDEX) Spring Field and contains 42 463 sources down to a sensitivity limit of |$\sim$| 1.55 mJy beam|$^{-1}$|⁠. Only one rQSO (1153+5651) lies in the LoLSS DR1 sky coverage and is also detected (⁠|$10\, \rm arcsec$| matching radius); despite the low angular resolution, the LoLSS data point is in agreement with the higher frequency data for this source and so we used this data in the fitting for an additional constraint.

2.4 Supplementary archival radio data

For our final radio SEDs (see Section 3) we also plot additional archival radio data including the e-MERLIN |$0{_{.}^{\prime\prime}} 2$| radio fluxes (Section 2.4.1) and radio data from surveys with a much lower angular resolution compared to those listed above (Section 2.4.2). Although part of the selection criteria ensured that the flux offset between FIRST (at |$5\, \rm arcsec$|⁠) and NVSS (at |$45\, \rm arcsec$|⁠) was less than two, biases may still be introduced in the SED fitting if lower resolution data are introduced due to potential large-scale radio lobes missed in the higher resolution surveys or additional contributions from additional background sources within the beam. Furthermore, the e-MERLIN observations have a considerably higher angular resolution (> 12|$\times$|⁠) compared to the other radio data utilized in this paper, and so might miss diffuse extended structures. Therefore, we do not include the following radio data in the SED fitting. Instead, these data provide another constraint on radio variability, which we discuss in Appendix  A, in addition to providing radio sizes (from e-MERLIN) and helping to visually asses the reliability of the SED fits.

2.4.1 e-MERLIN L-band

The sample used in this paper was previously observed using the e-MERLIN interferometer with the L-band receivers (1.23–1.74 GHz) at a |$0{_{.}^{\prime\prime}} 2$| resolution (PI: D. Rosario; Project ID: CY7220). The e-MERLIN flux densities were obtained from R21, who fit either one or multiple Gaussian components to the central |$2\, \rm arcsec$| of the e-MERLIN image, depending on a visual examination of the fitting residuals. For the first round of fitting, a single Gaussian model was used which was initialized to match the shape of the restoring beam. After inspecting the resulting residuals, seven sources (six used in this paper, see Section 2.1) were found to have emission extended beyond the beam which were then subject to a second round of fitting with > 2 Gaussian components (e-MERLIN extended sources; see Tables 1 and 3). For the remaining 32, a further analysis was performed, comparing the semimajor axis of the single-component Gaussian fits with that of the restoring beam to search for any additional extended source structure at the resolution limit of the e-MERLIN images. R21 found no statistical difference between the semimajor axis of the Gaussian component and the restoring beam, demonstrating that these sources are unresolved in the e-MERLIN images (e-MERLIN unresolved sources).

The flux ratio between FIRST and e-MERLIN can provide information on whether diffuse extended emission (at the same frequency) from radio lobes that has been resolved out; R21 found no significant differences between the flux densities for either the cQSOs or rQSOs (see section 3.2 in R21). They also note that for the flux estimates for diffuse e-MERLIN sources are likely to have larger errors and may be overestimated, which might explain the |$\approx$| 2.5 |$\times$| larger e-MERLIN flux compared to the FIRST flux for QSO 0946+2548. We plot the e-MERLIN flux densities as open stars on the SEDs (see Supplementary material and Fig. 2).

The e-MERLIN radio sizes were calculated by taking the maximum of either the largest separation between the two Gaussian subcomponents or the major-axis width of the largest single component. For sources that only have a single core component, its major-axis is used as a limit on the size (i.e. e-MERLIN unresolved sources). These were then converted to physical sizes using the angular diameter distances R21.

2.4.2 Lower angular resolution radio surveys

The Westerbork Northern Sky Survey (WENNS; Rengelink et al. 1997) is a 322.5–327.5 MHz radio survey, covering the whole of the sky above a declination 30|$^{\circ }$| at a resolution of |$54\, \rm arcsec$|⁠. The |$5\sigma$| combined catalogue contains 229 420 sources from the mini and main surveys down to a sensitivity limit of 18 mJy. Using a |$25\, \rm arcsec$| matching radius we find that five rQSOs and three cQSOs are detected.

The Rapid ASKAP Continuum Survey (RACS; McConnell et al. 2020) is the first large-area survey with the Australian Square Kilometer Array Pathfinder (ASKAP; Johnston et al. 2007; Hotan et al. 2021) which aims to observe the entire southern sky at a frequency of 700–1800 MHz with |$\sim$| |$25\, \rm arcsec$| resolution. The first release |$5\sigma$| catalogue (Hale et al. 2021) contains 2123 638 sources at a central frequency of 887.5 MHz (288 MHz bandwidth) down to a sensitivity limit of 0.25–0.3 mJy beam|$^{-1}$|⁠. Using a |$15\, \rm arcsec$| matching radius we find that seven rQSOs and seven cQSOs are detected.

The 20 cm NRAO VLA Sky Survey (NVSS; Condon et al. 1998) is a 1.4 GHz at radio survey covering the whole sky above a declination of |$-40^{\circ }$| at |$45\, \rm arcsec$| resolution. The |$5\sigma$| catalogue contains 1773 484 sources down to a sensitivity limit of |$\sim$| 2.5 mJy beam|$^{-1}$|⁠. By selection (see Section 2.1), all the QSOs in our sample are detected in NVSS.

2.5 Radio-loudness and luminosity

In this paper we utilize |$L_{\rm 1.4\, GHz}$|⁠, which is calculated from the FIRST integrated fluxes, using the methodology described in Alexander et al. (2003), assuming a uniform radio spectral index of |$\alpha$| |$=$| |$-0.5$| for the K-correction (Fig. 1). We also explored how |$L_{\rm 1.4\, GHz}$| changes once adopting our more robust individual values of |$\alpha$|⁠, obtained from our SED fitting (Section 3.2; median |$\alpha _{\rm PL}$| |$=$| |$-0.4$|⁠), and found a median absolute difference of 0.19 dex for the power-law sources.⁸ To explore the ‘radio-loudness’ of our samples, we adopted the same parameter as that first used in K19, defined as the dimension-less quantity:

By this definition, the radio-loud/radio-quiet threshold is |$\mathcal {R}=-4.2$|⁠, which is broadly consistent with the canonical definition often defined as the 5 GHz-to-2500 Å ratio (e.g. Kellermann et al. 1989). We utilize the 6 |$\mu$|m luminosity rather than the optical luminosity since this is less susceptible to obscuration from dust (see K19 for full details).⁹ K19 showed that the differences in the radio properties between rQSOs and cQSOs arose around radio-loudness values of |$-5$| < |$\mathcal {R}$| < |$-3.7$| (i.e. radio-quiet/radio-intermediate values; also demonstrated in Fawcett et al. 2020, 2021; Rosario et al. 2020). Therefore, in this study we explored our results when splitting at |$\mathcal {R}$| |$=$| |$-3.7$|⁠, referring to QSOs with a |$\mathcal {R}$| < |$-3.7$| as ‘radio-quiet/radio-intermediate’ and |$\mathcal {R}$| > |$-3.7$| as ‘extreme radio-loud’ (see the lower panel of Fig. 1).

2.6 Radio spectral fitting and characterization

In this paper we aim to characterize the radio SEDs of our sample of QSOs, in order to determine how many are best fit by a typical synchrotron power law (either continuous or broken) or display a curved, peaked SED. Furthermore, we can compare the dusty red QSOs and typical blue QSOs to understand whether they have different radio spectral properties. Three example SED fits are displayed in Fig. 2, with one source favouring a flat power-law model, one source favouring a steep broken power-law model, and the other favouring a peaked model.

Therefore, in order to model the radio spectral properties of our sample, we adopted three different spectral models to fit the uGMRT plus archival radio data (e.g. Patil et al. 2022; Kerrison et al. 2024). To fit the three models to the data we used the emcee¹⁰ package in Python (Foreman-Mackey et al. 2013).

The first model we used is a standard non-thermal power-law model:

where a is the amplitude of the synchrotron spectrum, |$\nu$| is the frequency, and |$\alpha _{\rm PL}$| is the spectral index. The bounds for |$\alpha _{\rm PL}$| were set at |$\pm 3$|⁠.

The second model we used is a broken power-law (BPL) model to characterize a standard optically thin synchrotron source with additional energy losses at higher frequencies (due to synchrotron and inverse-Compton losses; e.g. Panessa et al. 2022):

where |$\alpha _{\rm low}$| and |$\alpha _{\rm high}$| are the spectral indices of the power-law below and above the break frequency |$\nu _{\rm break}$|⁠, respectively. Since energy losses steepen the spectrum, we required that |$\alpha _{\rm high}$| < |$\alpha _{\rm low}$| and set the bounds for both power-low indices to be |$-3$| < |$\alpha _{\rm PL}$| < 0. We required |$\nu _{\rm break}$| to be within the frequency bounds of the radio data utilized in the SED fitting. There is some degeneracy between the two power-law indices and |$\nu _{\rm break}$|⁠, which is reflected in the resulting Monte–Carlo errors.

The third model we used is the commonly adopted generic curved model (e.g. Snellen et al. 1998; Callingham et al. 2017; Shao et al. 2020; Wołowska et al. 2021; Kerrison et al. 2024), which is able to characterize a peaked-spectrum source:

where |$S_{\rm peak}$| is the flux density at the peak frequency, |$\nu _{\rm peak}$| and |$\alpha _{\rm thick}$| and |$\alpha _{\rm thin}$| are spectral indices in the optically thick and thin regime, respectively (e.g. Snellen et al. 1998). The bounds for |$\alpha _{\rm thick}$| and |$\alpha _{\rm thin}$| were set at < 4 and > |$-4$|⁠, respectively, following Kerrison et al. (2024). Therefore, |$\alpha _{\rm thick}$| corresponds to the lower frequency (⁠|$\nu$| < |$\nu _{\rm peak}$|⁠), inverted part of the peaked spectrum and |$\alpha _{\rm thin}$| corresponds to the higher frequency (⁠|$\nu$| > |$\nu _{\rm peak}$|⁠), steep part of the spectrum. When |$\alpha _{\rm thick}$| |$=$| 2.5 this model reduces to the homogeneous, synchrotron self-absorbed case. It should be noted that this function sometimes fails to accurately fit the SED far away from the peak (Snellen et al. 1998). However, due to our limited frequency range, this is unlikely to significantly affect the fits.

The best-fitting model was selected by utilizing the Bayesian information criterion (BIC; Schwarz 1978). Nominally, the fit with the minimum BIC value is taken to be the best fitting model. However, to favour the curve model over the other models, we required strong evidence that neither the PL nor BPL models were a sufficient fit to the data by imposing a |$\Delta$|(BIC) > 10 between the curve model and other models. Furthermore, due to the degeneracy between the BPL model and both the curve and PL models, we required a |$\Delta$|(BIC) > 5 between the BPL and other models to favour the BPL model. The final radio SEDs fits were further categorized as one of the following (for the BPL model, the |$\alpha$| adopted is |$\alpha _{\rm high}$| if the majority of data points have a |$\nu$| > |$\nu _{\rm break}$| and |$\alpha _{\rm low}$| vice versa):

• Steep: Best fit with the PL or BPL models [equations (2) and (3)] with |$\alpha$| < |$-0.5$|⁠,

• Flat: Best fit with the PL or BPL models with |$-0.5\lt $| |$\alpha$| < 0.5 (consistent with e.g. Urry & Padovani 1995; Kerrison et al. 2024),

• Inverted: Best fit with the PL model with |$\alpha _{\rm PL}$| |$\gt 0.5$|⁠,

• Peaked: Best fit with the curve model (equation 4) and a spectral turnover from the fitted model within the frequency range of the data utilized (0.144–3 GHz),

• Curved: Best fit with the curve model and a spectral turnover from the model that is predicted to be outside of frequency range of the data utilized.

After visually inspecting the resulting SED fits and BIC values, we found that four QSOs (0951+5253, 1046+3427, 1203+4510, and 1630+3847) have an ‘upturned’ radio SED, meaning that none of the models were a good fit to the data. For all four QSOs, these upturned SEDs were found to be driven by large-scale extended radio emission in LoTSS that is either not detected or not taken into account in the flux extraction for the higher frequency data. For example, QSO 1046+3427 displays large-scale radio emission in all frequency bands (see Fig. B2) that was not associated with the core emission in the FIRST catalogue but was associated as one source in LoTSS (see Appendix  B for more discussion on missed large-scale radio emission in two QSOs). Some studies attempt to model these sources by using a linear combination of equations (2) and (4) (e.g. Kerrison et al. 2024). However, in this study we are restricted to only five data points which would be less than the number of free parameters when combining these two equations. Therefore, for the four upturned QSOs we removed the LoTSS data point and refitted their radio SEDs in order to explore the compact emission (⁠|$\lesssim$| few kpc), consistent with the rest of the sources (an example is displayed in the bottom panel of Fig. 2). We found that all upturned sources have a good fit without the LoTSS data point, demonstrating that their SEDs are comprised of two components; a steep low frequency component and a flat/peaked high frequency component. We comment on these sources further in Section 4.4.

An additional QSO (1251+4317) has a good SED fit, apart from the VLASS data point which is over a magnitude brighter than that predicted by the best-fitting model to the other data. Interestingly, this source also has the largest offset between the measured e-MERLIN flux (⁠|$0{_{.}^{\prime\prime}} 2$|⁠) and FIRST flux (⁠|$5\, \rm arcsec$|⁠), with the e-MERLIN flux 2.4|$\times$| brighter than FIRST, despite the much higher angular resolution of e-MERLIN. This could suggest that this sources is undergoing rapid core variability. However, in Appendix  A we explore the VLASS variability on a 2-yr time-scale and do not find any significant differences for any of the sources. Other possibilities include the superposition of multiple radio components or a young restarted core (e.g. Nyland et al. 2020).

If the number of data points in the fitting is less than or equal to the number of parameters in the best-fitting model then we indicate the parameters on these fits as ‘unconstrained’. However, we are still able to distinguish which is the best-fitting model and therefore include these sources in the raw numbers for each best-fitting model (these sources are indicated in Table 3). We note that our fitting is restricted to the frequency range of the radio data used. Therefore, a radio SED best fit by a PL or BPL may be better fit by a curve model if lower/higher frequency radio data was utilized. We take this into account when exploring trends with |$\nu _{\rm peak}$| in Section 4.2. We present the results of this fitting approach in Section 3.2.

2.7 Dust extinction fitting

In order to explore the differences in the radio SEDs as a function of dust extinction, we quantified the amount of dust extinction in each of the QSOs. To do this we fit a blue QSO spectral template with varying amounts of dust extinction to the SDSS optical spectra, following the same method as Fawcett et al. (2022). We briefly outline the approach here.

We used the blue QSO VLT/X-shooter composite from Fawcett et al. (2022) as our unreddened template. We then masked the emission lines and smoothed the composite with a Gaussian filter. This blue template was then fitted to the SDSS spectra using a least-squares minimization code, which varied the dust extinction [|$E(B-V)$|], using a simple PL dust extinction curve (⁠|$A_V$| |$\propto$| |$\lambda ^{-1}$|⁠; |$R_V$| |$=$| 4) ranging from |$-0.5$| < |$E(B-V)$| < 2.5,¹¹ and also the normalization, avoiding emission line regions. The resulting |$E(B-V)$| values for each QSO is displayed in Tables 1 and 3. We note that using the standard extinction law of the Small Magellanic Cloud (SMC; Prevot et al. 1984) results in similar |$E(B-V)$| values, with a median difference of 0.015 mag and 0.008 mag for the cQSOs and rQSOs, respectively. The maximum difference in fitted |$E(B-V)$| between the two curves is 0.04 mag. Due to the fitting often underestimating the error, an additional 10 per cent error on the |$E(B-V)$| has been included.

The top left panel of Fig. 3 displays the overall distribution of |$E(B-V)$| values, demonstrating a small overlap in the rQSOs and cQSOs at an |$E(B-V)$| < 0.1 mag (see Fawcett et al. 2023 for a comparison between optical colour and dust extinction) and a tail towards dustier QSOs at an |$E(B-V)$| > 0.1 mag. Therefore, despite the initial sample selection based on |$g-i$| colour, for the rest of the paper we instead focus on the more physical quantity of dust extinction, displaying both the |$E(B-V)$| and |$R_V$|-dependent |$A_V$| values on the relevant figures. We also compare the ‘dusty’ QSOs to the ‘non-dusty’ QSOs by applying a separation at an |$E(B-V)$| |$=$| 0.1 mag (⁠|$A_V$| |$=$| 0.4 mag). Due to previous work that found no significant differences in the accretion properties between SDSS-selected red and blue QSOs, we do not expect varying accretion rates across the sample to significantly affect the results (Fawcett et al. 2022; see also Fig. 8).

3 RESULTS

With the combination of dedicated uGMRT observations and archival radio data, we aim to explore the differences (if any) between the radio SEDs of SDSS QSOs as a function of obscuration. In Section 3.1 we first present the uGMRT Band-3–4 radio spectral slopes and in Section 3.2 we describe the basic results from the SED fitting, exploring trends with dust extinction and radio-loudness.

3.1 uGMRT spectral slopes

To characterize the radio spectral slopes of our sample to first order, we can calculate the 2-point radio spectral slope utilizing only the uGMRT Band-3 (400 MHz) and Band-4 (650 MHz) data. This provides an empirical measurement without any model assumptions. Additionally, since these data were taken during the same week and from the same instrument, any variability issues will be minimized. Furthermore, comparing the radio spectral slopes obtained with only two data points to our fully characterized SEDs (Section 3.2), we can determine what information is gained from the inclusion of additional radio data.

The left and right panels of Fig. 3 display the uGMRT spectral slope versus |$E(B-V)$| (see Section 2.7) and radio-loudness (⁠|$\mathcal {R}$|⁠; see Section 2.5), respectively. There appears to be a tentative negative trend between the spectral slope and the amount of dust extinction (i.e. the redder the QSO, the steeper the radio spectral slope); the Spearman’s rank correlation coefficient reveals a weak trend (⁠|$\rho$| |$=$| |$-0.3$|⁠, p |$=$| 0.07). Splitting the sample by dust extinction, we find that the average spectral slope for the dusty QSOs [|$E(B-V)$| > 0.1 mag; |$A_V$| > 0.4 mag] is much steeper than that for the non-dusty QSOs [|$E(B-V)$| < 0.1 mag], with an |$\alpha$| |$=$| |$-0.71$| and |$\alpha$| |$=$| |$-0.05$|⁠, respectively.¹² We also find that 6/25 (24|$\pm$|6 per cent) of the non-dusty QSOs are inverted (⁠|$\alpha$| > 0.5) compared to zero dusty QSOs.¹³ A positive trend (⁠|$\rho$| |$=$| 0.38, p |$=$| 0.02) is found between the spectral slope and radio-loudness (i.e. more radio-quiet QSOs appear to be slightly steeper). Splitting the sample by radio-loudness, we find that the average radio spectral slope for the radio-quiet/intermediate QSOs (⁠|$\mathcal {R}$| < |$-3.7$|⁠) is much steeper than that for the extreme radio-loud QSOs (⁠|$\mathcal {R}$| > |$-3.7$|⁠; |$\alpha$| |$=$| |$-0.70$| and |$\alpha$| |$=$| |$-0.08$|⁠, respectively).

Utilizing only the uGMRT data we have found trends between the steepness of the radio spectral slope with dust extinction and radio-loudness. Therefore, this two-point fitting shows some initial tentative differences between the underlying shape of the radio SEDs in the QSOs as a function of dust extinction. In the following sections, we include the additional archival data (see Section 2.3) in order to model and characterize the full radio SED from 0.144 MHz–3 GHz.

3.2 Radio SEDs of red and blue QSOs

We fit the radio SEDs of the QSOs with a PL model (equation 2), a broken PL (BPL) model (equation 3), and a curve model (equation 4) in order to characterize the radio spectral shape, following the methodology described in Section 2.6. Additionally, each QSO was assigned a subclassification depending on either the value of |$\alpha$| in the case of a PL or BPL model best fit, or whether there was a spectral turnover in the model within the frequency range of the radio data utilized (see Section 2.6 for subclassification definitions). Example radio SED fits are display in Fig. 2, with the full sample displayed in the online Supplementary material. A breakdown of the fitting parameters for each QSO is displayed in Table 3.

We found that the majority of QSOs are best fit by either the PL or BPL model (23/38; 61|$\pm$|8 per cent¹⁴; 15 PL and 8 BPL). Out of these PL/BPL sources, 9/23 (39|$\pm$|9 per cent) have steep (⁠|$\alpha _{\rm PL}$| < |$-0.5$|⁠) radio spectral indices, indicating optically thin and evolved synchrotron emission from either radio jets or shocks caused by outflows and/or jets interacting with the surrounding ISM (Faucher-Giguère & Quataert 2012; Nims et al. 2015). The other PL/BPL sources 14/23 (61|$\pm$|11 per cent) have a flat radio spectrum (⁠|$-0.5$| < |$\alpha _{\rm PL}$| < 0.5), indicating either a self-absorbed radio nucleus or an unresolved jet-base, both on very compact scales (Rybicki & Lightman 1979; Eckart et al. 1986). No QSOs displayed a rare inverted PL spectrum.¹⁵ Of the sources not well-described by a PL or BPL model, 15/38 (39|$\pm$|7 per cent), were best fit with the curve model, with 14 displaying a clear peak within the frequency of data utilized in the fitting. A peak in the radio SED is expected to be caused by either synchrotron self-absorption (e.g. Snellen et al. 2000) or free–free absorption by ionized gas surrounding the radio emission (e.g. Stawarz et al. 2008).

We find that dusty QSOs [those with |$E(B-V)$| > 0.1 mag] are more likely to exhibit steep spectral slopes (⁠|$\alpha$| < |$-0.5$|⁠) than the non-dusty QSOs [|$E(B-V)$| < 0.1 mag], with 6/13 (46|$\pm$|12 per cent) and 3/25 (12|$\pm$|4 per cent), respectively. On the other hand, a higher fraction of the non-dusty sources are curved sources compared to the dusty QSOs; 12/25 (48|$\pm$|9 per cent) and 3/13 (23|$\pm$|8 per cent), respectively. There are no significant differences in the fraction of sources that are flat; 4/13 (31|$\pm$|10 per cent) and 10/25 (40|$\pm$|9 per cent) for the dusty and non-dusty QSOs, respectively. This could suggest that the radio emission in non-dusty QSOs is more likely to be core dominated, potentially due to compact jets, whereas the radio emission in dusty QSOs is more likely to be extended, either due to outflow-driven shocks and/or more extended jets. We explore this further in Section 4. Isolating the eight BPL sources, which indicate high frequency synchrotron losses likely due to evolved radio emission, we find a higher fraction of the dusty QSOs are best fit by a BPL compared to the non-dusty QSOs; 5/13 (38|$\pm$|11 per cent) and 3/25 (12|$\pm$|4 per cent), respectively. All five of the dusty QSOs best fit with a BPL model are also extended in the |$0{_{.}^{\prime\prime}} 2$| e-MERLIN 1.4 GHz images, which is expected for an evolved, optically thin synchrotron component (e.g. Laor et al. 2019). The other extended e-MERLIN source is non-dusty [|$E(B-V)$| |$=$| 0.04 mag] and is best fit by a flat PL.

Exploring the fitting results further, Fig. 4 displays radio-loudness versus |$E(B-V)$| for the QSOs, highlighting the sources best fit by either the PL/BPL or curve model. Splitting the samples at a radio-loudness value of |$\mathcal {R}$| |$=$| |$-3.7$| (motivated by K19, who found red QSOs showed differences in their radio properties at a |$\mathcal {R}$| < |$-3.7$|⁠) we find that the majority of radio-quiet/radio-intermediate QSOs (⁠|$\mathcal {R}$| < |$-3.7$|⁠) are best fit by the PL/BPL models (14/17; 82|$\pm$|13 per cent) compared to the curve model (3/17; 18|$\pm$|6 per cent). For the extreme radio-loud QSOs (⁠|$\mathcal {R}$| > |$-3.7$|⁠), there is an even split between QSOs best fit by the PL/BPL models (9/21; 43|$\pm$|10 per cent) or the curve model (12/21; 57|$\pm$|11 per cent), although we find 80|$\pm$|14 per cent of the curved sources lie in this region. Overall, this implies that the curved radio SEDs are more likely to be associated with radio-loud QSOs.

Splitting by dust extinction [|$E(B-V)$| |$=$| 0.1 mag], we find a higher fraction of dusty sources are best fit by the PL or BPL model (10/13; 77|$\pm$|15 per cent) compared to the curve model (3/13; 23|$\pm$|7 per cent). We find no significant differences between the fraction of the non-dusty sources best fit by either the PL/BPL models or the curve model; 13/25 (52|$\pm$|10 per cent) and 12/25 (48|$\pm$|10 per cent), respectively. Therefore, radio-loudness appears to be the most important parameter in determining whether the radio SED of a QSO is peaked or not, with more radio-quiet and dusty QSOs preferring a PL/BPL fit.

Fig. 5 displays the radio spectral index from the PL or BPL model versus the measured dust extinction (for the BPL model, the |$\alpha$| adopted is |$\alpha _{\rm high}$| if the majority of data points have a |$\nu$| > |$\nu _{\rm break}$| and |$\alpha _{\rm low}$| vice versa). We removed the three PL/BPL QSOs with unconstrained parameters due to lack of degrees of freedom (see Section 2.6 and Table 3). The horizontal dashed line indicates the boundary in |$\alpha _{\rm PL}$| between the steep and flat sub-classifications. The steep spectrum PL/BPL sources are, on average, considerably redder compared to the flat sources [median |$E(B-V)$| |$=$| 0.19 and |$E(B-V)$| |$=$| 0.05 mag for steep and flat, respectively]. Furthermore, applying the Spearman’s rank correlation coefficient reveals a tentative negative trend between |$\alpha _{\rm PL}$| and |$E(B-V)$| (⁠|$\rho$| |$=$| |$-0.37$|⁠, p |$=$| 0.11; also found in Fig. 3). This suggests that there is a connection between the steepness of the radio spectral slope and the dust extinction in QSOs. Removing the QSO with the highest dust extinction [|$E(B-V)$| |$\sim$| 0.6 mag] reduces the significance of the trend (⁠|$\rho$| |$=$| |$-0.27$|⁠, p |$=$| 0.26), demonstrating the need for a larger sample of QSOs with high values of dust extinction to robustly test this result. However, even after removing the extreme dusty QSO, the median dust extinction for the steep sources remains considerably higher than that of the flat sources [|$E(B-V)$| |$=$| 0.17 and |$E(B-V)$| |$=$| 0.05 mag for steep and flat, respectively]. We note that the non-dusty inverted sources found in Fig. 3 are no longer present after the inclusion of the archival data due to these sources now best fit by the curve model. The majority of the e-MERLIN extended sources tend to be dusty and steep, suggesting that the tentative connection between |$E(B-V)$| and |$\alpha _{\rm PL}$| might be driven by small-scale extended radio emission (i.e. |$\lesssim$|2 kpc). We discuss the potential mechanisms behind this connection in Sections 4.2 and 4.3.

4 DISCUSSION

In this paper, we have used a combination of dedicated uGMRT observations (Band-3 and 4) and archival radio data to explore the radio SEDs of a sample of QSOs across a range of dust extinctions, with additional size and morphological information from previous |$0{_{.}^{\prime\prime}} 2$| e-MERLIN 1.4 GHz imaging (presented in R21). On the basis of our analyses, we have found evidence that dusty QSOs are more likely to exhibit a steep synchrotron radio spectrum compared to non-dusty QSOs, which are more likely to be peaked. In the following sections, we discuss this key result (Section 4.1), what is the most likely origin of the radio emission in dusty QSOs and whether this is different to typical blue QSOs (Sections 4.2 and 4.3), and explore any evidence that dusty QSOs are in an earlier evolutionary phase (Section 4.4).

4.1 Steep radio power laws are more common in dusty QSOs

Splitting the QSOs into two bins of reddening [boundary: |$E(B-V)$| |$=$| 0.1 mag], we found that a higher fraction of dusty QSOs are steep compared to the non-dusty QSOs; 6/13 (46|$\pm$|12 per cent) and 3/25, 12|$\pm$|4 per cent, respectively. Furthermore, comparing the steepness of the radio spectral slope to the dust extinction, we found a weak but significant correlation (⁠|$\rho$| |$=$| |$-0.37$|⁠, p |$=$| 0.11), where the dusty QSOs are more likely to have steep spectral slopes (Fig. 5), also confirmed from only exploring the uGMRT Band-3–4 radio spectral slopes (Fig. 3). We also found a much higher fraction of radio-quiet/radio-intermediate sources (⁠|$\mathcal {R}$| < |$-3.7$|⁠) are fit by the PL/BPL models compared to the curve model (82|$\pm$|13 per cent and 18|$\pm$|6 per cent, respectively; Fig. 4). This suggests a causal connection between the radio spectral index and the dust extinction in QSOs, whereby dusty QSOs are more likely to be radio-quiet/intermediate and display a steep spectrum, which could arise from the interaction between winds and/or jets with the surrounding dusty ISM. These results are consistent with previous work, which found that the optical spectra of radio-quiet QSOs with steep radio spectral indices are more likely to display significant reddening (Baker & Hunstead 1995; Gaskell et al. 2004). We found a higher fraction of the non-dusty QSOs [|$E(B-V)$| < 0.1 mag] had peaked radio spectra compared to the dusty QSOs; 12/25 (48|$\pm$|9 per cent) and 3/13 (23|$\pm$|8 per cent), respectively, which indicates a young or frustrated radio jet. Additionally, the majority of the peaked sources are extremely radio-loud (⁠|$\mathcal {R}$| > |$-3.7$|⁠; Fig. 4; 80|$\pm$|14 per cent). We found no difference in the fraction of flat spectrum sources for the dusty and non-dusty QSOs.

Overall, these results suggest that the radio emission in non-dusty blue QSOs is more likely due to radio-loud compact jets. On the other hand, dusty QSOs are more likely to host outflows which shock the surrounding dusty medium, resulting in more radio-quiet, small-scale extended radio emission. However, both radio jets and shocks due to winds and/or jets interacting with the ISM can cause extended radio emission with steep spectral slopes. Therefore, it is still not clear whether jets or winds (or both) are the origin of the radio emission in radio-quiet/intermediate dusty QSOs; we discuss this further in Sections 4.2 and 4.3. Although radio emission from star-formation can also produce a steep spectral slope (⁠|$\alpha$| < |$-0.5$|⁠; e.g. Kimball et al. 2011; Condon et al. 2013; Calistro Rivera et al. 2017), we do not discuss this in the following sections due to the significant amount of evidence that suggests star formation is not responsible for the enhanced radio emission observed in dusty QSOs in this luminosity regime, i.e. |$L_{\rm 1.4\, GHz}\approx 10^{25-26}$| W Hz|$^{-1}$| (Fawcett et al. 2020; Rosario et al. 2020; Calistro Rivera et al. 2021, 2024; Yue et al. 2024).

4.2 Jets as the origin of the radio emission

Traditionally, radio jets have been associated with the most radio-loud QSOs (e.g. Kellermann et al. 1989). However, recent work has shown that radio jets may be more ubiquitous in QSOs, with observations of low-powered radio jets in traditionally radio-quiet sources, often associated with driving multiphase outflows (Jarvis et al. 2019; Macfarlane et al. 2021; Venturi et al. 2021; Girdhar et al. 2022; Calistro Rivera et al. 2024; Yue et al. 2024). In our sample, the QSOs with a peaked radio SED tend to be radio-loud (Fig. 4) and do not follow a strong trend with dust extinction compared to the PL/BPL sources. This could be evidence that the origin of the radio emission in these QSOs is small-scale jets, since a peaked spectrum is often associated with either a jet ‘frustrated’ by the ISM or a young jet (O’Dea et al. 1991). One way to test this theory is to compare the peaked QSOs to known jetted objects, such as GPS and CSS sources (O’Dea & Baum 1997). GPS and CSS sources are also characterized by a peaked radio SED and have compact radio emission which is thought to be jet-dominated; whether the radio jets are young or frustrated is still debated (Breugel et al. 1984; O’Dea et al. 1991; Bicknell et al. 2018).

To test whether our QSOs are consistent with GPS/CSS-like objects, we plotted the full QSO sample on the radio luminosity versus linear size plot (Fig. 6), utilizing the size constraints from e-MERLIN. We can then compare how the radio sizes compare to common radio populations: GPS, CSS, and Fanaroff–Riley Class I (FRI; brighter radio core compared to the lobes) or Class II (FRII; brighter radio lobes compared to the core; Fanaroff & Riley 1974). All of the e-MERLIN extended QSOs have radio luminosities and sizes consistent with CSS sources. However, due to the upper limits for the e-MERLIN unresolved sources, we are unable to determine if these sources are more consistent with the more compact GPS sources.

Although the QSOs in our sample have similar size constraints and radio luminosities to GPS/CSS sources, we require additional information in order to robustly determine if they are consistent with jetted sources. Another parameter we can explore is the peak frequency of the spectral turnover from the radio SED fits. For GPS/CSS sources, the frequency of the spectral turnover has been found to strongly anticorrelate with the projected linear size of the radio emission, as expected for synchrotron self-absorption (Fanti et al. 1990). For example, O’Dea & Baum (1997) found the simple relationship:

(or |$\nu _{\rm peak}$| |$\propto$| |$D^{-0.65}$|⁠) where D is the linear size of the radio emission. In Fig. 7 we plot the peak frequency versus projected linear size obtained from the e-MERLIN 1.4 GHz images (⁠|$0{_{.}^{\prime\prime}} 2$|⁠; |$\sim$| kpc-scale). We also plot on the steep PL/BPL sources since these may have a radio SED that peaks beyond our frequency range (i.e. |$\nu _{\rm peak}$| < 144 MHz, if LoTSS data are available). None of the 15 QSOs that were best fit by the curve model are extended in e-MERLIN and therefore have size upper limits. Out of the nine QSOs best fit with a steep PL or BPL, only four have constrained e-MERLIN sizes. The other QSOs were unresolved at the resolution of e-MERLIN and therefore are plotted as upper limits for the projected linear size (see Section 2.4.1 for more details on calculating the e-MERLIN sizes). For the QSOs that did not show a clear turnover within the frequency range of the radio data utilized, an upper or lower limit on the peak frequency corresponding to either the maximum or minimum frequency data point is plotted, respectively, depending on whether the overall shape of the spectrum tends to be steep or inverted. We find a few sources, both PL/BPL and curved, which are on the edge of the GPS/CSS distribution and therefore are potentially inconsistent with GPS/CSS sources, although we cannot robustly conclude this due to upper limits. Most of the other QSOs best fit by the curve model (other than the one QSO, 1103+5849, that did not display a peak in the SED and therefore also has a limit on the peak frequency) appear to be consistent with GPS/CSS-like sources within their upper limits.

To robustly conclude or rule out whether any of these sources are GPS/CSS-like we require either: (1) higher spatial resolution radio data to constrain the projected linear size, and/or (2) lower/higher frequency radio data to either constrain a turnover or, if there is still no turnover present in the spectrum, push the upper limit for the steep/inverted sources further away from the relation. For example, two of the steep QSOs do not have LoTSS data and so their upper limit is fixed at the uGMRT Band-3 data point (400 MHz). In a future study, we will use |$0{_{.}^{\prime\prime}} 05$| e-MERLIN C-band data (5 GHz; PI: V. Fawcett; Project ID: CY18008) for the same sample of QSOs studied in this paper to improve the size constraints by a factor of four, which will help us to determine whether these QSOs are GPS/CSS-like objects.

In Fig. 4 we found that the peaked QSOs, which tend to be non-dusty, are mostly radio-loud and do not follow any trends with |$E(B-V)$|⁠. This could suggest that these sources are similar to GPS/CSS sources, powered by small-scale radio-loud jets. On the other hand, it is not clear whether the PL/BPL sources (in particular, the steep sources which tend to be more dusty) are also powered by compact radio jets. Since these sources tend to be more radio-quiet and display a correlation with the amount of dust extinction, it could be possible that the predominant radio emission mechanism in these sources is shocks due to outflows, which may be driven by radio jets or accretion disc winds (e.g. Haidar et al. 2024). This is in agreement with Calistro Rivera et al. (2024), who found an enhancement in the [O iii] outflow velocities for radio-intermediate red QSOs, but not in the radio-faint regime, suggesting that dust was the key to producing large-scale ionized outflows and radio-emitting shocks. Additionally, although larger C iv outflows were found in the red QSOs compared to the blue, no differences in the C iv velocity distribution were found as a function of radio-loudness. Calistro Rivera et al. (2024) concluded that the mechanism producing the enhanced radio emission in red QSOs is not directly associated with the accretion disc, but on circumnuclear narrow line region scales such as jet-driven outflows, although this does not rule out wind-driven outflows (see Section 4.3). Simulations also predict that the interplay between compact radio jets and circumnuclear environment is expected to produce both outflows and shocks (Bicknell et al. 2018; Mukherjee et al. 2018; Young et al. 2024). This is also in agreement with recent work that has found an intrinsic, continuous connection between the amount of dust extinction in QSOs and the radio detection fraction (Fawcett et al. 2023; Calistro Rivera et al. 2024; Petley et al. 2024). These results could suggest that there are not distinct mechanisms producing the radio emission in dusty and non-dusty QSOs, but instead a continuous evolution such as jets of different sizes and/or ages (Hardcastle & Croston 2020). However, this scenario is also consistent with wind-driven shocks driving the radio emission in these QSOs, which we discuss in the following section.

4.3 Winds as the origin of the radio emission

Shocks from radiatively driven winds due to AGN radiation on dust (e.g. Ishibashi & Fabian 2015; Ishibashi, Fabian & Maiolino 2018; Leftley et al. 2019; Venanzi, Hönig & Williamson 2020; Arakawa et al. 2022) has been suggested by multiple authors to explain the origin of compact radio emission (Nims et al. 2015; Hönig 2019; Calistro Rivera et al. 2021; Yamada et al. 2024). For example, exploring a sample of ‘Extremely Red Quasars’ (ERQs; Hamann et al. 2017), selected to have red optical–MIR colours, Hwang et al. (2018) similarly found that the radio spectral slopes of the ERQs were steeper compared to typical blue QSOs. They also found that the radio emission and [O iii] velocity widths were well correlated, concluding that wind-driven shocks were the mostly likely origin of the steep radio spectra in the ERQs. Another interesting QSO subpopulation, broad absorption line QSOs (BALQSOs; Weymann et al. 1991) which are selected to have broad absorption troughs bluewards of the C iv emission line, indicative of powerful outflows, have been found to display redder optical spectra compared to non-BALQSOs (Morabito et al. 2019; Petley et al. 2022). Similarly, they are more radio-detected in the radio-quiet regime compared to non-BALQSOs (Morabito et al. 2019), with radio emission that can be explained by a wind-shock model (Petley et al. 2024). The similarities between BALQSOs and dusty QSOs could suggest that winds are also more prevalent in the dusty QSO population.

Exploring the 5–8.4 GHz radio spectral slope for a sample of 25 radio-quiet Palomar–Green (PG) QSOs (Boroson & Green 1992), Laor et al. (2019) found a significant correlation between the Eddington ratio (⁠|$L/L_{\rm Edd}$|⁠) and the steepness of the radio spectral slope. Since a high |$L/L_{\rm Edd}$| is often associated with a radiatively driven wind (Baskin & Laor 2005), they suggested that this correlation is due to wind-driven shocks in radio-quiet QSOs. However, they did not find a similar correlation between |$L/L_{\rm Edd}$| and the steepness of the radio spectral slope in radio-loud QSOs, which could suggest two different mechanisms are at play for the radio-quiet and radio-loud QSOs (see also Mehdipour & Costantini 2019). Comparing the 1.4–3 GHz radio spectral slopes of a sample of MIR-selected red and blue QSOs, Glikman et al. (2022) found that the red QSOs had on average steeper radio spectral slopes (⁠|$\alpha$| |$=$| |$-0.70\pm 0.05$| and |$\alpha$| |$=$| |$-0.34\pm 0.04$| for the red and blue QSOs, respectively). Due to the significantly higher Eddington ratios in their red sample compared to the blue, they followed the same argument as in Laor et al. (2019), suggesting that the steeper slopes are due to AGN-driven winds shocking the surrounding dust and gas. On the other hand, Rosario et al. (2020) analysed the 144 MHz–1.4 GHz radio spectral slopes of a sample of SDSS selected red and blue QSOs and found no significant differences, with both populations peaking around |$\alpha$| |$=$| |$-0.7$|⁠. In Fawcett et al. (2022) we found that SDSS red QSOs have on average slightly lower |$L/L_{\rm Edd}$| values compared to blue QSOs, albeit with a small sample. However, other studies have found either no differences in the |$L/L_{\rm Edd}$| values for red and blue QSOs (Calistro Rivera et al. 2021) or much higher |$L/L_{\rm Edd}$| values for red QSOs compared to typical blue QSOs (Urrutia et al. 2012; Kim & Im 2018; Kim et al. 2024a).

To test whether accretion-driven outflows could be responsible for the relationship between the steepness of the radio spectral slope and dust extinction, we calculated |$L/L_{\rm Edd}$| for our sample utilizing the Rakshit, Stalin & Kotilainen (2020) catalogue. Bolometric luminosities (⁠|$L_{\rm bol}$|⁠) are provided in Rakshit et al. (2020), but these are inferred from the rest-frame UV-optical continuum measurements and have not been corrected for dust extinction (see e.g. Calistro Rivera et al. 2021; Kim et al. 2023). Therefore we derive |$L_{\rm bol}$| from the 6 |$\mu$|m luminosity by applying |$L_{\rm bol}$| |$=$| |$BC_{\rm 6\mu m}$| |$\times$| |$L_{\rm 6\mu m}$|⁠, where |$BC_{\rm 6\mu m}$| |$=$| 8 from Richards et al. (2006). We inferred |$L_{\rm Edd}$| from the black hole mass (M|$_{\rm BH}$|⁠), which we estimated using the viral relation. We adopted the M|$_{\rm BH}$| calibration from Shen & Liu (2012) with the broad Mg ii full width at half-maximum and |$L_{\rm 3000}$| from Rakshit et al. (2020). We first corrected the |$L_{\rm 3000}$| for dust extinction using the |$E(B-V)$| values calculated in Section 2.7. The uncertainties for |$L/L_{\rm Edd}$| are calculated from the Mg ii emission line fitting errors presented in Rakshit et al. (2020). The top panel of Fig. 8 displays the resulting values for |$L/L_{\rm Edd}$| compared to the steepness of the PL/BPL spectral slope; we find no significant correlation between the two properties (⁠|$\rho$| |$=$| 0.14, p |$=$| 0.54), although our sample includes very radio-loud sources. This suggests that accretion-driven winds are unlikely to be the main driver behind the correlation between radio spectral steepness and dust extinction. However, Temple et al. (2023) found that accretion-driven winds, measured from the blueshift of the C iv|$\lambda$|1550 emission line, were correlated with both Eddington ratio and the black hole mass. Since our Eddington ratios are inferred from the black hole mass, future work utilizing an independent measurement on Eddington ratio, such as C iv distance (Richards et al. 2011; Rankine et al. 2020),¹⁶ is required to robustly test any correlations with |$\alpha$|⁠. We find that one of the reddest QSOs (1153+5651) displays a very steep slope, is extended in e-MERLIN, and has a high Eddington ratio. For this source, it is possible that the origin of the radio emission is due to a wind-driven shock.

To further explore whether winds can explain the observed radio SEDs of our sample, we can test a popular wind-shock model presented in Nims et al. (2015). They first assume that a wind travelling at |$0.1c$| through a typical ambient medium can reach |$\gt 1$| kpc scales over the lifetime of the QSO. They then calculate the thermal and non-thermal emission that would be produced by the shock, including the associated synchrotron emission which they predict to be in the form:

Assuming an |$L_{\rm wind}$| |$=$| |$0.05 L_{\rm bol}$|⁠, we calculate the values of |$\eta _{\rm wind}$| required to produce the observed radio luminosities of our sample. This method produces an average value of |$\eta _{\rm wind}$| to be 4 per cent; this is in agreement with previous estimates based on observations of ionized outflows, which suggest an |$\eta _{\rm wind}$| |$\lesssim$| 7 per cent (Liu et al. 2013; Sun, Greene & Zakamska 2017). However, six QSOs (16|$\pm$|4 per cent), which are also the most radio-loud (⁠|$\mathcal {R}$| > |$-3.2$|⁠), require an |$\eta _{\rm wind}$| > 7 per cent, which is likely unfeasible (see bottom panel of Fig. 8). These sources therefore require an additional radio mechanism than an accretion-disc wind alone.

An additional signature for wind-driven shocks is the detection of ionized outflows, which can be traced by the broadening or blueshifting of forbidden lines. The most popular of these lines, due to being typically bright, is the [O iii]|$\lambda$|5007 emission line; due to the redshift range of our sample (1.0 < z < 1.55) we are unable to explore the [O iii] ionized outflow properties in our sample. Previous studies have found an enhancement in the strength of [O iii] outflows in reddened QSOs (e.g. Calistro Rivera et al. 2021; Stepney et al. 2024), although other studies have found no significant differences (Fawcett et al. 2022). To robustly determine whether enhanced ionized outflows are present in dusty QSOs we require NIR spectroscopy, which can observe [O iii] out to higher redshifts; for example, the future Multi-Object Optical-to-NIR Spectrograph (MOONS; Cirasuolo et al. 2020).

In conclusion, outflow-driven shocks are a likely mechanism to explain the relationship between dust extinction and radio spectral steepness for the radio-quiet/intermediate sources. A basic wind-shock model can account for the radio luminosities in the majority of these sources; however, it is not possible to completely rule out compact/low power radio jets as the outflow driving mechanism. In the following section we discuss how these results add to the overall discussion on whether dusty QSOs represent a phase in the evolution of QSOs.

4.4 Are dusty QSOs in a younger evolutionary phase?

Although we cannot conclusively say what is the mechanism driving the radio emission in dusty QSOs (i.e. winds or jets), the higher fraction of PL/BPL sources (Section 3.2) that are more likely to be radio-quiet/intermediate (Fig. 4) and the connection between the steepness of the PL slope and dust extinction (Fig. 5) suggests that outflow-driven (wind and/or jet) shocks are the most plausible origin of the radio emission in dusty QSOs. This would be consistent with the dusty ‘blow-out’ phase (e.g. Hopkins et al. 2006; Glikman et al. 2007), since the shocks will heat the surrounding dust, destroying it and eventually revealing an unobscured blue QSO.

Further evidence that dusty QSOs are in a younger evolutionary phase arises from the LoTSS extended sources. Out of the seven QSOs with extended LoTSS radio emission, five are non-dusty [|$E(B-V)$| |$\le$| 0.06 mag], with sizes consistent with FRII sources (Fig. 6). Inspecting the LoTSS images for these sources, we find that some indeed do look like FRII objects with clear extended, symmetric bright lobes (0951+5253, 1203+4510, and 1630+3847). Others are less clear, with the presence of one bright lobe but not a second symmetric lobe (1531+4528 and 1046+3427). The remaining are only slightly extended in LoTSS, with no clear visible lobes (1304+3206 and 1428+2916). By selection, all of our sources display no extended core emission in the FIRST images (two sources do display large-scale extension, see Appendix  B), suggesting this extended low frequency emission may be evidence of relic radio emission. Indeed, An & Baan (2012) suggested that there is an evolutionary track in the radio luminosity–size plane, whereby compact GPS/CSS jets evolve over time into large-scale FRII-like jets (see Fig. 1 therein). This is in agreement with non-dusty QSOs representing an older phase in QSO evolution compared to the dusty QSOs, with large, older jet-like structures that are only observable at the lower radio frequencies.

Furthermore, three non-dusty QSOs compared to only one dusty QSO originally had upturned radio SEDs. After removing the LoTSS data point, all three of these QSOs display a peaked spectrum at |$\sim$| 650–3000 MHz frequencies, consistent with a young radio source (Kukreti et al. 2023; Kukreti & Morganti 2024), and a steep radio spectrum at 144–400 MHz, consistent with relic radio emission from previous AGN episodes (also known as ‘restarted’ AGN; e.g. Baum et al. 1990; Hancock et al. 2010; Kharb et al. 2016; Silpa et al. 2020; Jurlin et al. 2024; Nair et al. 2024). This further suggests that these QSOs have already undergone a blow-out phase, with an older radio jet that has evolved (e.g. An & Baan 2012) and managed to escape the host-galaxy. The fact that a higher fraction of these QSOs are non-dusty supports this scenario, since we expect the dust responsible for the reddening in QSOs to get either destroyed or dispersed during a blow-out phase. Alternatively, it may be evidence that these QSOs are undergoing multiple cycles of activity on short time-scales (see for example Nyland et al. 2020; Harrison & Ramos Almeida 2024). However, we also find a higher fraction of the non-dusty QSOs that display a peaked spectrum (Section 3.2), which are often associated with young sources (e.g. Kukreti et al. 2023; Kukreti & Morganti 2024). This is in disagreement with the radio emission in non-dusty QSOs being older, although there might have been previous, older relic activity that has already faded and is no longer detectable with the low frequency radio data. A larger uGMRT sample and/or a wider statistical study is required to robustly tie down the fraction of QSOs with peaked radio SEDs and low frequency extended radio emission.

5 CONCLUSIONS

We have constructed 4–5 band radio SEDs from new uGMRT observations (at 400 and 650 MHz), in combination with archival radio data, to analyse and compare the radio SEDs of 38 QSOs at 1.0 < z < 1.5 (⁠|$L_{\rm 1.4\, GHz}$| |$\sim$| |$2\times 10^{25-26}$| W Hz|$^{-1}$| and |$L_{\rm 6\mu m}$| |$\sim$| |$10^{45-46}$| erg s|$^{-1}$|⁠). Each SED was fit by either a PL, BPL, or curve model over a frequency range of 0.144–3 GHz. With the additional constraints from our |$0{_{.}^{\prime\prime}} 2$| e-MERLIN 1.4 GHz imaging, we have explored trends between the radio spectral shape, dust obscuration, radio loudness, and morphology. From our analyses we found that:

• The majority of the QSOs have radio SEDs that are best fit by a standard or broken power law: overall, 61|$\pm$|8 per cent of QSOs are best fit by the PL/BPL models, with the remaining 39|$\pm$|7 per cent best fit by the the curved model. Out of the PL/BPL sources, 39|$\pm$|9 per cent are steep (⁠|$\alpha$| < |$-0.5$|⁠), indicating standard synchrotron emission, and 61|$\pm$|11 per cent are flat, which either indicates a self-absorbed core or unresolved jet-base, both on very compact scales (Section 3.2).

• There is a trend between dust obscuration and the steepness of the radio spectral slope: splitting the QSOs best fit by the PL/BPL models into two bins of dust extinction [boundary: |$E(B-V)$| |$=$| 0.1 mag], we found that a higher fraction of the dusty QSOs are steep (46|$\pm$|12 per cent compared to 12|$\pm$|4 per cent; Section 3.2). Furthermore, we found a correlation between the amount of dust extinction and the steepness of the spectral slope, which is marginally statistically significant (Fig. 5) and confirms the result from exploring only the uGMRT Band-3–4 radio spectral slopes (Fig. 3). This is in conjunction with the previously known result that the dustier sources are also more likely to have extended 1.4 GHz radio structures on a few kpc scales (R21).

• A higher fraction of the non-dusty QSOs display a peaked spectrum, which are consistent with radio-loud jetted systems: 48|$\pm$|9 per cent of the non-dusty QSOs (⁠|$E(B-V)$| < 0.1 mag) are peaked compared to the dusty QSOs (23|$\pm$|8 per cent). A peak in the radio spectral shape is associated with either a frustrated or young, compact jet. This is also in agreement with the fact that the majority of the curved sources are extremely radio-loud (⁠|$\mathcal {R}$| > |$-3.7$|⁠; Fig. 4). Plotting the turnover frequency versus the linear projected size obtained from the e-MERLIN imaging, we found that the majority of the QSOs best fit by the curve model are broadly consistent with GPS/CSS sources (Figs 6 and 7). However, the majority of our sample are unresolved in the e-MERLIN imaging and therefore only have upper limits on the projected linear sizes.

• The radio emission in dusty QSOs is likely due to jet- or wind- driven shocks: we found that the radio emission calculated from a wind-shock model is unable to recreate the radio luminosities observed in 39|$\pm$|7 per cent of our sample, which represent the most radio-loud objects (Section 4.3). One of the reddest QSOs that displays a very steep radio spectral slope also has a very high Eddington ratio, which could suggest that accretion-driven winds are responsible for the radio emission in this object (Fig. 8). However, for the other QSOs it is currently unclear whether jets or winds are responsible for driving the radio emission. Due to the connection between the steepness of the radio spectral slope, dust extinction, and radio loudness, with dusty objects tending to be more radio-quiet and display steeper spectral slopes, we argue that the radio emission in dusty QSOs is likely dominated by outflow-driven shocks (either launched by winds or low-powered jets) in the surrounding ISM.

Overall, our results are consistent with previous studies which suggest an intrinsic connection between the amount of dust in a QSO and the production of radio emission (Fawcett et al. 2023). This connection is likely due to outflow-driven shocks in the surrounding ISM and appears to be more important in the radio-quiet/intermediate regime, where powerful jets are unlikely to be dominating (Fig. 8). This is consistent with dusty QSOs residing in a younger evolutionary phase compared to bluer QSOs, whereby the shocks destroy/disperse the surrounding dust/gas, eventually revealing a typical unobscured blue QSO. Furthermore, we find a higher fraction of non-dusty QSOs with extended LoTSS emission, which could be due to relic radio emission from a previous AGN episode. Therefore, this may be an indication that non-dusty QSOs are on average older than dusty QSOs, although a larger sample is needed to robustly test this (Section 4.4).

Ultimately, in order to confirm whether outflow-driven shocks are responsible for the radio emission in dusty QSOs, we need to analyse the radio SEDs of a larger sample of QSOs, pushing to redder and more radio-quiet objects to robustly test the observed trends. Furthermore, future high-resolution e-MERLIN C-band (5 GHz) imaging (⁠|$0{_{.}^{\prime\prime}} 05$|⁠; |$\sim$| 500 pc scales) will help us to constrain the radio morphology of these QSOs, revealing if they display jet-like structures and comparing their sizes to GPS/CSS sources. Complementary statistical studies exploring the radio spectral slopes of red and blue QSOs (e.g. Sargent et al., in preparation) will also aid our understanding of how the dust and radio properties in QSOs connect.

ACKNOWLEDGEMENTS

We would like to thank the anonymous referee for their constructive comments.

VAF and CMH acknowledge funding from an United Kingdom Research and Innovation grant (code: MR/V022830/1). DMA thanks the Science Technology Facilities Council (STFC) for support from the Durham consolidated grant (ST/T000244/1). PK, JB, and SG acknowledge the support of the Department of Atomic Energy, Government of India, under the project 12-R&D-TFR-5.02-0700. JP acknowledges support from the CAS–NWO programme for radio astronomy with project number 629.001.024, which is financed by the NWO. JP acknowledges support from the ERC Starting Grant ClusterWeb 804208. SS acknowledges financial support from Millennium Nucleus NCN23_002 (TITANs) and Comité Mixto ESO-Chile. CS acknowledges STFC for studentship support (ST/Y509346/1).

We thank the staff of the GMRT that made these observations possible. GMRT is run by the National Centre for Radio Astrophysics of the Tata Institute of Fundamental Research.

DATA AVAILABILITY

The uGMRT data underlying this article are available via the GMRT online archive facility: https://naps.ncra.tifr.res.in/goa/data/search. The uGMRT Band-3 and Band-4 radio images utilized in this paper can be found online at https://doi.org/10.25405/data.ncl.27854934.v2. The archival radio data are publicly available online (see Section 2.3).

An electronic table containing the extracted uGMRT fluxes, uGMRT image details, SED classification and best-fitting parameters, and the SDSS identifier for the red and blue QSO samples can be found in the online Supplementary material.

Footnotes

REFERENCES

APPENDIX A: RADIO VARIABILITY

A potential bias introduced when combining multiple archival radio surveys is variability. It has been long established that QSOs vary across multiple frequencies, including radio, on day-to-year time-scales (Wagner & Witzel 1995; Barvainis et al. 2005). Therefore, utilizing different radio data that span a large time-scale could result in skewed radio SEDs that do not capture the true nature of the source.

Although it is not possible to have a complete understanding of the variability of these sources due to the wide range of time-scales, we can compare the VLASS Epoch 1 and 2 data to gain some insight into short-term variability. The VLASS Epoch 1 and 2 catalogues contain two epochs of data that were taken |$\sim$| 2 yr apart. Matching to both catalogues (⁠|$0{_{.}^{\prime\prime}} 5$| matching radius for both) resulted in all 38 QSOs detected. Fig. A1 displays a comparison between the VLASS Epoch 1 and 2 integrated flux densities for the red and blue QSOs. Overall, there is little variability, with a maximum offset between the VLASS flux densities of both the red and blue QSOs of |$\sim$| |$1.6\times$|⁠. Similarly, Sargent et al. (in preparation) explored the variability of the in-band VLASS spectral slopes versus the FIRST–VLASS spectral slopes of a larger sample of SDSS red and blue QSOs and found no significant differences. Therefore, it is unlikely that radio variability will strongly affect the interpretation of our SED fitting overall; however, the SED fits for some individual sources may be impacted by variability. Interestingly, the source that have a VLASS flux much higher than the fitted radio SED model (1251+4317; see Section 3.2) does not show signs of significant variability based on the VLASS Epoch 1 to Epoch 2 integrated fluxes.

Another indication that there is no strong radio variability in our sources come from inspecting the lower resolution archival radio data, displayed on SEDs (e.g. Fig. 2), but not used in the fitting (see Table 2). We find that despite the |$\sim$| 25 yr timespan of the various data, the majority of the data points from these lower resolution surveys are consistent with the model fitted to the uGMRT, FIRST, LoTSS, and VLASS data, and also cover the same frequency range. The only exceptions are for the upturned sources, where the lower resolution, lower frequency radio data picks up the relic radio emission, as expected (e.g. QSO 1046+3427, which displays a much higher WENNS |$54\, \rm arcsec$| data point). For these sources, the lower frequency LoTSS point is removed and the SED refitted.

APPENDIX B: LARGE-SCALE RADIO EMISSION

By selection, the QSOs in our sample were chosen to have no visually extended radio emission at 1.4 GHz in the FIRST images (see Section 2.1). To confirm this and to justify our use of integrated fluxes rather than peak fluxes in our SED fitting, in Fig. B1 we plot the peak versus integrated flux density from FIRST. We find a very tight correlation which is expected for unresolved sources.

However, from inspecting the radio images of our sample we found seven sources with extended radio emission in LoTSS that appear compact in FIRST (see Section 4.4). To understand whether this is relic radio emission, which is not be present in the higher frequency data, or missed extended emission we visually inspected 3 arcmin cutouts for all the QSOs in all the bands utilized in the SED fitting. We found two sources (1046+3427 and 1057+3119) that displayed extended radio emission in all bands that appear as multiple entries in the FIRST and VLASS catalogues (Fig. B2). The extended LoTSS radio emission in QSO 1046+3427 was treated as associated in the Hardcastle et al. (2023) catalogue, and therefore originally had an upturned radio SED. The LoTSS radio emission for QSO 1057+3119, on the other hand, appears as three entries in the catalogue, and therefore the radio SED is flat due to the core emission. The other five QSOs with extended LoTSS emission either show extension in multiple bands, suggesting diffuse structures that increase in brightness with decreasing frequency, or display extension in only LoTSS and faintly in uGMRT Band-3, which could be evidence of relic emission.

APPENDIX C: COMPARISON TO THE RADIO SPECTRAL INDEX QUADRANT

The majority of radio spectral studies are limited to two or three data points, due to the large publicly available radio catalogues. A popular approach to classify the shape of the radio spectrum, rather than just obtaining a two-point spectral index, is to use the radio spectral index quadrant plot (e.g. Mahony et al. 2016; Patil et al. 2022; Kukreti et al. 2023; Sinha, Mangla & Datta 2023; Fig. C1). Typically, a higher frequency radio spectral slope (⁠|$\alpha _{\rm high}$|⁠) is plotted against a lower frequency radio spectral slope (⁠|$\alpha _{\rm low}$|⁠), with each quadrant representing a different global spectral shape. For example, both a steep |$\alpha _{\rm high}$| and |$\alpha _{\rm low}$| would give a globally steep spectral index, and likewise, if both |$\alpha _{\rm high}$| and |$\alpha _{\rm low}$| were inverted, then the global spectral index would be inverted. The top left quadrant (inverted |$\alpha _{\rm high}$| and steep |$\alpha _{\rm low}$|⁠) would produce an upturned spectrum and the bottom right quadrant (steep |$\alpha _{\rm high}$| and inverted |$\alpha _{\rm low}$|⁠) would produce a peaked spectrum. A region around |$\alpha _{\rm high}$| |$=$| 0 and |$\alpha _{\rm low}$| |$=$| 0 is usually defined to represent flat spectrum sources. To assess the reliability of this method we can compare our robust spectral classifications from the SED modelling to what we would have concluded with only three data points (e.g. LoTSS, FIRST, and VLASS).

Fig. C1 displays the radio spectral index quadrant plot for our sample, with each marker representing a different spectral classification from our SED modelling. Unfortunately, due to the lack of LoTSS coverage we can only compare the classifications of 27 QSOs in our sample. Defining our flat region to be within |$|\alpha _{\rm high}|$| < 0.5 and |$|\alpha _{\rm low}|$| < 0.5, we find the |$\sim$|18/27 of the QSOs fall in the correct quadrant (not including one borderline source). Inspecting the other nine QSO classifications we find six sources (1057+3315, 1103+5849, 1222+3723, 1323+3948, 1554+2859, and 1602+4530) are on the borderline of the classifications. The other three sources (1251+4317, 1304+3206, and 1531+4528) have one data point that is unusually higher than the other data points (either VLASS or FIRST), which may be due to variability. In conclusion, for the majority of sources the radio spectral quadrant provides an effective method for determining the global shape of the radio SED. For a small fraction of unusual or borderline cases, additional data points are required to fully describe the radio SED.