# Exploring the radio loudness of SDSS quasars with spectral stacking

M. I. Arnaudova, D. J. B. Smith, M. J. Hardcastle, S. Das, A. Drake, K. Duncan, G. Gürkan, M. Magliocchetti, L. K. Morabito, J. W. Petley, et al.

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

## To cite this version:

M. I. Arnaudova, D. J. B. Smith, M. J. Hardcastle, S. Das, A. Drake, et al. Exploring the radio loudness of SDSS quasars with spectral stacking. Monthly Notices of the Royal Astronomical Society, 2024, 528, pp.4547-4567. 10.1093/mnras/stae233. insu-04821674

HAL Id: insu-04821674

https://insu.hal.science/insu-04821674v1

Submitted on 5 Dec 2024

HAL is a multi-disciplinary open access archive for the deposit and dissemination of scientific research documents, whether they are published or not. The documents may come from teaching and research institutions in France or abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est destinée au dépôt et à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

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

Distributed under a Creative Commons Attribution 4.0 International License
# Exploring the radio loudness of SDSS quasars with spectral stacking

M. I. Arnaudova<sup>1,1*</sup> D. J. B. Smith,<sup>1</sup> M. J. Hardcastle<sup>1</sup>,<sup>1</sup> S. Das<sup>1</sup>,<sup>1</sup> A. Drake,<sup>1</sup> K. Duncan<sup>1,2</sup>

G. Gürkan,<sup>1,3,4</sup> M. Magliocchetti,<sup>5</sup> L. K. Morabito<sup>1,6,7</sup> J. W. Petley,<sup>6</sup> S. Shenoy<sup>1</sup> and C. Tasse<sup>8,9</sup>

Accepted 2024 January 16. Received 2024 January 15; in original form 2023 November 8

## ABSTRACT

We use new 144 MHz observations over 5634 deg<fnref n="2" /> from the LOFAR (Low Frequency Array) Two-metre Sky Survey (LoTSS) to compile the largest sample of uniformly selected, spectroscopically confirmed quasars from the 14th data release of the Sloan Digital Sky Survey (SDSS-DR14). Using the classical definition of radio loudness, $R = \log (L_{1.4GHz}/L_i)$, we identify 3697 radio-loud (RL) and 111 132 radio-quiet (RQ) sources at $0.6 < z < 3.4$. To study their properties, we develop a new rest-frame spectral stacking algorithm, designed with forthcoming massively multiplexed spectroscopic surveys in mind, and use it to create high signal-to-noise composite spectra of each class, matched in redshift and absolute *i*-band magnitude. We show that RL quasars have redder continuum and enhanced [O II] emission than their RQ counterparts. These results persist when additionally matching in black hole mass, suggesting that this parameter is not the defining factor in making a quasi-stellar object (QSO) RL. We find that these features are not gradually varying as a function of radio loudness, but are maintained even when probing deeper into the RQ population, indicating that a clear-cut division in radio loudness is not apparent. Upon examining the star formation rates (SFRs) inferred from the [O II] emission line, with the contribution from active galactic nucleus removed using the [Ne V] line, we find that RL quasars have a significant excess of star formation relative to RQ quasars out to $z = 1.9$ at least. Given our findings, we suggest that RL sources either preferably reside in gas-rich systems with rapidly spinning black holes, or represent an earlier obscured phase of QSO evolution.

**Key words:** techniques: spectroscopic – galaxies: active – quasars: general – radio continuum: galaxies.

## 1 INTRODUCTION

The most luminous manifestations of active galactic nuclei (AGNs) are quasi-stellar objects (QSOs), also known as quasars, whose bolometric luminosity can reach up to $10^{47-48}$ erg s⁻<fnref n="1" /> (e.g. Rakshit, Stalin & Kotilainen 2020; Shen et al. 2020). About 5 per cent–10 per cent of these optically selected sources are found to emit strongly in the radio band, likely due to the presence of relativistic jets (Urry & Padovani 1995), while the remaining 90 per cent are only weak radio sources, whose emission could be purely a result of star formation (SF, e.g. Kimball et al. 2011; Condon et al. 2013; Kellermann et al. 2016). This division into radio-loud (RL) and radio-quiet (RQ) quasars raises the question of whether these two types of objects represent physically distinct populations or different evolutionary stages of a single one.

To provide an answer to this question, studies have looked into the distribution of the radio-loudness parameter ($R$; the ratio of radio-to-optical flux density or luminosity). Some claim that this distribution is bimodal (e.g. Kellermann et al. 1989; Ivezić et al. 2002; White et al. 2007), whilst others present evidence against this bimodality (e.g. Cirasuolo et al. 2003a, b; Baloković et al. 2012; Gürkan et al. 2019; Macfarlane et al. 2021). Part of the reason for these contradictory results could be due to the definition of $R$. For example, the use of optical and radio information leads to inhomogeneous samples as a result of different selection effects. In addition, the $R$ parameter is calculated using different bands, depending on data availability, which may not give consistent results (Ivezić et al. 2002). Finally, both the optical and radio emission could be contaminated by the host galaxy, while the radio emission could be further complicated by the jet power’s dependence on the environment, time, and Doppler boosting (e.g. Liu, Jiang & Gu 2006; Gürkan et al. 2019; Radcliffe et al. 2021).

Another debated issue involves the source of radio emission in RQ quasars. In star-forming galaxies, the radio emission is

<footnote n="1"><sup>1</sup>Centre for Astrophysics Research, University of Hertfordshire, College Lane, Hatfield AL10 9AB, UK</footnote>

<footnote n="2"><sup>2</sup>Institute for Astronomy, University of Edinburgh Royal Observatory, Blackford Hill, Edinburgh EH9 3HJ, UK</footnote>

<footnote n="3"><sup>3</sup>CSIRO Space and Astronomy, ATNF, PO Box 1130, Bentley, WA 6102, Australia</footnote>

<footnote n="4"><sup>4</sup>Thüringer Landessternwarte, Sternwarte 5, D-07778 Tautenburg, Germany</footnote>

<footnote n="5"><sup>5</sup>INAF – IAPS, Via Fosso del Cavaliere 100, I-00133 Rome, Italy</footnote>

<footnote n="6"><sup>6</sup>Department of Physics, Centre for Extragalactic Astronomy, Durham University, Durham DH1 3LE, UK</footnote>

<footnote n="7"><sup>7</sup>Department of Physics, Institute for Computational Cosmology, University of Durham, South Road, Durham DH1 3LE, UK</footnote>

<footnote n="8"><sup>8</sup>GEPI and ORN, Observatoire de Paris, Université PSL, CNRS, 5 Place Jules Janssen, F-92190 Meudon, France</footnote>

<footnote n="9"><sup>9</sup>Department of Physics and Electronics, Rhodes University, PO Box 94, Grahamstown 6140, South Africa</footnote>

<footnote>*E-mail: m.i.arnaudova@gmail.com</footnote>
associated with SF through free-free emission from H II regions and synchrotron radiation from electrons accelerated to relativistic speeds by supernova remnants (Condon 1992). This leads to the question of whether SF in the host galaxy is sufficient to account for the observed radio emission from RQ quasars. Some studies find that SF is enough (e.g. Kimball et al. 2011; Condon et al. 2013), whilst others argue that it must come from AGN, in the form of small-scale jets, AGN-driven winds or disc coronal activity (e.g. Laor & Behar 2008; Zakamska et al. 2016; White et al. 2015; Symeonidis et al. 2016; White et al. 2017; Morabito et al. 2022).

The nature of jet production in RL quasars is also not clear. Following the work of Blandford & Znajek (1977), some propose that the black hole (BH) spin plays a vital role in powering radio jets (e.g. Wilson & Colbert 1995; Sikora, Stawarz & Lasota 2007; Tchekhovskoy, Narayan & McKinney 2010). However, due to the extreme difficulty in measuring this quantity, it is challenging to test this model observationally. Another potential physical parameter involved in determining the jet power and thus the distinction between RL/RQ quasars is the BH mass. While some authors find that radio loudness strongly depends on BH mass (e.g. Gu, Cao & Jiang 2001; Dunlop et al. 2003; McLure & Jarvis 2004; Best et al. 2005; Metcalf & Magliocchetti 2006; Chiaberge & Marconi 2011), others have found only a weak dependence or no dependence at all (e.g. Ho 2002; Shankar et al. 2010; Gürkan et al. 2019; Macfarlane et al. 2021).

To investigate these problems, Gürkan et al. (2019) combined a sample of optically selected quasars from the fourteenth data release of the Sloan Digital Sky Survey (SDSS DR14; Pâris et al. 2018) with radio observations from the first data release of the LOFAR (Low Frequency Array) Two-metre Sky Survey (LoTSS DR1; Shimwell et al. 2017) over the Hobby-Eberly Telescope Dark Energy Experiment spring field (HETDEX; Hill et al. 2008) and the LOFAR Herschel Astrophysical Terahertz Large Area Survey/North Galactic field (H-ATLAS/NGP; Hardcastle et al. 2016). With the high sensitivity and wide area coverage of LoTSS, the authors were able to study the dependence of the radio-loudness parameter on galaxy properties such as redshift, bolometric luminosity, radio luminosity, BH mass, and Eddington ratio. They found that quasars exhibit a wide continuum of radio properties, with no clear signatures of a bimodality. Given these results, the authors favoured the scenario where both AGN jets and SF contribute to the radio emission in quasars such that there is no RL/RQ dichotomy, but rather a smooth transition between the dominance of these two processes.

Recently, Macfarlane et al. (2021) built upon these results by developing a numerical model of the radio flux densities of quasars, in which the radio emission of every quasar consists of two components: AGNs (jets) and SF. This model, coupled with Monte Carlo simulations, allowed the authors to create quasar mock samples and compare them with observations. Their results were found to be in excellent agreement with the observed radio flux distributions of ~42 000 SDSS quasars as measured in LoTSS DR1 across several redshift and absolute *i*-band magnitude ranges. This is consistent with a model in which jet production is present in all quasars with a different powering efficiency such that it leads to a smooth transition between the RQ and RL quasar regimes.

Our work takes a different approach by developing a spectral stacking algorithm and using it with the much larger LoTSS DR2 sample. With its extensive coverage of 5634 deg², LoTSS DR2 provides a much larger observational volume, resulting in a substantial increase of RL quasars. Employing our stacking techniques on this expanded data set allows us to systematically create composite spectra for each radio class in a given parameter regime (i.e. redshift, *i*-band

luminosity, and BH mass). This approach enables us to thoroughly explore the continuum and emission-line properties of quasars, while also determining the influence of key physical parameters such as BH mass and Eddington ratio. Therefore, although radio loudness may not correspond to a physical property of QSOs, it can be useful for identifying sources with powerful jets, and high signal-to-noise (S/N) ratio spectra provide an excellent way to investigate their properties. This stacking algorithm is designed with the upcoming WEAVE-LOFAR survey (Smith et al. 2016) in mind, which will provide over a million spectra of LoTSS targets selected at 144 MHz. Such a tool is necessary as a result of the survey's radio selection criteria, which produces samples dominated by AGN and/or ongoing SF. The spectra of such sources are rich in emission lines which allow us to robustly determine the redshifts, but a continuum detection is not always available (e.g. for faint star-forming galaxies). Stacking such sources together, however, allows us to statistically detect the continuum and thus recover spectral features otherwise indistinguishable in individual detections. Furthermore, stacking sources of different demographics will enable us to create a large library of high-resolution templates that will help improve the WEAVE-LOFAR survey's redshift estimates.

