# Exploring the radio-loudness of SDSS quasars with spectral stacking

M. I. Arnaudova,1 D. J. B. Smith1, M. J. Hardcastle1, S. Das1, A. Drake1, K.
Duncan2, G. G$\ddot{\rm{u}}$rkan 1,3,4 E-mail<EMAIL_ADDRESS>M.
Magliocchetti5, L. K. Morabito6,7, J. W. Petley6, S. Shenoy1 and C. Tasse8,9
1Centre for Astrophysics Research, University of Hertfordshire, College Lane,
Hatfield AL10 9AB, UK
2Institute for Astronomy, Royal Observatory, Blackford Hill, Edinburgh, EH9
3HJ, UK
3Astronomy, ATNF, PO Box 1130, Bentley, WA 6102, Australia
4Thüringer Landessternwarte, Sternwarte 5, 07778 Tautenburg, Germany
5INAF - IAPS, Via Fosso del Cavaliere 100, 00133, Rome, Italy
6Centre for Extragalactic Astronomy, Department of Physics, Durham University,
Durham DH1 3LE, UK
7Institute for Computational Cosmology, Department of Physics, University of
Durham, South Road, Durham DH1 3LE, UK
8GEPI $\&$ ORN, Observatoire de Paris, Université PSL, CNRS, 5 Place Jules
Janssen, 92190 Meudon, France
9Department of Physics $\&$ Electronics, Rhodes University, PO Box 94,
Grahamstown, 6140, South Africa

(Accepted XXX. Received YYY; in original form ZZZ)

###### Abstract

We use new 144 MHz observations over 5634 deg2 from the LOFAR 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_{\rm{1.4GHz}}/L_{i})$, we identify 3,697 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 QSO radio-loud. 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 AGN 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 radio-loud sources
either preferably reside in gas-rich systems with rapidly-spinning black
holes, or represent an earlier obscured phase of QSO evolution.

###### keywords:

quasars: general – galaxies: active – radio continuum: galaxies – techniques:
spectroscopic

††pubyear: 2024††pagerange: Exploring the radio-loudness of SDSS quasars with
spectral stacking–17

##  1 Introduction

The most luminous manifestations of active galactic nuclei (AGN) are quasi-
stellar objects (QSOs), also known as quasars, whose bolometric luminosity can
reach up to 1047-48 erg s-1 (e.g. Rakshit et al. 2020; Shen et al. 2020).
About 5-10$\%$ 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$\%$ are only weak radio sources, whose
emission could be purely a result of star formation (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
et al. 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 associated with star formation
through free-free emission from H${\rm{II}}$ regions and synchrotron radiation
from electrons accelerated to relativistic speeds by supernova remnants
(Condon, 1992). This leads to the question of whether star formation 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 1994; Sikora
et al. 2007; Tchekhovskoy et al. 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 et al. 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; Pâris et al. 2018) with radio observations from the
first data release of the LOFAR Two-metre Sky Survey (LoTSS; Shimwell et al.
2017) over the HETDEX spring field (Hill et al., 2008) and the LOFAR
H-ATLAS/NGP survey (Hardcastle et al., 2016). With the high sensitivity and
wide areal 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 $\sim$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 square degrees, LoTSS DR2 provides a much larger
observational volume, resulting in a substantial increase of RL quasars.
Employing our stacking techniques on this expanded dataset allows us to
systematically create composite spectra for each radio class in a given
parameter regime (i.e. redshift, $i$-band luminosity and black hole 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 black hole 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 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 star formation. 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 signal-to-noise (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$CDM cosmology with
$\Omega_{\Lambda}$ = 0.7, $\Omega_{M}$ = 0.3 and $H_{0}$ = 70 km s-1 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^{2}$ as shown in Figure 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
targetting 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 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 and 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_{\mathrm{bol}}$), derived virial black
hole mass ($M_{\mathrm{BH}}$) and Eddington ratio ($R_{\mathrm{edd}}$)
estimates. This was done using the publicly available spectral fitting code
PyQSOFit (Guo et al., 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
$\leq z$ < 1.9) and $L_{1350}$ ($z\geq$ 1.9) using the bolometric corrections
from Richards et al. (2006). For the virial black hole 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\geq 1.9$) using the calibrations from Vestergaard & Peterson
(2006), and the Mg ii line ($0.8\leq 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 second data release of the
LOFAR Two-metre Sky Survey (LoTSS DR2; Shimwell et al. 2022). This data
release covers 5634 deg2 of the northern sky (see footprint in Figure 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 summarise,
for smaller, more compact radio sources, the likelihood ratio method is
applied (e.g. Richter 1975; De Ruiter et al. 1977; Sutherland & Saunders
1992). This statistical approach uses both colour and magnitude information as
described by Williams et al. (2019) in order to maximise 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 multi-wavelength images. Finally, for everything in between a
decision tree is used to combine the two methods.

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.

###  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 (MOC) map as seen in Figure 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_{\mathrm{BH}}$ and $L_{\mathrm{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 LF
measurements) which require a uniform selection, and so is ideal for our
purposes.

To separate quasars into radio-loud and radio-quiet, we adopt the standard
definition of the radio-loudness parameter ($R$): the logarithm of the ratio
of 1.4 GHz radio luminosity ($L_{1.4\mathrm{GHz}}$) 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.4\mathrm{GHz}}$ for radio-detected
QSOs, we use the integrated 144 MHz flux density ($S_{144\mathrm{MHz}}$) from
the LoTSS DR2 catalogue, a radio spectral index of
$\alpha_{\mathrm{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.4\mathrm{GHz}}$, where $S_{144\mathrm{MHz}}$ 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_{\mathrm{opt}})\log(1+z).$ (1)

