Exploring the radio-loudness of SDSS quasars with spectral stacking

We use new 144 MHz observations over 5634 deg$^2$ 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 [OII] 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 [OII] emission line, with the contribution from AGN removed using the [NeV] 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.


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 10 47−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 (; the ratio of radio to ★ E-mail: m.arnaudova@herts.ac.uk 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 .For example, the use of optical and radio information leads to inhomogeneous samples as a result of different selection effects.In addition, the  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 HII 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 radioloudness parameter on galaxy properties such as redshift, bolometric luminosity, radio luminosity, BH mass and Eddington ratio.They found that quasars exhibit a wide continuum of radio properties, with no clear signatures of a bimodality.Given these results, the authors favoured the scenario where both AGN jets and SF contribute to the radio emission in quasars such that there is no RL/RQ dichotomy, but rather a smooth transition between the dominance of these two processes.
Recently, Macfarlane et al. (2021) built upon these results by developing a numerical model of the radio flux densities of quasars, in which the radio emission of every quasar consists of two components: AGNs (jets) and SF.This model, coupled with Monte Carlo simulations, allowed the authors to create quasar mock samples and compare them with observations.Their results were found to be in excellent agreement with the observed radio flux distributions of ∼42,000 SDSS quasars as measured in LoTSS DR1 across several redshift and absolute -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, -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 ΛCDM cosmology with Ω Λ = 0.7, Ω  = 0.3 and  0 = 70 km s −1 Mpc −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 -band absolute magnitudes brighter than   [ = 2] = −22.0over a wide redshift range of 0.065 <  < 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 ′′ and an average resolving power of Δ/ ≈ 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 ′′ , 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 (  [ = 2] = −20.5)at redshift ranges of 2.15 <  < 3.5 and 0.9 <  < 2.2, respectively (for further details see Ross et al. 2012 andMyers 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 valueadded catalogue by Rakshit et al. (2020) which includes continuum and line property measurements, including bolometric luminosity ( bol ), derived virial black hole mass ( BH ) and Eddington ratio ( 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 ( < 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  5100 ( < 0.7),  3000 (0.7 ≤  < 1.9) and  1350 ( ≥ 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  line ( < 0.8) and the C iv line ( ≥ 1.9) using the calibrations from Vestergaard & Peterson (2006), and the Mg ii line (0.8 ≤  < 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.

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 deg 2 of the northern sky (see footprint in Figure 1) at a resolution of 6 ′′ and a median rms sensitivity of 83 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.

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  BH and  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 homogeneouslyselected 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 (): the logarithm of the ratio of 1.4 GHz radio luminosity ( 1.4GHz ) to optical -band luminosity (  ), where  = 1 marks the boundary (e.g.Baloković et al. 2012).To calculate the -corrected  1.4GHz for radio-detected QSOs, we use the integrated 144 MHz flux density ( 144MHz ) from the LoTSS DR2 catalogue, a radio spectral index of  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 upper limit for  1.4GHz , where  144MHz is taken to be 5× 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 -band magnitude from the DR14Q catalogue which is -corrected to  = 2 and use the following conversion given by Richards et al. (2006) to obtain a -corrected estimate to  = 0: (1) Here  = 2 and we adopt an optical spectral index of  opt = -0.5.To obtain the -band luminosity, we use the simple relationship between luminosity and magnitude, where the solar luminosity (3.827 × 10 26 W) and solar -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 upper limit puts them in the radioloud 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  in terms of the  144MHz for comparison).It is evident that in both cases the conventional demarcation line ( = 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.

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 ,   and ,   ,  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  and   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  and   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 -value is greater than a significance level of  = 0.05.This results in a sample of 77,532 RQ QSOs matched in  and   to our RL sample (hereafter the 2D-matched or  −   sample).
Similarly, the 3D matching process to create a sample matched in ,   and  BH (hereafter the 3D-matched or  −   −  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 ,   and  BH distributions, along with contours of the 2D-distributions of our 3Dmatched sample of 23,716 RQ QSOs with the whole RL and RQ samples.

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 observedframe by using the re-calibrated reddening data,  ( − ), from Schlegel, Finkbeiner & Davis (1998), along with the Milky Way reddening curve from Fitzpatrick (1999) for an extinction-to-reddening ratio of   = 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 Å/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.

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 andTremonti et al. 2004).This galaxy template has a broad wavelength range extending from 90 Å to 160 m, a resolving power similar to that of SDSS across the range from 3200 to 9500 Å and contains the main nebular emission lines targeted by both WEAVE-LOFAR and SDSS.These features enable us to simulate mock samples of galaxy spectra based on SDSS characteristics across a wide range of redshifts, providing us with a robust way of testing the spectral stacking algorithm and further demonstrating its flexibility in handling different type of spectra (see Figure A1 for an example of stacking spectra as a function of BPT classes).
Following a range of tests including continuum recovery at low brightness and performance related to different sampling of spectra (the main motivations for implementing it for the WEAVE-LOFAR survey), we have identified two implications of the stacking method.First of all, it is not always possible for the algorithm to choose a 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 deredshifting 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, -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: where  ,input and  ,sim are the flux density of the input and simulated composite, respectively and  ,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  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.

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 <  < 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  −   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 <   < −24 and divide them in redshift bins of size Δ = 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- 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 <  < 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- 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  > 4.9 and its use demonstrates the flexibility of our stacking algorithm in dealing with non-SDSS spectra.
Despite the difference in selection, we find qualitatively a good agreement in the continuum shape with the various templates.The difference observed redwards of 5000Å amongst the composites is likely caused by various degrees of host galaxy contamination.The SDSS composite from Vanden Berk et al. (2001) includes intrinsically faint sources at low redshifts, making it subject to 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, 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- 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-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.

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  −   and  −   −  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.

High S/N Comparison of RL and RQ QSOs
Starting with our  −   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"  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  and   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 where N is the number of pixels in each stacked spectrum.Having  done this 1000 times, we obtain a  2  distribution under the null hypothesis for each redshift range.By fitting this distribution with its parametric form, we can estimate the probability of obtaining our observed results under the null hypothesis.
The results of the null hypothesis test for this 2D matched sample are presented in the upper left corner of each panel in Figure 7 and in Table 1.For all redshift bins, we find -values smaller than a significance level of  = 0.05, indicating that the two samples are not drawn from the same parent population.

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  BH ≳ 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  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 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  −   −  BH sample, i.e. including the M BH in the process, even though it reduces the sample size.We note that by matching in   and  BH that we are also effectively controlling for the Eddington ratio.
The resulting RL/RQ comparison is found to exhibit similar fea-   2020) for the radio-loud (dark blue) and radio-quiet population matched in redshift and -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≤  < 1.9) and the C IV line ( ≥ 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.tures as for the  −   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  −   −  BH sample).We use the null hypothesis test as described in the previous section to determine the significance of this result.Similarly to the  −   sample, we are able to rule out the null hypothesis with  < 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  BH is on average of little consequence (except for sample size), we proceed with the rest of our analysis using only the  −   sample, since it is about three times larger.

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 popu-lations, 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 ( 1 ) as we cannot confidently assign the rest of the sources to the correct  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  and   as demonstrated by Figure 9, which presents the redshiftluminosity plane.We can see that the more radio-loud bins ( 2 ,  3 and  4 ) in comparison to  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.
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  1 (see results for  1 - 2 / 3 / 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 ( 1D ) and radio-undetected sources ( 1U ) within the lowest radio-loudness bin ( 1 ).We find that we can accept the null hypothesis that the two samples are drawn from the same parent distribution ( ≥ 0.05) for all, apart from the lowest redshift bin, where  = 0.01 indicating ambiguity.This ambiguous result, however, is not enough to explain -values as low as 10 −41 as found for the  1 - 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  ⊙ , 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  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  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.

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.

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: (4) The power-law continuum component has been added to represent the emission from the accretion disc and is defined as: where  0 and  1 are the normalisation and power-law slope, and  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  ∼ 3000Å.In PyQSOFit it is modelled by the function given by Grandi (1982) for the case of optically thick clouds as: where  BE is the position of the Balmer edge,  BE is the flux at the Balmer edge,   (  ) is the Planck function at the electron temperature   and   is the optical depth.Here,  BE ,   and   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 pseudocontinuum.If not properly subtracted, this pseudo-continuum could contaminate the continuum and emission line measurements.Here, PyQSOFit models it as: where  0 is the normalisation constant,  1 is a constant describing the Gaussian broadening and  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)  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: where   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.

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  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 ( centre ) and the SDSS -band filter curve and the median apparent -band magnitude of the individual spectra to obtain a flux density in units of erg s −1 cm −2 Å −1 .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 radioquietest bin ( 2 / 3 / 4 - 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 ) 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.

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.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  and   , with the division set at  = 0.2 (i.e. the separation between  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 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 classifica- tion ( > 0.2) also exhibits the same trend, further indicating that the exact value of  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 ≈ 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).

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. 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 <  < 1.9.These estimates are then converted to a radio luminosity at 150 MHz by using the mass independent SFR- 150MHz relation from Smith et al. (2021).Comparing these results to the RL sample, we find that up to 5 per cent in any given redshift bin, could have a radio luminosity consistent with contribution from RQ starbursts.However, taking into account the large difference between our doubly corrected [O ii] SFR and the values obtained here for the starburst galaxies, it is unlikely that the radio emission is due to such extreme star-forming processes.Furthermore, considering the fact that we are investigating the median stacks of the populations, we believe that if starburst galaxies are present in our samples they will not noticeably impact our results.

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 moderatevelocity 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  = −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 <  < 0.8), and up to 13 per cent for the highest range (1.7 <  < 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.

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. 2020Fawcett et al. , 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  6 of Hardcastle et al. (2009), where the solid line indicates the regression line, whereas the dashed lines show ±3× the obtained scatter.
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 shortlived 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. < 1).It is therefore clear that this association alone cannot explain our results (certainly we cannot equate the RL QSOs with the 'red QSO' class, etc).
To make further progress on what controls the radio-loudness of QSOs, we need more information.This could come from the subarcsecond 144 MHz imaging that is now becoming possible with LOFAR (and the morphological information that it can provide; e.g.Morabito et al. 2022), along with larger statistical samples from new and forthcoming facilities such as WEAVE (Dalton et al. 2012;Jin et al. 2023), the Dark Energy Spectroscopic Instrument (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.

SUMMARY AND FUTURE PROSPECTS
In this work we have used the second data release of the LO-FAR 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 <  < 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  BH is higher for RL QSOs.Furthermore, using a sample matched in ,   and  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  and   .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 radioloudness shows that all more radio-loud bins ( 2 ,  3 ,  4 ) differ from the radio-quietest ( 1 ), with features consistent with the classical RL and RQ division.
-These changes cannot be explained by the addition of radioundetected sources, as we have found that the radio-detected and radio-undetected quasars in  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 <  < 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.

Figure 1 .
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.

Figure 2 .
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  as denoted.

Figure 3 .
Figure 3.A corner plot summary of the redshift, absolute -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 ,   and  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 .
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.

Figure 5 .
Figure5.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.

Figure 6 .Figure 7 .
Figure 6.A comparison between various composite spectra as indicated by the legend, where the thick lines indicate the flux density, whereas the shaded area corresponds to ± the uncertainty.The upper panel presents the full wavelength coverage, where each composite is normalised to the RQ composite at 3000 − 3600Å.The lower left panel presents a zoom-in window indicating the differences bluewards of Ly, whereas the lower right panel -the discrepancy above 5000Å which is likely caused by host galaxy contamination.

Figure 8 .
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

Figure 9 .
Figure 9.The  −   plane for the different radio-loudness bins.The lowest radio-loudness bin ( 1 ) is shown in light blue, whereas the bins with increasing R are presented in dark blue throughout the three panels.

Figure 10 .
Figure 10.The emission line flux ratio of a given pair, as indicated in the legend, as a function of redshift.Each panel represent an emission line under comparison as denoted in the upper right corner, which must be significantly detected (>3) for inclusion.The dashed line in each panel represents the line of equality, facilitating visual comparison of the emission line flux ratios between the different pairs.

Figure 11 .
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 12 .
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 Table6ofHardcastle et al. (2009), where the solid line indicates the regression line, whereas the dashed lines show ±3× the obtained scatter.

Figure A1 .
Figure A1.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.
This paper has been typeset from a T E X/L A T E X file prepared by the author.

Figure A2 .
Figure A2.A Monte Carlo simulation representing the comparison between RL and RQ quasars under the null hypothesis that the RL population is drawn at random from the RQ population.The upper panel of each bin presents the composite spectra of the RL test sample (dark blue) and the RQ test sample (light blue), whereas the lower panel indicates the residual in units of propagated uncertainty (grey).The -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 A3 .Figure A4 .
Figure A3.A high S/N comparison between the RL and RQ population matched in ,   and  BH .As in Figure7, 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 figureA2, -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 A5 .
Figure A5.An example of the spectral fitting procedure with PyQSOFit.The top panel presents the composite spectra of the RQ population in the redshift range of 0.6 <  < 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.

Table 1 .
The -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 ( 2 ,  3 and  4 ) and the radio-quietest sample ( 1 ), and the radio-detected ( 1D ) and the radio-undetected ( 1U ) sample for the radio-quietest bin defined in section 4.3 per a given redshift bin.

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.