The structure of this paper is as follows. Section 2 provides details of the spectroscopic and radio data used in this study, along with the methodology for the sample selection and matching process. In Section 3, we describe the spectral stacking technique employed for comparing the radio classes of QSOs. Subsequently, in Sections 4 and 5, we create high S/N composite spectra of QSOs and investigate potential factors contributing to the observed effects between them. Finally, Section 6 discusses possible explanations and Section 7 summaries our main results. Throughout this work, we use vacuum wavelengths and a flat lambda-cold dark matter cosmology with $\Omega_{\Lambda} = 0.7$, $\Omega_{M} = 0.3$, and $H_{0} = 70 \text{ km s}^{-1} \text{ Mpc}^{-1}$.

## 2 DATA

### 2.1 Sloan Digital Sky Survey

The spectroscopic data used in this work are taken from the fourteenth data release of the Sloan Digital Sky Survey Quasar Catalogue (SDSS-DR14Q), which is fully described by Pâris et al. (2018). This catalogue includes all spectroscopically confirmed quasars from SDSS-I/II, SDSS-III/BOSS, and SDSS-IV/eBOSS programmes, resulting in a sample of 526 356 objects over a region of 9376 deg² as shown in Fig. 1.

SDSS-I/II contains 79 847 quasars with *i*-band absolute magnitudes brighter than $M_i[z = 2] = -22.0$ over a wide redshift range of $0.065 < z < 5.46$. The targeting algorithm used to obtain these sources has been updated throughout the years and more information about it is provided by Richards et al. (2002) and Schneider et al. (2010). The spectra are produced by a pair of multi-object spectrographs (SDSS) that have a total of 640 fibres with an entrance diameter of 3 arcsec and an average resolving power of $\Delta\lambda/\lambda \approx 2000$ across the vacuum wavelength range from 3800 to 9200 Å.

The spectra for SDSS-III/BOSS and SDSS-IV/eBOSS, on the other hand, are produced by the upgraded pair of spectrographs used for the Baryon Oscillation Spectroscopic Survey (BOSS). These spectrographs benefit from having a total of 1000 fibres with a smaller entrance diameter of 2 arcsec, an extended vacuum wavelength coverage of 3600–10 400 Å, and a resolving power higher in the red channel and lower in the blue channel compared to the SDSS (Smee et al. 2013). These two programs use multiple target selection
<figure><img src="image_3.png" /><figcaption>**Figure 1.** Sky coverage of SDSS-DR14Q and LoTSS DR2. The grey region indicates all the sources from the SDSS-DR14Q catalogue, whereas the light blue region – those located in the LoTSS DR2 footprint, where cross-matching has been performed.</figcaption></figure>

algorithms to detect much fainter quasars ($M_i[z = 2] = -20.5$) at redshift ranges of $2.15 < z < 3.5$ and $0.9 < z < 2.2$, respectively (for further details, see Ross et al. 2012; Myers et al. 2015). As a result of these differences and the considerable size of 446 781 for the latter two programs, we use only the BOSS spectra in our analysis.

Complementary to this SDSS-DR14Q catalogue is the value-added catalogue by Rakshit et al. (2020) which includes continuum and line property measurements, including bolometric luminosity ($L_{bol}$), derived virial black hole mass ($M_{BH}$), and Eddington ratio ($R_{edd}$) estimates. This was done using the publicly available spectral fitting code PYQSOFIT (Guo, Shen & Wang 2018), which uses two independent sets of eigenspectra – pure galaxy (Yip et al. 2004a) and pure quasar (Yip et al. 2004b) – to decompose the spectrum into host galaxy and quasar contribution. This decomposition is particularly important in low-redshift quasars ($z < 0.8$), where the stellar contribution can be significant (e.g. Yue et al. 2018; Shen et al. 2019; Rakshit et al. 2020). Briefly, the continuum of the host-free spectrum is modelled by a combination of power law, Fe II and Balmer continuum models, whereas the emission lines are fitted with Gaussian distributions. The bolometric luminosity is computed from $L_{5100}$ ($z < 0.7$), $L_{3000}$ ($0.7 \le z < 1.9$), and $L_{1350}$ ($z \ge 1.9$) using the bolometric corrections from Richards et al. (2006). For the virial BH mass multiple measurements are included in the catalogue depending on the availability of strong emission lines and various calibrations. In this work, we choose to use the ‘fiducial’ estimate calculated based on the H$\beta$ line ($z < 0.8$) and the C IV line ($z \ge 1.9$) using the calibrations from Vestergaard & Peterson (2006), and the Mg II line ($0.8 \le z < 1.9$) using the calibration from Vestergaard & Osmer (2009). These measurements are also used to calculate the Eddington luminosity and thus the Eddington ratio, which is used as a proxy for the accretion rate.

## 2.2 LOFAR Two-metre Sky Survey

The radio data used in this work are taken from the LoTSS DR2 (Shimwell et al. 2022). This data release covers $5634 \text{ deg}^2$ of the northern sky (see footprint in Fig. 1) at a resolution of 6 arcsec and a median rms sensitivity of 83 $\mu$Jy beam$^{-1}$, providing over 4 million radio sources, the vast majority of which have never been detected at radio wavelengths before. In addition to the large area and high sensitivity, LoTSS benefits from using the low radio central frequency of 144 MHz which amongst other advantages reduces the effects from Doppler boosting in jetted sources.

In preparation for the upcoming WEAVE-LOFAR survey, Hardcastle et al. (2023) have created a preliminary cross-matched catalogue containing 296 921 SDSS counterparts by using a combination of statistical methods and visual classification in a similar manner to Kondapally et al. (2021). To summarize, for smaller, more compact radio sources, the likelihood ratio method is applied (e.g. Richter 1975; De Ruiter, Willis & Arp 1977; Sutherland & Saunders 1992). This statistical approach uses both colour and magnitude information as described by Williams et al. (2019) in order to maximize the number of identifications, as well as increase the robustness of the cross-matching. For larger sources with extended radio emission, the web-based interface for visual inspection – the Radio Galaxy Zoo: LOFAR, is needed. Here, the user can perform identification and association for a given radio source by making use of available multiwavelength images. Finally, for everything in between a decision tree is used to combine the two methods.

## 2.3 Final sample

Starting with a sample of 446 781 BOSS quasars from the DR14Q catalogue, we identify those located in the LoTSS DR2 area with the use of a Multi-Order Coverage map as seen in Fig. 1. This results in 265 578 objects, out of which 33 968 are part of the LoTSS DR2 optical cross-matched catalogue. Next, we remove spectra with bad plate quality and a non-zero quality flag for $M_{BH}$ and $L_{bol}$ to ensure reliability of results. Finally, to mitigate any selection bias, we use the target selection flags to discard objects targeted by SDSS solely based on their radio or X-ray emission, as well as any found to have time variability. The final sample consists of 189 680 quasars, out of which 123 742 are part of the BOSS/eBOSS homogeneously selected "CORE" sample. The CORE sample is created by using a single target selection algorithm, the Extreme Deconvolution (XDQSO; Bovy et al. 2011), which is designed to meet science goals (such as clustering and luminosity function measurements) which require a uniform selection, and so is ideal for our purposes.

To separate quasars into RL and RQ, we adopt the standard definition of the radio-loudness parameter ($R$): the logarithm of the ratio of 1.4 GHz radio luminosity ($L_{1.4GHz}$) to optical $i$-band luminosity ($L_i$), where $R = 1$ marks the boundary (e.g. Baloković et al. 2012). To calculate the $k$-corrected $L_{1.4GHz}$ for radio-detected QSOs, we use the integrated 144 MHz flux density ($S_{144MHz}$) from the LoTSS DR2 catalogue, a radio spectral index of $\alpha_{rad} = -0.7$ and the sources’ spectroscopic redshifts reported in the DR14Q catalogue. The same calculation is applied to the radio-undetected QSOs in order to obtain a $5\sigma$ upper limit for $L_{1.4GHz}$, where $S_{144MHz}$ is taken to be 5 times the local rms value taken directly from the LoTSS rms maps at the coordinates of the quasars given by SDSS DR14Q.

For the optical luminosity, we take the absolute $i$-band magnitude from the DR14Q catalogue which is $k$-corrected to $z = 2$ and use the following conversion given by Richards et al. (2006) to obtain a $k$-corrected estimate to $z = 0$:

$$M_i[z = 0] = M_i[z = 2] + 2.5(1 + \alpha_{\text{opt}}) \log(1 + z). \quad (1)$$

Here, $z = 2$ and we adopt an optical spectral index of $\alpha_{opt} = -0.5$. To obtain the $i$-band luminosity, we use the simple relationship between luminosity and magnitude, where the solar luminosity ($3.827 \times 10^{26}$ W) and solar $i$-band absolute magnitude (4.58) were used. This gives us a total of 3697 RL and 111 132 RQ quasars. The remaining sources are discarded as their $5\sigma$ upper limit puts them in the RL regime, such that we are unable to classify them with confidence.