Here $z=2$ and we adopt an optical spectral index of $\alpha_{\mathrm{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 3,697 RL and 111,132 RQ quasars. The remaining sources are
discarded as their $5\sigma$ upper limit puts them in the radio-loud 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 Figure 2 (where we have also included a second
horizontal axis to show $R$ in terms of the $L_{\rm{144MHz}}$ 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
Figure 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.

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 radio-loud and radio-quiet QSOs, whereas the dashed black
lines are additional divisions we make to define classes as a function of $R$
as denoted.

###  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_{i}$ and $z$, $M_{i}$, $M_{\rm{BH}}$ 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_{i}$ 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 normalise the RL
2D-histogram by the total number of RL sources, and multiply the normalised
counts by the maximum number of RQ QSOs for which we can say that they are
drawn from the same $z$ and $M_{i}$ 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 $\alpha=0.05$. This
results in a sample of 77,532 RQ QSOs matched in $z$ and $M_{i}$ to our RL
sample (hereafter the 2D-matched or $z-M_{i}$ sample).

Similarly, the 3D matching process to create a sample matched in $z$, $M_{i}$
and $M_{\rm{BH}}$ (hereafter the 3D-matched or $z-M_{i}-M_{\rm{BH}}$ sample)
involves creating a normalised 3D-histogram of the RL population, and again
using multiple K-S tests to determine the maximum number of RQ sources. The
effectiveness of this matching procedure can be seen in Figure 3, where we
compare the individual $z$, $M_{i}$ and $M_{\mathrm{BH}}$ 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

Figure 3: A corner plot summary of the redshift, absolute $i$-band magnitude
and black hole mass distributions for the RL and RQ QSOs. The density contours
represent the 10, 30, 50, 70 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_{i}$ and $M_{\rm{BH}}$
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. 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
normalised, resampled spectra that went in to the stack. The labelled dotted
blue lines indicate prominent emission lines visible in the stacked spectrum.

###  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_{V}=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 value of the rest-frame
wavelengths sampled, along with the choice of sampling input in the algorithm
(e.g. 1 $\text{\AA}/\text{pixel}$). 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 normalise 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 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 normalised 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 normalised 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 realisations with the same size as the original sample. The
median method is then performed for each realisation to estimate the flux
density distribution at a given bin and thus determine the standard deviation.

Figure 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 (MPA-JHU; Brinchmann et al. 2004; Kauffmann et al. 2003 and
Tremonti et al. 2004). This galaxy template has a broad wavelength range
extending from 90 Å to 160 $\mu$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 Figure 13 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 normalisation 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 normalised
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_{\mathrm{input-
sim}}=\frac{F_{\lambda,\mathrm{input}}-F_{\lambda,\mathrm{sim}}}{\sigma_{\lambda,\mathrm{sim}}},$
(2)

where $F_{\lambda,\mathrm{input}}$ and $F_{\lambda,\mathrm{sim}}$ are the flux
density of the input and simulated composite, respectively and
${\sigma_{\lambda,\mathrm{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_{\mathrm{input-sim}}$ into 30
wavelength slices (or less, depending on pixel availability) and for each one
we fit a normal distribution. The best-fit parameters are then interpolated
with a spline to produce wavelength dependent corrections as indicated in
Figure 5. The best-fit of the mean is used to correct for the offset as a
result of the normalisation 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.

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 panel 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.

###  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_{i}$ 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 normalisation process, which requires a common wavelength range
free from strong spectral features. To account for this, we select all sources
with $-26<M_{i}<-24$ and divide them in redshift bins of size $\Delta 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 Figure 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 radio-loud) 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 (FOCAS; 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 (LRIS;
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 and it is therefore likely that the flux
density has a normalisation bias with uncertainties systematically
underestimated. However, creating and using such a template is beneficial as
this is the largest sample of radio-loud quasars at $z>4.9$ and its use
demonstrates the flexibility of our stacking algorithm in dealing with non-
SDSS spectra.

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 $\pm$ the uncertainty. The upper panel presents the full
wavelength coverage, where each composite is normalised to the RQ composite at
$3000-3600\text{\AA}$. The lower left panel presents a zoom-in window
indicating the differences bluewards of Ly$\alpha$, whereas the lower right
panel - the discrepancy above 5000Å which is likely caused by host galaxy
contamination.

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 significant host contamination (e.g. Glikman et al. 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$\alpha$,
on the other hand, can be explained by the different 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 Figure 6, where no IGM correction is applied as well.
Interestingly, the RQ composite appears to be less affected by the
Ly$\alpha$-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.

Figure 7: A comparison between the RL and RQ population across ten 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.

##  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 signal-to-noise (S/N) composite spectra for the
classically defined RL and RQ QSOs using the $z-M_{i}$ and
$z-M_{i}-M_{\rm{BH}}$ samples to establish whether there are any physical
differences between them and investigate the potential influence of the black
hole 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 section 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_{\rm{BH}}$ from Rakshit et al. (2020) and contain comparable
numbers of RL sources.

The resulting comparisons for all ten redshift bins are presented in Figure 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 sub-sample 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 ten redshift bins for
which we find a "reduced chi-squared value" of

$\chi^{2}_{\nu}=\frac{\sum_{i=1}^{N}\frac{(F_{RL,i}-F_{RQ,i})^{2}}{(\sigma_{RL,i}^{2}+\sigma_{RQ,i}^{2})}}{N},$
(3)

where N is the number of pixels in each stacked spectrum. Having done this
1000 times, we obtain a $\chi^{2}_{\nu}$ 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 Figure 7 and in 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.

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 radio-loud bins ($R_{2}$, $R_{3}$ and $R_{4}$) and the radio-quietest sample ($R_{1}$), and the radio-detected ($R_{1\rm{D}}$) and the radio-undetected ($R_{1\rm{U}}$) sample for the radio-quietest bin defined in section 4.3 per a given redshift bin. Redshift Range | (RL - RQ)2D | (RL - RQ)3D |  $R_{1}$ \- $R_{2}$ |  $R_{1}$ \- $R_{3}$ |  $R_{1}$ \- $R_{4}$ |  $R_{1\rm{U}}$ \- $R_{\rm{1D}}$
---|---|---|---|---|---|---
0.6–0.8 | 0.01 | 0.03 | $5.51\times 10^{-41}$ | $7.90\times 10^{-15}$ | $2.26\times 10^{-6}$ | 0.01
0.8–1.1 | $3.39\times 10^{-18}$ | $1.48\times 10^{-14}$
1.1–1.4 | $2.31\times 10^{-53}$ | $3.014\times 10^{-51}$ | $4.73\times 10^{-45}$ | $1.43\times 10^{-18}$ | $2.16\times 10^{-18}$ | 0.31
1.4–1.7 | $2.83\times 10^{-26}$ | $1.88\times 10^{-25}$
1.7–1.9 | $2.88\times 10^{-45}$ | $3.40\times 10^{-28}$ | $2.38\times 10^{-38}$ | $1.92\times 10^{-23}$ | $1.37\times 10^{-17}$ | 0.05
1.9–2.1 | $7.43\times 10^{-26}$ | $1.67\times 10^{-13}$
2.1–2.4 | $1.03\times 10^{-39}$ | $5.15\times 10^{-24}$ | $6.61\times 10^{-30}$ | $4.59\times 10^{-29}$ | $1.25\times 10^{-11}$ | 0.08
2.4–2.7 | $1.45\times 10^{-36}$ | $2.02\times 10^{-24}$
2.7–3.0 | $3.21\times 10^{-24}$ | $6.22\times 10^{-23}$ | 0.001 | $3.28\times 10^{-9}$ | $5.95\times 10^{-5}$ | 0.86
2.7–3.4 | $3.90\times 10^{-16}$ | $1.38\times 10^{-10}$

###  4.2 Black Hole Mass Dependence

As discussed in the introduction, several studies have found the black hole
mass to be a defining factor in the radio-loudness dichotomy, where RL QSOs
are found to harbour black holes with $M_{\rm{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 Figure 8 we present the
black hole mass distributions for each radio class using the virial black hole
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_{\rm{BH}}$ ranges,
irrespective of the different black hole mass calibrations. In addition, the
mean black hole 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 Figure 7, we re-do the stacking
procedure with the $z-M_{i}-M_{\rm{BH}}$ sample, i.e. including the MBH in the
process, even though it reduces the sample size. We note that by matching in
$M_{i}$ and $M_{\rm{BH}}$ that we are also effectively controlling for the
Eddington ratio.

Figure 8: The black hole mass estimates provided from Rakshit et al. (2020)
for the radio-loud (dark blue) and radio-quiet 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$\leq z$ < 1.9) and the
C IV line ($z\geq$ 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.

The resulting RL/RQ comparison is found to exhibit similar features as for the
$z-M_{i}$ 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_{i}-M_{\rm{BH}}$ 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_{i}$ 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 black-hole mass and accretion rate is required to explain the difference
between the RL and RQ population. As matching in $M_{\rm{BH}}$ is on average
of little consequence (except for sample size), we proceed with the rest of
our analysis using only the $z-M_{i}$ 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
population represent a bimodality or a transition between populations, we
further separate our sample into four radio-loudness bins as indicated in
Figure 2. As with our previous approach, we include radio-undetected sources
only in the lowest radio-loudness bin ($R_{1}$) 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 Figure 2). This leads to
some challenges when matching in $z$ and $M_{i}$ as demonstrated by Figure 9,
which presents the redshift-luminosity plane. We can see that the more radio-
loud bins ($R_{2}$, $R_{3}$ and $R_{4}$) in comparison to $R_{1}$ 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 maximise 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 signal-to-noise levels of the stacked spectra as before.

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.

The results of applying our spectral stacking algorithm and the null
hypothesis test on this sample division for five redshift bins show that all
radio-loud bins exhibit statistically significant differences when compared to
$R_{1}$ (see results for $R_{1}$-$R_{2}/R_{3}/R_{4}$ 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 radio-loud bins, we further conducted a comparison
between radio-detected ($R_{\rm{1D}}$) and radio-undetected sources
($R_{\rm{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\geq 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 sub-arcsecond 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 radio-loud and radio-quiet 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_{\rm{conti}}=f_{\rm{pl}}+f_{\rm{BC}}+f_{\rm{Fe\,\textsc{ii}}}+f_{\rm{poly}}.$
(4)

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

$f_{\rm{pl}}=a_{0}(\lambda/\lambda_{0})^{a_{1}},$ (5)

where $a_{0}$ and $a_{1}$ are the normalisation 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\text{\AA}$. In PyQSOFit it is modelled by the function given by Grandi
(1982) for the case of optically thick clouds as:

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

where $\lambda_{\rm{BE}}$ is the position of the Balmer edge, $F_{\rm{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_{\rm{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_{\rm{Fe\,\textsc{ii}}}=c_{0}F_{\rm{Fe\,\textsc{ii}}}(\lambda,c_{1},c_{2}),$
(7)

where $c_{0}$ is the normalisation constant, $c_{1}$ is a constant describing
the Gaussian broadening and $c_{2}$ represents the wavelength shift 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 et al. (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_{\rm{poly}}=\sum^{i=3}_{i=0}b_{i}(\lambda-\lambda_{0})^{i},$ (8)

where $b_{i}$ are the model free parameters.

After subtracting the best-fit 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-1) 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.

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. Emission Line | Wavelength | Emission | Number of
---|---|---|---
Complex | Range | Line | Gaussian
Name | (Å) | Name | components
H${\beta}$ | 4640-5100 | H${\beta}$ broad | 2
|  | H${\beta}$ narrow | 1
|  |  $[$O III$]$ (4959 Å) | 1
|  |  $[$O III$]$ (5007 Å) | 1
H${\gamma}$ | 4200-4440 | H${\gamma}$ broad | 1
|  | Hγ narrow | 1
|  | O III (4363 Å) | 1
$[$O ii$]$ | 3650-3800 |  $[$O ii$]$ (3728 Å) | 1
$[$Ne v$]$ | 3380-3480 |  $[$Ne v$]$ (3426 Å) | 1
Mg ii] | 2700-2900 | Mg ii] broad | 2
|  | Mg ii] narrow | 1
C iii] | 1810-1970 | C iii] broad | 2
|  | C iii] narrow | 1
C iv | 1500-1700 | C iv broad | 2
|  | C iv narrow | 1
|  | He ii (1640 Å) broad | 1
|  | He ii (1640 Å) narrow | 1
|  | O iii] (1663 Å) broad | 1
|  | O iii] (1663 Å) narrow | 1
Si iv \+ O iv] | 1290-1450 | Si iv \+ O iv] |
Ly$\alpha$ | 1150-1290 | Ly$\alpha$ broad | 2
|  | Ly$\alpha$ narrow | 1
|  | N v (1240 Å) | 1
O vi | 980-1100 | O vi | 1

###  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_{\nu}^{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 localised 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_{\rm{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-1 cm-2 $\text{\AA}^{-1}$.

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$\sigma$) 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.

Figure 10 presents the ratio of the emission line flux of the RL to RQ QSOs
(RL-RQ), along with the more radio-loud 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
radio-loud 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
ionised by the continuous radiation coming from the accretion disc. The [O ii]
emission line, on the other hand, is a well-known star formation 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 a [Ne v] detection and making use of the [O ii]-to-[Ne v] emission line
flux ratio, which in a typical AGN-dominated galaxy is found to be of order
unity as a result of excitation from NLR. The [Ne v] is chosen because of its
wavelength proximity to [O ii] and its high ionisation 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 [Ne v], we used PyQSOFit to fit a single Gaussian in
a spectral window of 150 Å centered on the [Ne v] central wavelength. The
emission line is considered as a non-detection if the measured emission line
flux is detected at a S/N<3\. Although in previous sections we investigated
how the optical properties varied in bins of radio-loudness, the need to
remove the [Ne v] 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_{i}$, with the division set at
$R=0.2$ (i.e. the separation between $R_{1}$ and the more radio-loud 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 [O ii] and [Ne v] in the same manner as before. A [Ne v]
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 [O ii] / [Ne v] = 1.05 obtained from the SDSS composite
spectrum from Vanden Berk et al. (2001), we are able to derive a measure of
the [O ii] emission line flux unaffected by AGN.

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 a SFR of
zero.

Figure 11 shows the [O ii] derived, doubly corrected SFR against redshift
obtained for the RL and RQ populations using the [O ii] - SFR conversion from
Kewley et al. (2004). We can see that the SFR of the RQ population is close to
zero, indicating that on average star-formation processes in the host galaxy
of RQ QSOs are not sufficient to produce a strong [O ii] emission line. In
contrast, the excess of [O ii] 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
(FIRST; Becker et al. 1995), but we have shown that their findings persist
with higher statistical significance and to a further $\approx 1.4$ Gyr of
cosmic history, closer to the peak epoch of cosmic star formation rate 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).

##  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 star
formation 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 radio-quiet
quasars hosted by starburst galaxies in our RL sample. The presence of intense
star formation in such hosts can lead to a larger degree of reddening (Gordon
et al. 1997), which would further explain the redder continuum observed in
section 4.1. In addition, this scenario does not require any differences
between the black hole mass and accretion rate of the QSOs themselves (i.e.
$R_{\rm{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 a SFR value for a typical starburst galaxy
with a stellar mass of $10^{11}$ M⊙ at a redshift equal to our central bin
values from Figure 11. This resulted in SFR values ranging from 127 - 488 M⊙
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_{150\rm{MHz}}$ 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 [Ne v] emission line to remove contamination
from EELR and correct for contribution from the AGN NLR. However, shocks in
the interstellar medium generated 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 km
s-1) are able to produce [Ne v]. 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 MHz to 178 MHz, in order to compare our results with
Hardcastle et al. (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 Figure 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.

Figure 12: The [O ii] emission line luminosity as a function of 178MHz total
radio luminosity for RL sources without significant [Ne v] 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.

###  6.3 Possible explanations

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

The first model involves a rapidly spinning black hole, coupled with a rich
gas supply required for accretion and jet production (e.g. Hardcastle et al.
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 fueling star
formation, these two processes are inherently interconnected. Therefore, the
question becomes how to get a rapidly spinning black hole and a rich gas
supply to make a RL quasar. We suggest two plausible scenarios: Firstly, 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 black hole, but also serves as a
trigger for star formation. 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 star
formation 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 black hole, 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 e.g. 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; Fawcett et al. 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 radio-quiet (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 sub-arcsecond 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 (DESI; Aghamousa et al.
2016a, b) and the Multi-object Optical and Near-IR spectrograph (MOONs;
Cirasuolo et al. 2014). In addition, improved black hole 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 (NuSTAR; 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 second data release of the LOFAR Two-metre Sky
Survey and the fourteenth data release of the Sloan Digital Sky Survey to
create the largest, uniformly-selected, spectroscopically-confirmed sample of
radio-loud and radio-quiet 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:

  * -

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 [O ii]
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.

  * -

The RL and RQ population are found to span similar black hole mass ranges with
no systematic trend showing that the mean $M_{\rm{BH}}$ is higher for RL QSOs.
Furthermore, using a sample matched in $z$, $M_{i}$ and $M_{\rm{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.

  * -

The observed differences between the RL and RQ population are not gradual.
Comparing composite spectra as a function of radio-loudness shows that all
more radio-loud bins ($R_{2}$, $R_{3}$, ${R_{4}}$) differ from the radio-
quietest ($R_{1}$), with features consistent with the classical RL and RQ
division.

  * -

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.

  * -

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 [O ii] 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 black hole mass or the accretion rate. As a result, we propose two
distinct models: one requiring that RL quasars have rapidly spinning black
holes 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 radio-loud 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 number 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 number ST/W507490/1. SS acknowledges support from the STFC
via studentship grant number 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 benefited 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 Sloan Digital Sky Survey website
(www.sdss.org.)

## References

  * Aghamousa et al. (2016a) Aghamousa A., et al., 2016a, arXiv preprint arXiv:1611.00036
  * Aghamousa et al. (2016b) Aghamousa A., et al., 2016b, arXiv preprint arXiv:1611.00037
  * Alexander & Hickox (2012) Alexander D. M., Hickox R. C., 2012, New Astronomy Reviews, 56, 93
  * Allen et al. (2008) Allen M. G., Groves B. A., Dopita M. A., Sutherland R. S., Kewley L. J., 2008, The Astrophysical Journal Supplement Series, 178, 20
  * Baloković et al. (2012) Baloković M., Smolčić V., Ivezić Ž., Zamorani G., Schinnerer E., Kelly B. C., 2012, The Astrophysical Journal, 759, 30
  * Barcons et al. (2017) Barcons X., et al., 2017, Astronomische Nachrichten, 338, 153
  * Becker et al. (1995) Becker R. H., White R. L., Helfand D. J., 1995, The Astrophysical Journal, 450, 559
  * Best & Heckman (2012) Best P., Heckman T., 2012, Monthly Notices of the Royal Astronomical Society, 421, 1569
  * Best et al. (2005) Best P., Kauffmann G., Heckman T., Brinchmann J., Charlot S., Ivezić Ž., White S., 2005, Monthly Notices of the Royal Astronomical Society, 362, 25
  * Blandford & Znajek (1977) Blandford R. D., Znajek R. L., 1977, Monthly Notices of the Royal Astronomical Society, 179, 433
  * Bovy et al. (2011) Bovy J., et al., 2011, The Astrophysical Journal, 729, 141
  * Brinchmann et al. (2004) Brinchmann J., Charlot S., White S., Tremonti C., Kauffmann G., Heckman T., Brinkmann J., 2004, Monthly Notices of the Royal Astronomical Society, 351, 1151
  * Bruzual & Charlot (2003) Bruzual G., Charlot S., 2003, Monthly Notices of the Royal Astronomical Society, 344, 1000
  * Calabrò et al. (2021) Calabrò A., et al., 2021, Astronomy & Astrophysics, 646, A39
  * Carnall (2017) Carnall A., 2017, arXiv preprint arXiv:1705.05165
  * Chiaberge & Marconi (2011) Chiaberge M., Marconi A., 2011, Monthly Notices of the Royal Astronomical Society, 416, 917
  * Chonis et al. (2016) Chonis T. S., et al., 2016, in Ground-based and Airborne Instrumentation for Astronomy VI. pp 1345–1372
  * Cirasuolo et al. (2003a) Cirasuolo M., Magliocchetti M., Celotti A., Danese L., 2003a, Monthly Notices of the Royal Astronomical Society, 341, 993
  * Cirasuolo et al. (2003b) Cirasuolo M., Celotti A., Magliocchetti M., Danese L., 2003b, Monthly Notices of the Royal Astronomical Society, 346, 447
  * Cirasuolo et al. (2014) Cirasuolo M., et al., 2014, in Ground-based and airborne instrumentation for astronomy V. pp 202–214
  * Condon (1992) Condon J., 1992, Annual review of astronomy and astrophysics, 30, 575
  * Condon et al. (2013) Condon J., Kellermann K., Kimball A. E., Ivezić Ž., Perley R., 2013, The Astrophysical Journal, 768, 37
  * Dalton et al. (2012) Dalton G., et al., 2012, in Ground-based and Airborne Instrumentation for Astronomy IV. pp 220–231
  * De Ruiter et al. (1977) De Ruiter H., Willis A., Arp H., 1977, Astronomy and Astrophysics Supplement Series, 28, 211
  * Dunlop et al. (2003) Dunlop J., McLure R., Kukula M., Baum S., O’Dea C., Hughes D., 2003, Monthly Notices of the Royal Astronomical Society, 340, 1095
  * Fawcett et al. (2020) Fawcett V., Alexander D., Rosario D., Klindt L., Fotopoulou S., Lusso E., Morabito L., Calistro Rivera G., 2020, Monthly Notices of the Royal Astronomical Society, 494, 4802
  * Fawcett et al. (2022) Fawcett V., Alexander D., Rosario D., Klindt L., Lusso E., Morabito L., Calistro Rivera G., 2022, Monthly Notices of the Royal Astronomical Society, 513, 1254
  * Fitzpatrick (1999) Fitzpatrick E. L., 1999, Publications of the Astronomical Society of the Pacific, 111, 63
  * Fynbo et al. (2012) Fynbo J., Krogager J.-K., Venemans B., Noterdaeme P., Vestergaard M., Møller P., Ledoux C., Geier S., 2012, The Astrophysical Journal Supplement Series, 204, 6
  * Glikman et al. (2006) Glikman E., Helfand D. J., White R. L., 2006, The Astrophysical Journal, 640, 579
  * Gloudemans et al. (2022) Gloudemans A., et al., 2022, Astronomy & Astrophysics, 668, A27
  * Gordon et al. (1997) Gordon K. D., Calzetti D., Witt A. N., 1997, The Astrophysical Journal, 487, 625
  * Grandi (1982) Grandi S., 1982, The Astrophysical Journal, 255, 25
  * Gu et al. (2001) Gu M., Cao X., Jiang D., 2001, Monthly Notices of the Royal Astronomical Society, 327, 1111
  * Guo et al. (2018) Guo H., Shen Y., Wang S., 2018, Astrophysics Source Code Library, pp ascl–1809
  * Gürkan et al. (2015) Gürkan G., et al., 2015, Monthly Notices of the Royal Astronomical Society, 452, 3776
  * Gürkan et al. (2019) Gürkan G., et al., 2019, Astronomy & Astrophysics, 622, A11
  * Hardcastle et al. (2007) Hardcastle M. J., Evans D., Croston J., 2007, Monthly Notices of the Royal Astronomical Society, 376, 1849
  * Hardcastle et al. (2009) Hardcastle M., Evans D., Croston J., 2009, Monthly Notices of the Royal Astronomical Society, 396, 1929
  * Hardcastle et al. (2016) Hardcastle M., et al., 2016, Monthly Notices of the Royal Astronomical Society, 462, 1910
  * Hardcastle et al. (2023) Hardcastle M., et al., 2023, arXiv preprint arXiv:2309.00102
  * Harris et al. (2016) Harris D. W., et al., 2016, The Astronomical Journal, 151, 155
  * Harrison et al. (2010) Harrison F. A., et al., 2010, in Space Telescopes and Instrumentation 2010: Ultraviolet to Gamma Ray. pp 189–196
  * Heckman et al. (2004) Heckman T. M., Kauffmann G., Brinchmann J., Charlot S., Tremonti C., White S. D., 2004, The Astrophysical Journal, 613, 109
  * Hill et al. (2008) Hill G. J., et al., 2008, arXiv preprint arXiv:0806.0183
  * Ho (2002) Ho L. C., 2002, The Astrophysical Journal, 564, 120
  * Ho (2005) Ho L. C., 2005, The Astrophysical Journal, 629, 680
  * Hopkins et al. (2008) Hopkins P. F., Hernquist L., Cox T. J., Kereš D., 2008, The Astrophysical Journal Supplement Series, 175, 356
  * Husemann et al. (2014) Husemann B., Jahnke K., Sánchez S., Wisotzki L., Nugroho D., Kupko D., Schramm M., 2014, Monthly Notices of the Royal Astronomical Society, 443, 755
  * Ivezić et al. (2002) Ivezić Ž., et al., 2002, The Astronomical Journal, 124, 2364
  * Jin et al. (2023) Jin S., et al., 2023, Monthly Notices of the Royal Astronomical Society, p. stad557
  * Kalfountzou et al. (2012) Kalfountzou E., Jarvis M., Bonfield D., Hardcastle M., 2012, Monthly Notices of the Royal Astronomical Society, 427, 2401
  * Kashikawa et al. (2002) Kashikawa N., et al., 2002, Publications of the Astronomical Society of Japan, 54, 819
  * Kauffmann et al. (2003) Kauffmann G., et al., 2003, Monthly Notices of the Royal Astronomical Society, 341, 33
  * Kellermann et al. (1989) Kellermann K., Sramek R., Schmidt M., Shaffer D., Green R., 1989, The Astronomical Journal, 98, 1195
  * Kellermann et al. (2016) Kellermann K., Condon J., Kimball A., Perley R., Ivezić Ž., 2016, The Astrophysical Journal, 831, 168
  * Kewley et al. (2004) Kewley L. J., Geller M. J., Jansen R. A., 2004, The Astronomical Journal, 127, 2002
  * Kewley et al. (2006) Kewley L. J., Groves B., Kauffmann G., Heckman T., 2006, Monthly Notices of the Royal Astronomical Society, 372, 961
  * Kimball et al. (2011) Kimball A. E., Kellermann K. I., Condon J. J., Ivezić Ž., Perley R. A., 2011, The Astrophysical Journal Letters, 739, L29
  * Klindt et al. (2019) Klindt L., Alexander D. M., Rosario D. J., Lusso E., Fotopoulou S., 2019, Monthly Notices of the Royal Astronomical Society, 488, 3109
  * Kondapally et al. (2021) Kondapally R., et al., 2021, Astronomy & Astrophysics, 648, A3
  * Laor & Behar (2008) Laor A., Behar E., 2008, Monthly Notices of the Royal Astronomical Society, 390, 847
  * Liu et al. (2006) Liu Y., Jiang D. R., Gu M. F., 2006, The Astrophysical Journal, 637, 669
  * Macfarlane et al. (2021) Macfarlane C., et al., 2021, Monthly Notices of the Royal Astronomical Society, 506, 5888
  * Madau & Dickinson (2014) Madau P., Dickinson M., 2014, Annual Review of Astronomy and Astrophysics, 52, 415
  * Maddox (2018) Maddox N., 2018, Monthly Notices of the Royal Astronomical Society, 480, 5203
  * Magliocchetti et al. (2020) Magliocchetti M., et al., 2020, Monthly Notices of the Royal Astronomical Society, 493, 3838
  * Matsuoka et al. (2015) Matsuoka Y., et al., 2015, The Astrophysical Journal, 811, 91
  * McLure & Jarvis (2004) McLure R. J., Jarvis M. J., 2004, Monthly Notices of the Royal Astronomical Society, 353, L45
  * Metcalf & Magliocchetti (2006) Metcalf R., Magliocchetti M., 2006, Monthly Notices of the Royal Astronomical Society, 365, 101
  * Morabito et al. (2022) Morabito L. K., et al., 2022, Monthly Notices of the Royal Astronomical Society, 515, 5758
  * Myers et al. (2015) Myers A. D., et al., 2015, The Astrophysical Journal Supplement Series, 221, 27
  * Oke et al. (1995) Oke J., et al., 1995, Publications of the Astronomical Society of the Pacific, 107, 375
  * Pâris et al. (2018) Pâris I., et al., 2018, Astronomy & Astrophysics, 613, A51
  * Radcliffe et al. (2021) Radcliffe J. F., Barthel P., Garrett M., Beswick R., Thomson A., Muxlow T., 2021, Astronomy & Astrophysics, 649, L9
  * Rakshit et al. (2020) Rakshit S., Stalin C., Kotilainen J., 2020, The Astrophysical Journal Supplement Series, 249, 17
  * Richards et al. (2002) Richards G. T., et al., 2002, The Astronomical Journal, 123, 2945
  * Richards et al. (2006) Richards G. T., et al., 2006, The Astrophysical Journal Supplement Series, 166, 470
  * Richter (1975) Richter G., 1975, Astronomische Nachrichten, 296, 65
  * Rigby et al. (2018) Rigby J., et al., 2018, The Astrophysical Journal, 853, 87
  * Rodighiero et al. (2011) Rodighiero G., et al., 2011, The Astrophysical journal letters, 739, L40
  * Ross et al. (2012) Ross N. P., et al., 2012, The Astrophysical Journal Supplement Series, 199, 3
  * Rowlands et al. (2012) Rowlands K., et al., 2012, Monthly Notices of the Royal Astronomical Society, 419, 2545
  * Salviander et al. (2007) Salviander S., Shields G. A., Gebhardt K., Bonning E., 2007, The Astrophysical Journal, 662, 131
  * Sanders et al. (1990) Sanders D., Soifer B., Neugebauer G., 1990, in Dynamics and Interactions of Galaxies: Proceedings of the International Conference, Heidelberg, 29 May–2 June 1989. pp 459–464
  * Schlegel et al. (1998) Schlegel D. J., Finkbeiner D. P., Davis M., 1998, The Astrophysical Journal, 500, 525
  * Schneider et al. (2010) Schneider D. P., et al., 2010, The Astronomical Journal, 139, 2360
  * Schreiber et al. (2015) Schreiber C., et al., 2015, Astronomy & Astrophysics, 575, A74
  * Selsing et al. (2016) Selsing J., Fynbo J. P., Christensen L., Krogager J.-K., 2016, Astronomy & Astrophysics, 585, A87
  * Shankar et al. (2010) Shankar F., Sivakoff G. R., Vestergaard M., Dai X., 2010, Monthly Notices of the Royal Astronomical Society, 401, 1869
  * Shen et al. (2019) Shen Y., et al., 2019, The Astrophysical Journal Supplement Series, 241, 34
  * Shen et al. (2020) Shen X., Hopkins P. F., Faucher-Giguère C.-A., Alexander D., Richards G. T., Ross N. P., Hickox R., 2020, Monthly Notices of the Royal Astronomical Society, 495, 3252
  * Shimwell et al. (2017) Shimwell T., et al., 2017, Astronomy & Astrophysics, 598, A104
  * Shimwell et al. (2022) Shimwell T., et al., 2022, Astronomy & Astrophysics, 659, A1
  * Sikora et al. (2007) Sikora M., Stawarz Ł., Lasota J.-P., 2007, The Astrophysical Journal, 658, 815
  * Smee et al. (2013) Smee S. A., et al., 2013, The Astronomical Journal, 146, 32
  * Smith et al. (2016) Smith D., et al., 2016, arXiv preprint arXiv:1611.02706
  * Smith et al. (2021) Smith D., et al., 2021, Astronomy & Astrophysics, 648, A6
  * Sutherland & Saunders (1992) Sutherland W., Saunders W., 1992, Monthly Notices of the Royal Astronomical Society, 259, 413
  * Symeonidis et al. (2016) Symeonidis M., Giblin B., Page M., Pearson C., Bendo G., Seymour N., Oliver S., 2016, Monthly Notices of the Royal Astronomical Society, 459, 257
  * Tchekhovskoy et al. (2010) Tchekhovskoy A., Narayan R., McKinney J. C., 2010, The Astrophysical Journal, 711, 50
  * Tremonti et al. (2004) Tremonti C. A., et al., 2004, The Astrophysical Journal, 613, 898
  * Tsuzuki et al. (2006) Tsuzuki Y., Kawara K., Yoshii Y., Oyabu S., Tanabé T., Matsuoka Y., 2006, The Astrophysical Journal, 650, 57
  * Unger et al. (1987) Unger S., Pedlar A., Axon D. J., Whittle M., Meurs E., Ward M., 1987, Monthly Notices of the Royal Astronomical Society, 228, 671
  * Urry & Padovani (1995) Urry C. M., Padovani P., 1995, Publications of the Astronomical Society of the Pacific, 107, 803
  * Vanden Berk et al. (2001) Vanden Berk D. E., et al., 2001, The Astronomical Journal, 122, 549
  * Véron-Cetty et al. (2004) Véron-Cetty M.-P., Joly M., Véron P., 2004, Astronomy & Astrophysics, 417, 515
  * Vestergaard & Osmer (2009) Vestergaard M., Osmer P. S., 2009, The Astrophysical Journal, 699, 800
  * Vestergaard & Peterson (2006) Vestergaard M., Peterson B. M., 2006, The Astrophysical Journal, 641, 689
  * Vestergaard & Wilkes (2001) Vestergaard M., Wilkes B. J., 2001, The Astrophysical Journal Supplement Series, 134, 1
  * Villar-Martín et al. (2011) Villar-Martín M., Tadhunter C., Humphrey A., Encina R. F., Delgado R. G., Torres M. P., Martínez-Sansigre A., 2011, Monthly Notices of the Royal Astronomical Society, 416, 262
  * White et al. (2007) White R. L., Helfand D. J., Becker R. H., Glikman E., de Vries W., 2007, The Astrophysical Journal, 654, 99
  * White et al. (2015) White S. V., Jarvis M. J., Häußler B., Maddox N., 2015, Monthly Notices of the Royal Astronomical Society, 448, 2665
  * White et al. (2017) White S. V., Jarvis M. J., Kalfountzou E., Hardcastle M. J., Verma A., Cao Orjales J. M., Stevens J., 2017, Monthly Notices of the Royal Astronomical Society, 468, 217
  * Williams et al. (2019) Williams W., et al., 2019, Astronomy & Astrophysics, 622, A2
  * Wilson & Colbert (1994) Wilson A. S., Colbert E. J., 1994, arXiv preprint astro-ph/9408005
  * Yip et al. (2004a) Yip C.-W., et al., 2004a, The Astronomical Journal, 128, 585
  * Yip et al. (2004b) Yip C., et al., 2004b, The Astronomical Journal, 128, 2603
  * Yue et al. (2018) Yue M., et al., 2018, The Astrophysical Journal, 863, 21
  * Zakamska et al. (2016) Zakamska N. L., et al., 2016, Monthly Notices of the Royal Astronomical Society, 455, 4191
  * Zhu et al. (2015) Zhu G. B., et al., 2015, The Astrophysical Journal, 815, 48

##  Appendix A Supplementary material

Figure 13: Results of applying our spectral stacking algorithm to sources
according to their location on the BPT-NII diagram. Left panel: The BPT-NII
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-NII diagram. Emission lines of interest
are labelled, and the legend indicates the number of spectra included in each
stack. Figure 14: 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 15: A
high S/N comparison between the RL and RQ population matched in $z$, $M_{i}$
and $M_{\rm{BH}}$. As in Figure 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 figure 14, $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 16: A high S/N
comparison between the radio-quietest ($R_{1}$) and more radio-loud ($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. Figure 17: 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 fit model (blue). The grey shaded region show the wavelength
windows used for the continuum fit. The lower panels show the best-fit model
of individual line complexes (blue), along with the decomposition into broad
(red) and narrow (green) line components.