The radio-loudness distribution of all sources, including a subset of radio detections, is presented in Fig. 2 (where we have also included a
<figure><img src="image_4.png" /><figcaption>**Figure 2.** The radio-loudness parameter for all sources (dark blue) and only radio detections (light blue). The solid black line denotes the standard division line for RL and RQ QSOs, whereas the dashed black lines are additional divisions we make to define classes as a function of *R* as denoted.</figcaption></figure>

second horizontal axis to show *R* in terms of the *L*<sub>144MHz</sub> for comparison. It is evident that in both cases the conventional demarcation line (*R* = 1) does not represent any particularly significant level for the sample, however we have repeated the analysis considering different values (as denoted in Fig. 2) finding that our results are qualitatively unchanged. Doing these tests gives additional advantages of further testing whether there is a gradual change between the two populations as found in previous studies (e.g. Gürkan et al. 2019; Macfarlane et al. 2021), as well as avoiding biases that may result from assuming a single radio spectral index.

## 2.4 The matching process

To make a robust comparison between quasar populations, we develop a method to create samples matched in 2D and 3D parameter space. For clarity, in what follows we describe the method for the RL and RQ population, where we create samples matched in *z*, *M*<sub>i</sub> and *z*, *M*<sub>i</sub>, *M*<sub>BH</sub> to use in Sections 4.1 and 4.2. However, it is important to note that this method can be adapted to various other scenarios, as demonstrated throughout this study.

The 2D matching process involves generating a 2D histogram of the RL QSOs, and identifying the number of RQ that fall within the same *z* and *M*<sub>i</sub> bin. For the one-to-one case, this entails choosing the same number of RQ counts as found for the RL 2D histogram, such that if one bin is populated by 100 RL sources then we randomly select only 100 RQ counterparts. However, as our RQ sample is considerably larger, we instead choose to normalize the RL 2D histogram by the total number of RL sources, and multiply the normalized counts by the maximum number of RQ QSOs for which we can say that they are drawn from the same *z* and *M*<sub>i</sub> distribution as the RL population. This is done by employing multiple Kolmogorov-Smirnov (K-S) tests until we accept the null hypothesis that the two samples are drawn from the same distribution if the *p*-value is greater than a significance level of *α* = 0.05. This results in a sample of 77 532 RQ QSOs matched in *z* and *M*<sub>i</sub> to our RL sample (hereafter the 2D matched or *z* − *M*<sub>i</sub> sample).

Similarly, the 3D matching process to create a sample matched in *z*, *M*<sub>i</sub>, and *M*<sub>BH</sub> (hereafter the 3D matched or *z* − *M*<sub>i</sub> − *M*<sub>BH</sub> sample) involves creating a normalized 3D histogram of the RL population, and again using multiple K-S tests to determine the maximum number

<figure><img src="image_5.png" /><figcaption>**Figure 3.** A corner plot summary of the redshift, absolute *i*-band magnitude and BH mass distributions for the RL and RQ QSOs. The density contours represent the 10 per cent, 30 per cent, 50 per cent, 70 per cent, and 90 per cent of the whole (dark blue) and 3D-matched (dashed, light blue) RQ sample, whereas the density map represents the corresponding percentages for the RL sample. The histograms share a common *y*-axis and indicate the individual *z*, *M*<sub>i</sub>, and *M*<sub>BH</sub> distribution, where the quoted *p*-values in the upper-right corner are obtained by performing a K-S test on the RL and 3D-matched RQ sample.</figcaption></figure>

of RQ sources. The effectiveness of this matching procedure can be seen in Fig. 3, where we compare the individual *z*, *M*<sub>i</sub>, and *M*<sub>BH</sub> distributions, along with contours of the 2D distributions of our 3D-matched sample of 23 716 RQ QSOs with the whole RL and RQ samples.

# 3 THE SPECTRAL STACKING METHOD

## 3.1 The stacking technique

To create high S/N composite spectra of quasars we use the following method. First, we correct for Galactic extinction in the observed frame by using the re-calibrated reddening data, *E*(*B* − *V*), from Schlegel, Finkbeiner & Davis (1998), along with the Milky Way reddening curve from Fitzpatrick (1999) for an extinction-to-reddening ratio of *R*<sub>*V*</sub> = 3.1.

Next, we shift the spectra to the rest frame by using the spectroscopic redshifts from the SDSS-DR14Q catalogue, and resample onto a common wavelength grid, predetermined by the minimum and maximum values of the rest-frame wavelengths sampled, along with the choice of sampling input in the algorithm (e.g. 1 Å pixel<sup>-1</sup>). The resampling is performed using the SPECTRES: Simple Spectral Resampling tool (Carnall 2017), where the old wavelength grid is cross-matched with the new one, such that if an old wavelength bin spreads between multiple new ones, the flux density value associated with it is distributed in the new grid based on the fractional coverage. This ensures that the flux density is conserved. To normalize the resampled spectra, we divide through each spectrum by its median, computed at the reddest possible end of the area where all spectra populate the common wavelength grid, where prominent emission lines are masked out. Doing this ensures that the overall spectral shape is preserved, in the sense that we are able to recover stellar population synthesis models once the uncertainties have been estimated and corrected for bias (see the next section for details). In
addition, by using the reddest possible common wavelength range, we minimize the effects of extinction.

Stacking the de-redshifted, resampled, and normalized spectra is generally done in the literature either by taking the weighted average, where the weights are given by the inverse variance of the spectra, or by taking the median (e.g. Rowlands et al. 2012; Zhu et al. 2015; Rigby et al. 2018; Calabrò et al. 2021). The weighted average method is found to produce a higher S/N composite spectrum, however, as explained in Calabrò et al. (2021), this is the result of a bias towards lower redshifts and/or individual spectra with higher S/N. On that account, we have chosen to build a composite spectrum by taking the median of all normalized flux density values that fall within a given wavelength bin. To calculate the uncertainty associated with this stack, we use bootstrapping to randomly resample the spectra at each wavelength bin, creating 1000 realizations with the same size as the original sample. The median method is then performed for each realization to estimate the flux density distribution at a given bin and thus determine the standard deviation.

Fig. 4 presents the result of implementing this method on a sample of RQ quasars, where the individual spectra are shown in grey, and the high-resolution stack in black. This shows how powerful spectral stacking really is in identifying spectral features otherwise undetected in individual spectra.

## 3.2 The stacking corrections

To evaluate the performance of the stacking procedure described in the previous section, we create a model galaxy spectrum by making use of the stellar population synthesis code from Bruzual & Charlot (2003), along with equivalent width measurements of nebular emission lines taken from the value added catalogue produced by the Max-Planck-Institute for Astrophysics-John Hopkins University (Brinchmann et al. 2004; Kauffmann et al. 2003; Tremonti et al. 2004). This galaxy template has a broad wavelength range extending from 90 Å to 160 µm, a resolving power similar to that of SDSS across the range from 3200 to 9500 Å and contains the main nebular emission lines targeted by both WEAVE-LOFAR and SDSS. These features enable us to simulate mock samples of galaxy spectra based on SDSS characteristics across a wide range of redshifts, providing us with a robust way of testing the spectral stacking algorithm and further demonstrating its flexibility in handling different type of spectra (see Fig. A1 for an example of stacking spectra as a function of BPT classes).

Following a range of tests including continuum recovery at low brightness and performance related to different sampling of spectra (the main motivations for implementing it for the WEAVE-LOFAR survey), we have identified two implications of the stacking method. First of all, it is not always possible for the algorithm to choose a normalization range without strong spectral features, which causes an offset at the blue end of the composite spectrum. Secondly, the spectra are combined in the rest frame, where the normalized spectra no longer have the same spectral resolution as a result of the de-redshifting process. This is not accounted for in the bootstrap method, leaving the uncertainties underestimated. To correct for these effects, we include in our stacking algorithm the following procedure.

For a given spectroscopic sample, we create a composite spectrum (hereafter the input composite spectrum) and use it as a template to simulate quasar spectra based on the characteristics of the original sample such as redshift, i-band magnitude, and inverse variance. Next, we stack the simulated spectra to produce a second spectrum (hereafter the simulated composite spectrum) and obtain the residual

in uncertainty units, which is defined as:

$$ \chi_{\text{input-sim}} = \frac{F_{\lambda, \text{input}} - F_{\lambda, \text{sim}}}{\sigma_{\lambda, \text{sim}}}, \qquad (2) $$

where $F_{\lambda, \text{input}}$ and $F_{\lambda, \text{sim}}$ are the flux density of the input and simulated composite, respectively and $\sigma_{\lambda, \text{sim}}$ is the bootstrap uncertainty for the simulated composite. Provided that the issues discussed above are not present, this quantity would be normally distributed with a mean of zero and a standard deviation of one, making it the source of our corrections. To obtain the necessary correction factors, we separate $\chi_{\text{input-sim}}$ into 30 wavelength slices (or less, depending on pixel availability) and for each one we fit a normal distribution. The best-fitting parameters are then interpolated with a spline to produce wavelength-dependent corrections as indicated in Fig. 5. The best fit of the mean is used to correct for the offset as a result of the normalization by subtracting a multiple of it and the uncertainty of the simulated composite from the flux density of the input composite. The best fit of the standard deviation, on the other hand, is used as a multiplication factor for the uncertainty of the input composite spectrum, such that resolution effects are now taken into account.

We note that this procedure does not account for different levels of extinction or diversity of samples. However, spectra of QSOs are reasonably homogeneous and the use of large statistical samples binned in a given parameter space is close to homogeneity.

## 3.3 Qualitative template comparison

To compare our spectral stacking techniques with the literature, we create RL and RQ composites spanning the broadest possible wavelength coverage. Generating such templates involves using sources at widely different redshifts (0.5 < z < 3.5), which affects our method in two ways. First of all, it has implications for the choice of sources which must be drawn from a similar z - Mᵢ space in order to reduce the impact of the Malmquist bias. Secondly, the stacking procedure needs to be performed on smaller redshift ranges as a result of the normalization process, which requires a common wavelength range free from strong spectral features. To account for this, we select all sources with -26 < Mᵢ < -24 and divide them in redshift bins of size Δz = 0.5 to create six composite spectra per quasar population using the stacking algorithm and the spectral corrections described in Sections 3.1 and 3.2. These stacks are then re-scaled to the composite at the lowest redshift bin and combined into a single template using the inverse-variance-weighted average.

The resulting RL and RQ composites can be seen in Fig. 6 where they are compared to the SDSS composite from Vanden Berk et al. (2001), the X-Shooter composite from Selsing et al. (2016), the BOSS composite from Harris et al. (2016), and to a new high-z template that we have constructed by applying our algorithm to the sample of 24 radio-bright quasars (21 out of which are classified as RL) at 4.9 < z < 6.6, discovered by Gloudemans et al. (2022). Building such a template presents new challenges, since unlike the homogeneous eBOSS spectra used for the RL/RQ templates, the high-z stack requires combining heteroscedastic data obtained from a range of facilities, including the Faint Object Camera and Spectrograph (Kashikawa et al. 2002) on the Subaru Telescope, the LRS2 instrument (Chonis et al. 2016) on the Hobby-Eberly Telescope and the Low-Resolution Imaging Spectrometer (Oke et al. 1995) on the Keck telescope. Specifically, and as discussed in Section 3.2, the simulation component of our stacking algorithm requires photometry with a bandpass that overlaps with the observed spectrum, and given that this was unavailable for the Gloudemans et al. (2022) sample, we were unable to conduct the simulations
<figure><img src="image_6.png" /><figcaption>**Figure 4.** An example composite spectrum of RQ QSOs. The average spectrum is indicated by the thick black line, while the grey region indicates the individual normalized, resampled spectra that went in to the stack. The labelled dotted blue lines indicate prominent emission lines visible in the stacked spectrum.</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>**Figure 5.** A diagnostic plot of an example composite spectrum of RQ QSOs. The upper panel shows the input (dark blue) and simulated (light blue) composite spectra, along with the number of QSOs that went into them at each rest-frame wavelength divided by the total number of QSOs (dashed green line). The second and third panels present the best-fitting parameters for the mean and standard deviation for each wavelength slice, while the lower panel presents the corrected residual in units of the propagated uncertainty.</figcaption></figure>

and it is therefore likely that the flux density has a normalization bias with uncertainties systematically underestimated. However, creating and using such a template is beneficial as this is the largest sample of RL quasars at $z > 4.9$ and its use demonstrates the flexibility of our stacking algorithm in dealing with non-SDSS spectra.

Despite the difference in selection, we find qualitatively a good agreement in the continuum shape with the various templates. The difference observed redwards of 5000 Å amongst the composites is likely caused by various degrees of host galaxy contamination. The SDSS composite from Vanden Berk et al. (2001) includes intrinsically faint sources at low redshifts, making it subject to
<figure><img src="image_8.png" /><figcaption>**Figure 6.** A comparison between various composite spectra as indicated by the legend, where the thick lines indicate the flux density, whereas the shaded area corresponds to ± the uncertainty. The upper panel presents the full-wavelength coverage, where each composite is normalized to the RQ composite at 3000 – 3600Å. The lower left panel presents a zoom-in window indicating the differences bluewards of Lyα, whereas the lower right panel – the discrepancy above 5000Å which is likely caused by host galaxy contamination.</figcaption></figure>

significant host contamination (e.g. Glikman, Helfand & White 2006; Fynbo et al. 2012). The selection by Selsing et al. (2016), on the other hand, is chosen specifically to circumvent this problem, suggesting that both the RL and RQ composite contain low levels of host galaxy contamination. This will be further investigated in Section 5.2. The disagreement bluewards of Lyα, on the other hand, can be explained by the different intergalactic medium (IGM) absorption corrections applied (or the lack thereof). We have not made any corrections, and so the intense drop in flux density for the high-z template is expected considering the redshift range of the sample. The RL composite appears to agree well with the results of Vanden Berk et al. (2001) as seen in the lower left panel of Fig. 6, where no IGM correction is applied as well. Interestingly, the RQ composite appears to be less affected by the Lyα-forest absorption as it is found to be in better agreement with the X-shooter and BOSS composite where each spectrum within the composites has been individually corrected.

## 4 COMPOSITE SPECTRA OF QUASARS

In this section, we employ our spectral stacking techniques to compare quasars as a function of radio-loudness at different redshift bins. In Sections 4.1 and 4.2, we construct high S/N composite spectra for the classically defined RL and RQ QSOs using the $z - M_i$ and $z - M_i - M_{BH}$ samples to establish whether there are any physical differences between them and investigate the potential influence of the BH mass. In Section 4.3, we re-bin our sample into four radio-loudness bins to investigate the presence of a bimodality among the QSOs or a more gradual transition between these populations.

### 4.1 High S/N comparison of RL and RQ QSOs

Starting with our $z - M_i$ sample, we use the spectral stacking algorithm described in Sections 3.1 and 3.2 to define high S/N composite spectra across ten redshift bins for each population. The redshift bins are chosen such that they are consistent with the different redshift ranges used to calculate the ‘fiducial’ $M_{BH}$ from Rakshit et al. (2020) and contain comparable numbers of RL sources.

The resulting comparisons for all 10 redshift bins are presented in Fig. 7. For each bin, we present an upper panel including the composite spectra representative of the RL (dark blue) and RQ population (light blue) and a lower panel indicating the residual in units of the propagated uncertainties (grey). We can see that there is a consistent picture emerging across all redshift bins considered – the RL QSOs appear to have on average a redder continuum with notable differences in a number of prominent emission lines.

To quantify the significance of the observed results, we perform the same stacking methods using only the RQ population in order to create Monte Carlo simulations under the null hypothesis that the RL sample is consistent with having been randomly drawn from the RQ sample. For each simulation, we take the full sample of RQ sources and randomly draw a subsample of size equal to the RL population (hereafter referred to as the RL test sample). The remaining RQ sources are then matched in $z$ and $M_i$ as described in Section 2.4 to create a RQ test sample. Finally, the stacking algorithm is applied on the two test samples to create composite spectra across the same 10 redshift bins for which we find a ‘reduced chi-squared value’ of

$$ \chi_v^2 = \frac{\sum_{i=1}^{N} \frac{(F_{RL,i} - F_{RQ,i})^2}{\sigma_{RL,i}^2 + \sigma_{RQ,i}^2}}{N}, \quad (3) $$
<figure><img src="image_9.png" /><figcaption>**Figure 7.** A comparison between the RL and RQ population across 10 redshift bins. The upper panel of each bin presents the composite spectra of RL QSOs (dark blue) and RQ QSOs (light blue), whereas the lower panel indicates the residual in units of propagated uncertainty (grey). In each panel, the redshift range is indicated in the upper centre of the plot, while the *p*-values of the null hypothesis test are presented in the legend in the upper left corner. Finally, the number of sources and median S/N is indicated in the legend in the top right corner.</figcaption></figure>

where *N* is the number of pixels in each stacked spectrum. Having done this 1000 times, we obtain a $\chi^2$ distribution under the null hypothesis for each redshift range. By fitting this distribution with its parametric form, we can estimate the probability of obtaining our observed results under the null hypothesis.

The results of the null hypothesis test for this 2D-matched sample are presented in the upper left corner of each panel in Fig. 7 and Table 1. For all redshift bins, we find *p*-values smaller than a significance level of $\alpha = 0.05$, indicating that the two samples are not drawn from the same parent population.

## 4.2 Black hole mass dependence

As discussed in the Introduction, several studies have found the BH mass to be a defining factor in the radio-loudness dichotomy, where

RL QSOs are found to harbour BHs with $M_{BH} \gtrsim 10^8$ (e.g. Dunlop et al. 2003, McLure & Jarvis 2004; Chiaberge & Marconi 2011). This, however, does not appear to be the case for our sample. In Fig. 8, we present the BH mass distributions for each radio class using the virial BH mass estimates from Rakshit et al. (2020). We can see that both the RL (dark blue) and RQ QSOs (light blue) are found to span similar $M_{BH}$ ranges, irrespective of the different BH mass calibrations. In addition, the mean BH masses do not show any systematic trend—the RL appear more massive for the $H\beta$ and Mg II, but less so for the C IV calibration (cf. McLure & Jarvis 2004). However, to determine whether this parameter plays a role in causing the differences observed in Fig. 7, we re-do the stacking procedure with the $z - M_i - M_{BH}$ sample, that is, including the $M_{BH}$ in the process, even though it reduces the sample size. We note that by matching in $M_i$ and $M_{BH}$ that we are also effectively controlling for the Eddington ratio.
<table><caption>**Table 1.** The *p*-values obtained from the null hypothesis test. The columns represent the results obtained from comparing: the RL and RQ populations with the 2D- and 3D-matched samples, the more RL bins (*R*<sub>2</sub>, *R*<sub>3</sub>, and *R*<sub>4</sub>) and the radio-quietest sample (*R*<sub>1</sub>), and the radio-detected (*R*<sub>1D</sub>) and the radio-undetected (*R*<sub>1U</sub>) sample for the radio-quietest bin defined in Section 4.3 per a given redshift bin.</caption><thead><tr><th>Redshift range</th><th>(RL - RQ)<sub>2D</sub></th><th>(RL - RQ)<sub>3D</sub></th><th>R<sub>1</sub> - R<sub>2</sub></th><th>R<sub>1</sub> - R<sub>3</sub></th><th>R<sub>1</sub> - R<sub>4</sub></th><th>R<sub>1U</sub> - R<sub>1D</sub></th></tr></thead><tbody><tr><td>0.6-0.8</td><td>0.01</td><td>0.03</td><td></td><td></td><td></td><td></td></tr><tr><td>0.8-1.1</td><td>3.39 × 10<sup>-18</sup></td><td>1.48 × 10<sup>-14</sup></td><td>5.51 × 10<sup>-41</sup></td><td>7.90 × 10<sup>-15</sup></td><td>2.26 × 10<sup>-6</sup></td><td>0.01</td></tr><tr><td>1.1-1.4</td><td>2.31 × 10<sup>-53</sup></td><td>3.014 × 10<sup>-51</sup></td><td></td><td></td><td></td><td></td></tr><tr><td>1.4-1.7</td><td>2.83 × 10<sup>-26</sup></td><td>1.88 × 10<sup>-25</sup></td><td>4.73 × 10<sup>-45</sup></td><td>1.43 × 10<sup>-18</sup></td><td>2.16 × 10<sup>-18</sup></td><td>0.31</td></tr><tr><td>1.7-1.9</td><td>2.88 × 10<sup>-45</sup></td><td>3.40 × 10<sup>-28</sup></td><td></td><td></td><td></td><td></td></tr><tr><td>1.9-2.1</td><td>7.43 × 10<sup>-26</sup></td><td>1.67 × 10<sup>-13</sup></td><td>2.38 × 10<sup>-38</sup></td><td>1.92 × 10<sup>-23</sup></td><td>1.37 × 10<sup>-17</sup></td><td>0.05</td></tr><tr><td>2.1-2.4</td><td>1.03 × 10<sup>-39</sup></td><td>5.15 × 10<sup>-24</sup></td><td></td><td></td><td></td><td></td></tr><tr><td>2.4-2.7</td><td>1.45 × 10<sup>-36</sup></td><td>2.02 × 10<sup>-24</sup></td><td>6.61 × 10<sup>-30</sup></td><td>4.59 × 10<sup>-29</sup></td><td>1.25 × 10<sup>-11</sup></td><td>0.08</td></tr><tr><td>2.7-3.0</td><td>3.21 × 10<sup>-24</sup></td><td>6.22 × 10<sup>-23</sup></td><td></td><td></td><td></td><td></td></tr><tr><td>2.7-3.4</td><td>3.90 × 10<sup>-16</sup></td><td>1.38 × 10<sup>-10</sup></td><td>0.001</td><td>3.28 × 10<sup>-9</sup></td><td>5.95 × 10<sup>-5</sup></td><td>0.86</td></tr></tbody></table>

<figure><img src="image_10.png" /><figcaption>**Figure 8.** The BH mass estimates provided from Rakshit et al. (2020) for the RL (dark blue) and RQ population matched in redshift and *i*-band magnitude (light blue). Each panel corresponds to the three different methods of calculation, which depend on the availability of a given emission line: the *H*β line (z<0.8), the Mg II line (0.8≤z < 1.9), and the C IV line (z ≥ 1.9). The sample means for the RL and RQ QSOs are presented in each panel as dashed dark blue and solid blue lines, respectively.</figcaption></figure>

The resulting RL/RQ comparison is found to exhibit similar features as for the *z* − *M*<sub>i</sub> sample, where both the continuum and the emission lines are found to differ between the two populations (see Appendix A for the RL/RQ comparison with the *z* − *M*<sub>i</sub> − *M*<sub>BH</sub> sample). We use the null hypothesis test as described in the previous section to determine the significance of this result. Similarly to the *z* − *M*<sub>i</sub> sample, we are able to rule out the null hypothesis with *p* < 0.05 for all redshift bins (see results for the 3D-matched sample in Table 1). Therefore, we conclude that additional information beyond the BH mass and accretion rate is required to explain the difference between the RL and RQ populations. As matching in *M*<sub>BH</sub> is on average of little consequence (except for sample size), we proceed with the rest of our analysis using only the *z* − *M*<sub>i</sub> sample, since it is about three times larger.

## 4.3 Comparison as a function of radio loudness

To investigate whether the differences observed between the RL and RQ populations represent a bimodality or a transition between populations, we further separate our sample into four radio-loudness bins as indicated in Fig. 2. As with our previous approach, we include

radio-undetected sources only in the lowest radio-loudness bin (*R*<sub>1</sub>) as we cannot confidently assign the rest of the sources to the correct *R* bin (i.e when the upper limits fall to the right of the demarcation lines in Fig. 2). This leads to some challenges when matching in *z* and *M*<sub>i</sub> as demonstrated by Fig. 9, which presents the redshift-luminosity plane. We can see that the more RL bins (*R*<sub>2</sub>, *R*<sub>3</sub>, and *R*<sub>4</sub>) in comparison to *R*<sub>1</sub> span a different parameter space. This means that we need to match to the lowest radio-loudness bin sample, unlike before where the matching was performed to the RL population in order to maximize the radio-loudest sample. As this resulted in a lower number of sources for each radio-loudness bin, we decided to reduce the number of redshift bins by doubling the bin width to maintain comparable S/N levels of the stacked spectra as before.

The results of applying our spectral stacking algorithm and the null hypothesis test on this sample division for five redshift bins show that all RL bins exhibit statistically significant differences when compared to *R*<sub>1</sub> (see results for *R*<sub>1</sub>−*R*<sub>2</sub>/*R*<sub>3</sub>/*R*<sub>4</sub> in Table 1 and Appendix A), providing no indication of a smooth transition between populations relative to the null hypothesis test. To address the possibility of a bias arising from our selection process, where only radio-detected sources were included in the more RL bins, we further
<figure><img src="image_11.png" /><figcaption>Figure 9. The $z - M_i$ plane for the different radio-loudness bins. The lowest radio-loudness bin ($R_1$) is shown in light blue, whereas the bins with increasing R are presented in dark blue throughout the three panels.</figcaption></figure>

conducted a comparison between radio-detected ($R_{1D}$) and radio-undetected sources ($R_{1U}$) within the lowest radio-loudness bin ($R_1$). We find that we can accept the null hypothesis that the two samples are drawn from the same parent distribution ($p \ge 0.05$) for all, apart from the lowest redshift bin, where $p = 0.01$ indicating ambiguity. This ambiguous result, however, is not enough to explain p-values as low as $10^{-41}$ as found for the $R_1-R_2$ comparison. Therefore, we conclude that the exclusion of non-detections is not significantly impacting our findings.

We must also consider the possible role of the host galaxy (e.g. Magliocchetti et al. 2020). Given our high S/N, we expect to be able to detect starlight up to redshifts normally inaccessible to SDSS; light from the hosts could therefore be contributing to our results. To test this, we examine the excess region of flux density for each pair per given redshift bin. Our findings indicate that the luminosity of this region is of order $10^{9-10} L_\odot$, suggesting that the host galaxy might indeed play a significant role in generating the observed effects.

Another possible explanation for the absence of a smooth transition between populations in terms of the hypothesis test results is that we are searching for differences only in the average spectrum in each bin. Since Macfarlane et al. (2021) showed that there exists a significant population of jetted QSOs even at the lowest R values, it is possible that the differences we see are a result of an increasing fraction of jetted AGN which become numerically dominant in all but the lowest R bin, and which therefore dominate the median stacked spectra. This hypothesis will become testable in the future as we are increasingly able to morphologically discern the presence of jets in forthcoming wide field subarcsecond 144 MHz imaging from the LOFAR international stations.

## 5 SPECTRAL PROPERTIES OF RL AND RQ QUASARS

In this section, we employ the PYQSOFIT spectral fitting code to obtain spectral properties from the quasar populations. In Section 5.2, we analyse the average emission-line properties of both RL and RQ quasars, as well as those categorized into different radio-loudness bins, while in Section 5.3, we focus on investigating the nature of the [O II] excess found for the RL regime.

### 5.1 The spectral fitting procedure

To investigate the continuum- and emission-line properties of the composite spectra of RL and RQ quasars, we use the PYQSOFIT spectral fitting code in a similar manner to previous works (e.g. Shen et al. 2019, Rakshit et al. 2020, Fawcett et al. 2022), where the continuum is fitted globally (i.e. the regions influenced by the presence of emission lines are masked out), whereas the emission-line complexes are fitted separately and locally.

The continuum of each composite spectrum is globally modelled by using a combination of power law, Balmer continuum, Fe II component, and third-order polynomial as:

$$f_{\text{conti}} = f_{\text{pl}} + f_{\text{BC}} + f_{\text{Fe II}} + f_{\text{poly}} \quad (4)$$

The power-law continuum component has been added to represent the emission from the accretion disc and is defined as:

$$f_{\text{pl}} = a_0(\lambda/\lambda_0)^{a_1}, \quad (5)$$

where $a_0$ and $a_1$ are the normalization and power-law slope, and $\lambda_0$ is a reference wavelength at 3000 Å.

The Balmer continuum component represents the sum of blended, higher-order Balmer lines that give rise to the well-known small blue bump at $\lambda \sim 3000$ Å. In PYQSOFIT, it is modelled by the function given by Grandi (1982) for the case of optically thick clouds as:

$$f_{\text{BC}} = F_{\text{BE}} B_{\lambda}(T_e)(1 - e^{-\tau_{\lambda}}(\lambda/\lambda_{\text{BE}})^3), \quad (6)$$

where $\lambda_{\text{BE}}$ is the position of the Balmer edge, $F_{\text{BE}}$ is the flux at the Balmer edge, $B_{\lambda}(T_e)$ is the Planck function at the electron temperature $T_e$, and $\tau_{\lambda}$ is the optical depth. Here, $F_{\text{BE}}$, $T_e$, and $\tau_{\lambda}$ are left as free parameters as in Fawcett et al. (2022).

The Fe II component is another essential part of modelling the continuum. At UV and optical wavelengths, AGN spectra contain numerous iron emission lines that blend together to form a pseudo-continuum. If not properly subtracted, this pseudo-continuum could contaminate the continuum and emission-line measurements. Here, PYQSOFIT models it as:

$$f_{\text{Fe II}} = c_0 F_{\text{Fe II}}(\lambda, c_1, c_2), \quad (7)$$

where $c_0$ is the normalization constant, $c_1$ is a constant describing the Gaussian broadening, and $c_2$ represents the wavelength shift
<table><caption>**Table 2.** The emission-line fitting information. The columns are as follows: the name of the line complexes, the wavelength range used for the fit, the emission lines present in each line complex, and the number of Gaussian components used for each line.</caption><thead><tr><th>Emission-line<br>complex<br>name</th><th>Wavelength<br>range<br>(Å)</th><th>Emission<br>line<br>name</th><th>Number of<br>Gaussian<br>components</th></tr></thead><tbody><tr><td rowspan="3">Hβ</td><td rowspan="3">4640–5100</td><td>Hβ broad</td><td>2</td></tr><tr><td>Hβ narrow</td><td>1</td></tr><tr><td>[O III] (4959 Å)</td><td>1</td></tr><tr><td rowspan="3">Hγ</td><td rowspan="3">4200–4440</td><td>[O III] (5007 Å)</td><td>1</td></tr><tr><td>Hγ broad</td><td>1</td></tr><tr><td>Hγ narrow</td><td>1</td></tr><tr><td>[O II]</td><td>3650–3800</td><td>O III (4363 Å)</td><td>1</td></tr><tr><td>[Ne V]</td><td>3380–3480</td><td>[O II] (3728 Å)</td><td>1</td></tr><tr><td rowspan="2">Mg II]</td><td rowspan="2">2700–2900</td><td>Mg II] broad</td><td>2</td></tr><tr><td>Mg II] narrow</td><td>1</td></tr><tr><td rowspan="2">C III]</td><td rowspan="2">1810–1970</td><td>C III] broad</td><td>2</td></tr><tr><td>C III] narrow</td><td>1</td></tr><tr><td rowspan="7">C IV</td><td rowspan="7">1500–1700</td><td>C IV broad</td><td>2</td></tr><tr><td>C IV narrow</td><td>1</td></tr><tr><td>He II (1640 Å) broad</td><td>1</td></tr><tr><td>He II (1640 Å) narrow</td><td>1</td></tr><tr><td>O III] (1663 Å) broad</td><td>1</td></tr><tr><td>O III] (1663 Å) narrow</td><td>1</td></tr><tr><td>Si IV + O IV]</td><td>1290–1450</td><td>Si IV + O IV]</td><td>1</td></tr><tr><td rowspan="3">Lyα</td><td rowspan="3">1150–1290</td><td>Lyα broad</td><td>2</td></tr><tr><td>Lyα narrow</td><td>1</td></tr><tr><td>NV (1240 Å)</td><td>1</td></tr><tr><td>O VI</td><td>980–1100</td><td>O VI</td><td>1</td></tr></tbody></table>

applied to the Fe II templates to match the data. For the UV part of the spectrum, we use the modified UV Fe II template by Shen et al. (2019) that combines the Vestergaard & Wilkes (2001) template for rest-frame wavelengths 1000–2200 Å, the Salviander et al. (2007) template for 2200–3090 Å, and the Tsuzuki et al. (2006) template for 3090–3500 Å. The optical Fe II template, on the other hand, is taken from Véron-Cetty, Joly & Véron (2004), covering the rest-frame wavelengths of 3535–7534 Å.

Finally, the third-order polynomial is used to account for the bending of the continuum, which is likely caused by the intrinsic dust reddening of the population, and is defined as:

$$f_{\text{poly}} = \sum_{i=0}^{i=3} b_i (\lambda - \lambda_0)^i, \qquad (8)$$

where $b_i$ are the model free parameters.

After subtracting the best-fitting model of the continuum, we fit the emission-line spectrum in log-space. The fitting is performed individually for each emission-line complex, where all the emission lines contained within a single complex are fit simultaneously. The narrow components (FWHM < 1200 km s⁻¹) are modelled by a single Gaussian, where the velocity offset and line width are tied for each line complex. The broad emission-line profiles, however, can be quite complex (e.g. double-peaked, with a flat top, or asymmetric) and thus are not well represented by a single Gaussian. In such cases, we use multiple Gaussian components depending on the line complexity. A full list of the line complexes containing the individual emission lines and number of Gaussian components used in the fit is provided in Table 2.

From the best-fitting models, we can obtain emission-line fluxes for each line in question. The uncertainties of these measurements are calculated using the Monte Carlo approach embedded within PYQSOFIT, where we choose to use 100 iterations.

## 5.2 Average emission-line properties

As mentioned in Section 2.1, PYQSOFIT contains an additional feature that separates a given quasar spectrum into host galaxy and pure quasar contribution. However, for this separation to be well represented, the continuum needs to be accurately fitted. Unfortunately, the host and QSO continuum models implemented in PYQSOFIT prove inadequate for our objectives, as demonstrated by poor fitting statistics (see example fit in Appendix A). They appear to be strongly dependent on the details of the wavelength range used, indicative of a local minimum problem. Moreover, these high $\chi^2$ values have additional implications as they prevent us from applying dust-extinction laws to investigate the underlying cause of the continuum differences observed in Section 4.1. We are, however, able to obtain acceptable fits for the emission lines. To improve their goodness of fit, we use an additional localized power-law continuum component for each line complex after subtracting the globally fitted continuum, following the method of Fawcett et al. (2022). In addition, we perform a flux calibration before the fit to obtain non-arbitrary units. This is done by shifting each composite spectrum to the observed frame by using its central redshift value ($z_{\text{centre}}$) and the SDSS i-band filter curve and the median apparent i-band magnitude of the individual spectra to obtain a flux density in units of erg s⁻¹ cm⁻² Å⁻¹.

Fig. 10 presents the ratio of the emission-line flux of the RL to RQ QSOs (RL–RQ), along with the more RL to radio-quietest bin ($R_2/R_3/R_4 - R_1$) for the most prominent emission lines. Across all considered pairs, we find that the majority of emission lines exhibit flux ratios close to unity. However, there is a strong excess for RL objects in the [Ne V] 3426 Å, [O II] 3728 Å, and the [O III] 4960/5008 Å (hereafter [Ne V], [O II], and [O III]) emission lines. These lines are of particular interest as [Ne V] is used to trace extended emission-line regions (EELR: see discussion in next section), whereas the [O III] emission line is generally associated with the AGN bolometric luminosity (Heckman et al. 2004; Best & Heckman 2012; Kalfountzou et al. 2012), as it is mostly ionized by the continuous radiation coming from the accretion disc. The [O II] emission line, on the other hand, is a well-known SF tracer in quasar host galaxies (e.g. Ho 2005; Kalfountzou et al. 2012; Matsuoka et al. 2015) as in contrast with the [O III], it is only weakly excited in the narrow-line region (NLR). However, as result of the low number of data points for the [O III] emission due to the available redshift range, we focus only on [Ne V] and [O II] in the following section.

## 5.3 Star formation rates

The [O II] forbidden line is a prominent feature that can be easily detected in moderate resolution spectra, making it a good SF tracer in the absence of $H\alpha$. When an AGN is present, however, this emission line can be severely contaminated by excitation from EELR which can span out to several kpc (e.g. Unger et al. 1987; Villar-Martín et al. 2011; Husemann et al. 2014; Maddox 2018). Therefore, to determine whether the elevated levels of [O II] found in Section 5.2 are due to SF, we employ the AGN correction technique described in Maddox (2018). Briefly, this method involves removing sources from the stacks with an [Ne V] detection and making use of the [O II]-to-[Ne V] emission-line flux ratio, which in a typical AGN-dominated
<figure><img src="image_12.png" /><figcaption>**Figure 10.** The emission-line flux ratio of a given pair, as indicated in the legend, as a function of redshift. Each panel represent an emission line under comparison as denoted in the upper right corner, which must be significantly detected (>3σ) for inclusion. The dashed line in each panel represents the line of equality, facilitating visual comparison of the emission-line flux ratios between the different pairs.</figcaption></figure>

galaxy is found to be of order unity as a result of excitation from NLR. The [NeV] is chosen because of its wavelength proximity to [OII] and its high ionization potential, making the effects of dust attenuation and star-forming regions negligible.

As the spectral measurements provided by Rakshit et al. (2020) did not contain any information regarding [NeV], we used PYQSOFIT to fit a single Gaussian in a spectral window of 150 Å centred on the [NeV] central wavelength. The emission line is considered as a non-detection if the measured emission-line flux is detected at an S/N<3. Although in previous sections we investigated how the optical properties varied in bins of radio loudness, the need to remove the [NeV] detections from the stacks (which involves about 20 per cent) reduces the sample size to a point that it is no longer possible to search for evolution in both redshift and radio loudness. Therefore, we continue the rest of the analysis with the RL and RQ populations. However, we also create a new sample matched in z and $M_*$, with the division set at $R = 0.2$ (i.e. the separation between $R_1$ and the more RL bins), to check whether using a slightly different classification would yield different results. From the newly computed composite spectra, we obtained emission-line measurements for [OII] and [NeV] in the same manner as before. An [NeV] detection is still present in the stacks as a result of the high S/N which is sufficient to detect the contribution from the NLR.

Using these measurements and the value of [OII]/[NeV] = 1.05 obtained from the SDSS composite spectrum from Vanden Berk et al. (2001), we are able to derive a measure of the [OII] emission-line flux unaffected by AGN.

Fig. 11 shows the [OII] derived, doubly corrected star formation rate (SFR) against redshift obtained for the RL and RQ populations using the [OII]–SFR conversion from Kewley, Geller & Jansen (2004). We can see that the SFR of the RQ population is close to zero, indicating that on average SF processes in the host galaxy of RQ QSOs are not sufficient to produce a strong [OII] emission line. In contrast, the excess of [OII] persists for the RL population even when correcting for contamination from EELR and NLR. The alternative classification ($R > 0.2$) also exhibits the same trend, further indicating that the exact value of $R$ is not important. Our results are consistent with those of Maddox (2018) who used SDSS DR7 with Faint Images of the Radio Sky at Twenty-cm (Becker, White & Helfand 1995), but we have shown that their findings persist with higher statistical significance and to a further ≈1.4 Gyr of cosmic history, closer to the peak epoch of cosmic SFR density (Madau & Dickinson 2014). This is because we benefit not only from the increased sensitivity of LoTSS and the reduced impact of line-of-sight effects (relative to the 1.4 GHz data), but also the extended wavelength coverage of the BOSS spectrograph (Smee et al. 2013).
<figure><img src="image_13.png" /><figcaption>**Figure 11.** The SFR obtained from the doubly corrected [O II] emission-line flux for the RL (dark blue) and RQ population (light blue) compared against the SFR obtained from Maddox (2018). The black dashed line indicates an SFR of zero.</figcaption></figure>

# 6 DISCUSSION

## 6.1 Could our results be caused by starburst galaxies contaminating our RL sample?

As discussed in the Introduction, Macfarlane et al. (2021) have shown that the radio flux density of quasars can be modelled by a combination of SF and radio jets. This naturally leads to the question of whether the [O II] excess found in Section 5.3 could be caused by having RQ quasars hosted by starburst galaxies in our RL sample. The presence of intense SF in such hosts can lead to a larger degree of reddening (Gordon, Calzetti & Witt 1997), which would further explain the redder continuum observed in Section 4.1. In addition, this scenario does not require any differences between the BH mass and accretion rate of the QSOs themselves (i.e. $R_{edd}$), which is in agreement with our results (see Section 4.2).

To explore this possibility, we used the main-sequence relation from Schreiber et al. (2015) and an offset of 0.6 dex based on the starburst criterion from Rodighiero et al. (2011) to obtain an SFR value for a typical starburst galaxy with a stellar mass of $10^{11} M_{\odot}$ at a redshift equal to our central bin values from Fig. 11. This resulted in SFR values ranging from 127 to 488 $M_{\odot} \text{ yr}^{-1}$ for the redshift range of $0.6 < z < 1.9$. These estimates are then converted to a radio luminosity at 150 MHz by using the mass independent SFR-$L_{150MHz}$ relation from Smith et al. (2021). Comparing these results to the RL sample, we find that up to 5 per cent in any given redshift bin, could have a radio luminosity consistent with contribution from RQ starbursts. However, taking into account the large difference between our doubly corrected [O II] SFR and the values obtained here for the starburst galaxies, it is unlikely that the radio emission is due to such extreme star-forming processes. Furthermore, considering the fact that we are investigating the median stacks of the populations, we believe that if starburst galaxies are present in our samples they will not noticeably impact our results.

## 6.2 Could shock excitation be playing a role?

In Section 5.3, we have used the [NeV] emission line to remove contamination from EELR and correct for contribution from the AGN NLR. However, shocks in the interstellar medium generated

<figure><img src="image_14.png" /><figcaption>**Figure 12.** The [O II] emission-line luminosity as a function of 178 MHz total radio luminosity for RL sources without significant [NeV] detections. The solid and dashed black lines are generated from the regression analysis results for the 3CRR radio sources in table 6 of Hardcastle et al. (2009), where the solid line indicates the regression line, whereas the dashed lines show $\pm 3 \times$ the obtained scatter.</figcaption></figure>

by radio jets or AGN-driven winds could also contribute to the [O II] emission. To assess their potential influence, Maddox (2018) used the MAPPINGS III shock and photoionization modelling code from Allen et al. (2008). The author found that there is a variety of shock conditions capable of exciting [O II], but only high velocity shocks ($>600 \text{ km s}^{-1}$) are able to produce [NeV]. This suggests that the doubly corrected [O II] could still be contaminated from moderate-velocity shocks. A prevalence of such shocks for the RL population may be able to explain the apparent enhancement of SFR and their redder appearance.

To investigate this scenario, we explore the connection between [O II] emission and 178 MHz luminosity, where we use the low-frequency radio emission as a proxy for jet power. Here, we use a radio spectral index of $\alpha = -0.7$ to convert from 144 to 178 MHz, in order to compare our results with Hardcastle, Evans & Croston (2009), who have found a positive correlation for a sample of 3CRR radio sources, which is thought to be of nuclear origin. The results for the RL sample from Section 5.3, which was used to calculate the [O II] SFRs, are presented in Fig. 12. We can see that the [O II] appears to be independent of the radio emission at any given redshift. However, there is a subset of sources that lie within the region defined by the 3CRR sample. This indicates the presence of some AGN-related influence, which can be as high as 36 per cent for the lowest redshift range ($0.6 < z < 0.8$), and up to 13 per cent for the highest range ($1.7 < z < 1.9$). But, as previously discussed our analysis relies on SFR derived from the median composite spectra of the populations. Furthermore, the correction for the AGN NLR is not factored into the individual RL sources, which could potentially explain the correlation observed in the 3CRR sources. Therefore, if shocks are present in the RL population, we believe that they alone cannot account for the observed differences between RL and RQ quasars.

## 6.3 Possible explanations

To explain the differences between the RL and RQ populations, we propose two distinct models: one centred on BH spin dynamics (Blandford & Znajek 1977; Wilson & Colbert 1995; Sikora et al. 2007) and another on the evolutionary scenario proposed for red and
blue quasars (e.g. Sanders, Soifer & Neugebauer 1990; Hopkins et al. 2008; Alexander & Hickox 2012; Klindt et al. 2019; Fawcett et al. 2020, 2022).

The first model involves a rapidly spinning BH, coupled with a rich gas supply required for accretion and jet production (e.g. Hardcastle, Evans & Croston 2007; Gürkan et al. 2015). Given that the gas needed for accretion, which leads to the radio loudness of the source, is also essential for fuelling SF, these two processes are inherently interconnected. Therefore, the question becomes how to get a rapidly spinning BH and a rich gas supply to make an RL quasar. We suggest two plausible scenarios: first, the traditional major merger event which results in a rapid inflow of gas and dust. This material not only transports angular momentum toward the central region, leading to the rapid spin-up of the BH, but also serves as a trigger for SF. Secondly, we can consider a scenario in which RL quasars are hosted by massive, gas-rich spiral galaxies. In this setting, a substantial reservoir of cold gas is available, supporting continuous SF and accretion. Furthermore, due to the ordered rotation characteristic of spiral galaxies, it may give rise to an efficient transfer of angular momentum onto the supermassive BH, eliminating the necessity for merger events.

The alternative model assumes once more a gas-rich merger triggering an AGN, but here we focus on different stages of evolution, following, for example, Hopkins et al. (2008). The first one is a relatively short-lived phase where the QSO is heavily obscured by high-column gas density and dust. Subsequently, the AGN generates powerful winds and/or outflows, which disperse the obscuring material. In our study, RL quasars could represent the former evolutionary stage, where the redder continuum and enhanced SFR would be explained by the obscuring and dense material, whereas RQ quasars may fall into the latter (unobscured) category. Here, we do not require any difference in either the accretion rate or in the BH spins between the two classes. This scheme is related to the one presented for red and blue QSOs discussed in previous studies (e.g. Klindt et al. 2019; Fawcett et al. 2020, 2022), where dividing the QSO population according to their optical colours (rather than radio loudness) gives a red class with a significant radio flux excess relative to the blue class. However, both the red and blue QSO classes contain RL and RQ sources, and the average QSO in both classes is consistent with being RQ (i.e. *R* < 1). It is therefore clear that this association alone cannot explain our results (certainly we cannot equate the RL QSOs with the 'red QSO' class, etc).

To make further progress on what controls the radio-loudness of QSOs, we need more information. This could come from the subarcsecond 144 MHz imaging that is now becoming possible with LOFAR (and the morphological information that it can provide; e.g. Morabito et al. 2022), along with larger statistical samples from new and forthcoming facilities such as WEAVE (Dalton et al. 2012; Jin et al. 2023), the Dark Energy Spectroscopic Instrument (Aghamousa et al. 2016a, b), and the Multi-object Optical and Near-IR spectrograph (Cirasuolo et al. 2014). In addition, improved BH spin estimates from X-ray observatories such as the European Space Agency's Athena X-ray observatory (Barcons et al. 2017) and NASA's Nuclear Spectroscopic Telescope Array (Harrison et al. 2010) will be crucial in deepening our understanding of this subject.

# 7 SUMMARY AND FUTURE PROSPECTS

In this work, we have used the the LoTSS DR2 and the SDSS DR14 to create the largest, uniformly selected, spectroscopically confirmed sample of RL and RQ quasars. To study their spectroscopic properties, we have developed a new spectral stacking code which

accounts for a range of biases including those that arise as a result of the redshifting process. Such a tool not only allows us to robustly compare quasar and galaxy populations, but also enables us to statistically recover the continuum properties of faint sources. This will become particularly important for radio-selected spectroscopic surveys such as the WEAVE-LOFAR survey (Smith et al. 2016), which will generate more than one million optical spectra, with continuum detections absent in a significant fraction. Using this algorithm to investigate the average properties of QSOs as a function of their radio loudness, we have found the following results:

(i) The high S/N composite spectra representative of the RL and RQ populations differ. RL QSOs are found to have on average a redder continuum with an [OII] emission-line excess across the redshift range of $0.6 < z < 3.5$. Such differences highlight the importance of creating high-resolution stacks of both populations to improve the redshift classification of future spectroscopic surveys.

(ii) The RL and RQ populations are found to span similar BH mass ranges with no systematic trend showing that the mean $M_{BH}$ is higher for RL QSOs. Furthermore, using a sample matched in $z, M_i$, and $M_{BH}$ to make a comparison between the average spectra of RL and RQ QSOs is found to give similar statistically consistent results as for a sample matched in $z$ and $M_i$. This suggests that neither the BH mass, nor the accretion rate are defining factors in a QSO's radio loudness.

(iii) The observed differences between the RL and RQ populations are not gradual. Comparing composite spectra as a function of radio loudness shows that all more RL bins ($R_2, R_3, R_4$) differ from the radio quietest ($R_1$), with features consistent with the classical RL and RQ divisions.

(iv) These changes cannot be explained by the addition of radio-undetected sources, as we have found that the radio-detected and radio-undetected quasars in $R_1$ are consistent with being drawn from the same parent distribution.

(v) We have shown that RL quasars have on average higher SFRs than their RQ counterparts at any given redshift $0.5 < z < 1.9$, extending this result to significantly earlier cosmic epochs than previously known. The elevated levels of [OII] emission that we use to infer the SFRs have been corrected for possible influence of AGN contamination (following the procedure of Maddox 2018) and we have shown that our results cannot be explained by contamination from starburst galaxies.

Our results show that there is no clear-cut division in radio loudness between RL and RQ quasars, and that the differences observed between them is not related to BH mass or the accretion rate. As a result, we propose two distinct models: one requiring that RL quasars have rapidly spinning BHs in conjunction with abundant gas reservoirs, or are representatives of an earlier obscured phase of QSO evolution.

With the advent of future spectroscopic surveys such as WEAVE-LOFAR, the number of RL QSOs will significantly increase. This will allow to us to investigate their spectral properties in greater detail and with higher significance. With higher S/N, however, we will need more sophisticated theoretical models to fit the composite spectra in order to disentangle their dust and host properties.

# ACKNOWLEDGEMENTS

We thank the anonymous referee for the helpful comments that have improved the manuscript. MIA acknowledges support from the UK Science and Technology Facilities Council (STFC) studentship under the grant ST/V506709/1. DJBS, MJH, and ABD
acknowledge support from the STFC under the grant ST/V000624/1. KJD acknowledges support from the STFC through an Ernest Rutherford Fellowship (grant no. ST/W003120/1). LKM is grateful for support from the Medical Research Council (MR/T042842/1). SD acknowledges support from the STFC via studentship grant no. ST/W507490/1. SS acknowledges support from the STFC via studentship grant no. ST/X508408/1. LOFAR is the Low Frequency Array designed and constructed by ASTRON. It has observing, data processing, and data storage facilities in several countries, which are owned by various parties (each with their own funding sources), and that are collectively operated by the ILT foundation under a joint scientific policy. The ILT resources have benefitted from the following recent major funding sources: CNRS-INSU, Observatoire de Paris and Université d'Orléans, France; BMBF, MIWFNRW, MPG, Germany; Science Foundation Ireland (SFI), Department of Business, Enterprise and Innovation (DBEI), Ireland; NWO, The Netherlands; The Science and Technology Facilities Council, UK; Ministry of Science and Higher Education, Poland; The Istituto Nazionale di Astrofisica (INAF), Italy. This research made use of the Dutch national e-infrastructure with support of the SURF Cooperative (e-infra 180169) and the LOFAR e-infra group. The Jülich LOFAR Long Term Archive and the German LOFAR network are both coordinated and operated by the Jülich Supercomputing Centre (JSC), and computing resources on the supercomputer JUWELS at JSC were provided by the Gauss Centre for Supercomputing e.V. (grant CHTB00) through the John von Neumann Institute for Computing (NIC). This research made use of the University of Hertfordshire high performance computing facility and the LOFAR-UK computing facility located at the University of Hertfordshire and supported by STFC [ST/V002414/1], and of the Italian LOFAR IT computing infrastructure supported and operated by INAF, and by the Physics Department of Turin University (under an agreement with Consorzio Interuniversitario per la Fisica Spaziale) at the C3S Supercomputing Centre, Italy.

## DATA AVAILABILITY

The data used in this work are publicly available and can be found at the LOFAR Surveys (www.lofar-surveys.org) and the SDSS website (www.sdss.org.)

## REFERENCES

* DESI Collaboration, 2016a, preprint (arXiv:1611.00036)

* DESI Collaboration, 2016b, preprint (arXiv:1611.00037)

* Alexander D. M., Hickox R. C., 2012, New Astron. Rev., 56, 93

* Allen M. G., Groves B. A., Dopita M. A., Sutherland R. S., Kewley L. J., 2008, ApJS, 178, 20

* Baloković M., Smolčić V., Ivezić Ž., Zamorani G., Schinnerer E., Kelly B. C., 2012, ApJ, 759, 30

* Barcons X. et al., 2017, Astron. Nachr., 338, 153

* Becker R. H., White R. L., Helfand D. J., 1995, ApJ, 450, 559

* Best P., Heckman T., 2012, MNRAS, 421, 1569

* Best P., Kauffmann G., Heckman T., Brinchmann J., Charlot S., Ivezić Ž., White S., 2005, MNRAS, 362, 25

* Blandford R. D., Znajek R. L., 1977, MNRAS, 179, 433

* Bovy J. et al., 2011, ApJ, 729, 141

* Brinchmann J., Charlot S., White S., Tremonti C., Kauffmann G., Heckman T., Brinkmann J., 2004, MNRAS, 351, 1151

* Bruzual G., Charlot S., 2003, MNRAS, 344, 1000

* Calabrò A. et al., 2021, A&A, 646, A39

* Carnall A., 2017, preprint (arXiv:1705.05165)

* Chiaberge M., Marconi A., 2011, MNRAS, 416, 917

* Chonis T. S. et al., 2016, in Evans C. J., Simard L., Takami H., eds, Proc. SPIE Conf. Ser. Vol. 9908, Ground-based and Airborne Instrumentation for Astronomy VI. SPIE, Bellingham, p. 1345

* Cirasuolo M., Magliocchetti M., Celotti A., Danese L., 2003a, MNRAS, 341, 993

* Cirasuolo M., Celotti A., Magliocchetti M., Danese L., 2003b, MNRAS, 346, 447

* Cirasuolo M. et al., 2014, in Ramsay S. K., McLean I. S., Takami H., eds, Proc. SPIE Conf. Ser. Vol. 9147, Ground-based and Airborne Instrumentation for Astronomy V. SPIE, Bellingham, p. 202

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

* Condon J., Kellermann K., Kimball A. E., Ivezić Ž., Perley R., 2013, ApJ, 768, 37

* Dalton G. et al., 2012, in McLean I. S., Ramsay S. K., Takami H., eds, Proc. SPIE Conf. Ser. Vol. 8446, Ground-based and Airborne Instrumentation for Astronomy IV. SPIE, Bellingham, p. 220

* De Ruiter H., Willis A., Arp H., 1977, A&AS, 28, 211

* Dunlop J., McLure R., Kukula M., Baum S., O'Dea C., Hughes D., 2003, MNRAS, 340, 1095

* Fawcett V., Alexander D., Rosario D., Klindt L., Fotopoulou S., Lusso E., Morabito L., Calistro Rivera G., 2020, MNRAS, 494, 4802

* Fawcett V., Alexander D., Rosario D., Klindt L., Lusso E., Morabito L., Calistro Rivera G., 2022, MNRAS, 513, 1254

* Fitzpatrick E. L., 1999, PASP, 111, 63

* Fynbo J., Krogager J.-K., Venemans B., Noterdaeme P., Vestergaard M., Møller P., Ledoux C., Geier S., 2012, ApJS, 204, 6

* Glikman E., Helfand D. J., White R. L., 2006, ApJ, 640, 579

* Cloudemans A. et al., 2022, A&A, 668, A27

* Gordon K. D., Calzetti D., Witt A. N., 1997, ApJ, 487, 625

* Grandi S., 1982, ApJ, 255, 25

* Gu M., Cao X., Jiang D., 2001, MNRAS, 327, 1111

* Guo H., Shen Y., Wang S., 2018, Astrophysics Source Code Library, record ascl:1809

* Gürkan G. et al., 2015, MNRAS, 452, 3776

* Gürkan G. et al., 2019, A&A, 622, A11

* Hardcastle M. J., Evans D., Croston J., 2007, MNRAS, 376, 1849

* Hardcastle M., Evans D., Croston J., 2009, MNRAS, 396, 1929

* Hardcastle M. et al., 2016, MNRAS, 462, 1910

* Hardcastle M. et al., 2023, A&A, 678, A151

* Harris D. W. et al., 2016, AJ, 151, 155

* Harrison F. A. et al., 2010, in Arnaud M., Murray S. S., Takahashi T., eds, Proc. SPIE Conf. Ser. Vol. 7732, Space Telescopes and Instrumentation 2010: Ultraviolet to Gamma Ray. SPIE, Bellingham, p. 189

* Heckman T. M., Kauffmann G., Brinchmann J., Charlot S., Tremonti C., White S. D., 2004, ApJ, 613, 109

* Hill G. J. et al., 2008, in Kodama T., Yamada T., Aoki K., eds, ASP Conf. Ser. Vol. 399, Panoramic Views of Galaxy Formation and Evolution. Astron. Soc. Pac., San Francisco, p. 115

* Ho L. C., 2002, ApJ, 564, 120

* Ho L. C., 2005, ApJ, 629, 680

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

* Husemann B., Jahnke K., Sánchez S., Wisotzki L., Nugroho D., Kupko D., Schramm M., 2014, MNRAS, 443, 755

* Ivezić Ž. et al., 2002, AJ, 124, 2364

* Jin S. et al., 2023, MNRAS, 1557

* Kalfountzou E., Jarvis M., Bonfield D., Hardcastle M., 2012, MNRAS, 427, 2401

* Kashikawa N. et al., 2002, PASJ, 54, 819

* Kauffmann G. et al., 2003, MNRAS, 341, 33

* Kellermann K., Sramek R., Schmidt M., Shaffer D., Green R., 1989, AJ, 98, 1195

* Kellermann K., Condon J., Kimball A., Perley R., Ivezić Ž., 2016, ApJ, 831, 168

* Kewley L. J., Geller M. J., Jansen R. A., 2004, AJ, 127, 2002

* Kewley L. J., Groves B., Kauffmann G., Heckman T., 2006, MNRAS, 372, 961

* Kimball A. E., Kellermann K. I., Condon J. J., Ivezić Ž., Perley R. A., 2011, ApJ, 739, L29
Klindt L., Alexander D. M., Rosario D. J., Lusso E., Fotopoulou S., 2019, MNRAS, 488, 3109

Kondapally R. et al., 2021, A&A, 648, A3

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

Liu Y., Jiang D. R., Gu M. F., 2006, ApJ, 637, 669

Macfarlane C. et al., 2021, MNRAS, 506, 5888

Madau P., Dickinson M., 2014, ARA&A, 52, 415

Maddox N., 2018, MNRAS, 480, 5203

Magliocchetti M. et al., 2020, MNRAS, 493, 3838

Matsuoka Y. et al., 2015, ApJ, 811, 91

McLure R. J., Jarvis M. J., 2004, MNRAS, 353, L45

Metcalf R., Magliocchetti M., 2006, MNRAS, 365, 101

Morabito L. K. et al., 2022, MNRAS, 515, 5758

Myers A. D. et al., 2015, ApJS, 221, 27

Oke J. et al., 1995, PASP, 107, 375

Pâris I. et al., 2018, A&A, 613, A51

Radcliffe J. F., Barthel P., Garrett M., Beswick R., Thomson A., Muxlow T., 2021, A&A, 649, L9

Rakshit S., Stalin C., Kotilainen J., 2020, ApJS, 249, 17

Richards G. T. et al., 2002, AJ, 123, 2945

Richards G. T. et al., 2006, ApJS, 166, 470

Richter G., 1975, Astron. Nachr., 296, 65

Rigby J. et al., 2018, ApJ, 853, 87

Rodighiero G. et al., 2011, ApJ, 739, L40

Ross N. P. et al., 2012, ApJS, 199, 3

Rowlands K. et al., 2012, MNRAS, 419, 2545

Salviander S., Shields G. A., Gebhardt K., Bonning E., 2007, ApJ, 662, 131

Sanders D., Soifer B., Neugebauer G., 1990, in Dynamics and Interactions of Galaxies: Proceedings of the International Conference, Heidelberg, 29 May–2 June 1989. p. 459

Schlegel D. J., Finkbeiner D. P., Davis M., 1998, ApJ, 500, 525

Schneider D. P. et al., 2010, AJ, 139, 2360

Schreiber C. et al., 2015, A&A, 575, A74

Selsing J., Fynbo J. P., Christensen L., Krogager J.-K., 2016, A&A, 585, A87

Shankar F., Sivakoff G. R., Vestergaard M., Dai X., 2010, MNRAS, 401, 1869

Shen Y. et al., 2019, ApJS, 241, 34

Shen X., Hopkins P. F., Faucher-Giguère C.-A., Alexander D., Richards G. T., Ross N. P., Hickox R., 2020, MNRAS, 495, 3252

Shimwell T. et al., 2017, A&A, 598, A104

Shimwell T. et al., 2022, A&A, 659, A1

Sikora M., Stawarz Ł., Lasota J.-P., 2007, ApJ, 658, 815

Smee S. A. et al., 2013, AJ, 146, 32

Smith D. et al., 2016, in Reylé C., Richard J., Cambrésy L., Deleuil M., Pécontal E., Tresse L., Vauglin I., eds, SF2A-2016: Proceedings of the Annual Meeting of the French Society of Astronomy and Astrophysics, p. 271

Smith D. et al., 2021, A&A, 648, A6

Sutherland W., Saunders W., 1992, MNRAS, 259, 413

Symeonidis M., Giblin B., Page M., Pearson C., Bendo G., Seymour N., Oliver S., 2016, MNRAS, 459, 257

Tchekhovskoy A., Narayan R., McKinney J. C., 2010, ApJ, 711, 50

Tremonti C. A. et al., 2004, ApJ, 613, 898

Tsuzuki Y., Kawara K., Yoshii Y., Oyabu S., Tanabé T., Matsuoka Y., 2006, ApJ, 650, 57

Unger S., Pedlar A., Axon D. J., Whittle M., Meurs E., Ward M., 1987, MNRAS, 228, 671

Urry C. M., Padovani P., 1995, PASP, 107, 803

Vanden Berk D. E. et al., 2001, AJ, 122, 549

Véron-Cetty M.-P., Joly M., Véron P., 2004, A&A, 417, 515

Vestergaard M., Osmer P. S., 2009, ApJ, 699, 800

Vestergaard M., Peterson B. M., 2006, ApJ, 641, 689

Vestergaard M., Wilkes B. J., 2001, ApJS, 134, 1

Villar-Martín M., Tadhunter C., Humphrey A., Encina R. F., Delgado R. G., Torres M. P., Martínez-Sansigre A., 2011, MNRAS, 416, 262

White R. L., Helfand D. J., Becker R. H., Glikman E., de Vries W., 2007, ApJ, 654, 99

White S. V., Jarvis M. J., Häußler B., Maddox N., 2015, MNRAS, 448, 2665

White S. V., Jarvis M. J., Kalfountzou E., Hardcastle M. J., Verma A., Cao Orjales J. M., Stevens J., 2017, MNRAS, 468, 217

Williams W. et al., 2019, A&A, 622, A2

Wilson A. S., Colbert E. J., 1995, ApJ, 438, 62

Yip C.-W. et al., 2004a, AJ, 128, 585

Yip C. et al., 2004b, AJ, 128, 2603

Yue M. et al., 2018, ApJ, 863, 21

Zakamska N. L. et al., 2016, MNRAS, 455, 4191

Zhu G. B. et al., 2015, ApJ, 815, 48
APPENDIX A: SUPPLEMENTARY MATERIAL

<figure><img src="image_15.png" /><figcaption>**Figure A1.** Results of applying our spectral stacking algorithm to sources according to their location on the BPT–N II diagram. Left panel: the BPT–N II diagram with dividing lines indicating regions populated by SFGs (blue), Seyferts (orange), and LINERs (green) as defined relative to the Kauffmann et al. (2003) and Kewley et al. (2006) dividing lines (which are solid and dashed, respectively). Right panel: composite spectra for each class colour-coded to match the regions of the BPT–N II diagram. Emission lines of interest are labelled, and the legend indicates the number of spectra included in each stack.</figcaption></figure>
<figure><img src="image_16.png" /></figure>

**Figure A2.** A Monte Carlo simulation representing the comparison between RL and RQ quasars under the null hypothesis that the RL population is drawn at random from the RQ population. The upper panel of each bin presents the composite spectra of the RL test sample (dark blue) and the RQ test sample (light blue), whereas the lower panel indicates the residual in units of propagated uncertainty (grey). The *p*-values for each of the null hypothesis tests are presented in the upper left corner, while the number of sources and median S/N are indicated in the top right corner of each panel.
<figure><img src="image_17.png" /></figure>

**Figure A3.** A high S/N comparison between the RL and RQ population matched in z, *M*<sub>i</sub>, and *M*<sub>BH</sub>. As in Fig. 7, the upper panel of each bin presents the composite spectra of RL QSOs (dark blue) and RQ QSOs (light blue), whereas the lower panel indicates the residual in units of propagated uncertainty (grey). As for Fig. A2, *p*-values for each of the null hypothesis tests are presented in the upper left corner, while the number of sources and median S/N are indicated in the top right corner of each panel.
<figure><img src="image_18.png" /><figcaption>**Figure A4.** A high S/N comparison between the radio quietest ($R_1$) and more RL ($R_2$, $R_3$ and $R_4$) parts of the QSO population. The upper panel of each bin presents the composite spectra of $R_1$ (black), $R_2$ (light blue), $R_3$ (dark blue), and $R_4$ (green), whereas the lower panel indicates the residual in units of propagated uncertainty for each comparison. As before, the *p*-values for each of the null hypothesis tests are presented in the upper left corner, while the number of sources and median S/N are indicated in the top right corner of each panel.</figcaption></figure>
<figure><img src="image_19.png" /></figure>

**Figure A5.** An example of the spectral fitting procedure with PYQSOFIT. The top panel presents the composite spectra of the RQ population in the redshift range of 0.6 < z < 0.8 (black) overlaid with the best-fitting model (blue). The grey-shaded region show the wavelength windows used for the continuum fit. The lower panels show the best-fitting model of individual line complexes (blue), along with the decomposition into broad (red) and narrow (green) line components.

This paper has been typeset from a TEX/LATEX file prepared by the author.