The radio loudness of SDSS quasars from the LOFAR Two-metre Sky Survey: ubiquitous jet activity and constraints on star formation

We examine the distribution of radio emission from ~42,000 quasars from the Sloan Digital Sky Survey, as measured in the LOFAR Two-Metre Sky Survey (LoTSS). We present a model of the radio luminosity distribution of the quasars that assumes that every quasar displays a superposition of two sources of radio emission: active galactic nuclei (jets) and star-formation. Our two-component model provides an excellent match to the observed radio flux density distributions across a wide range of redshifts and quasar optical luminosities; this suggests that the jet-launching mechanism operates in all quasars but with different powering efficiency. The wide distribution of jet powers allows for a smooth transition between the 'radio-quiet' and 'radio-loud' quasar regimes, without need for any explicit bimodality. The best-fit model parameters indicate that the star-formation rate of quasar host galaxies correlates strongly with quasar luminosity and also increases with redshift at least out to z~2. For a model where star-formation rate scales as $SFR \propto L_{bol}^\alpha (1+z)^\beta$, we find $\alpha = 0.47 \pm 0.01$ and $\beta = 1.61 \pm 0.05$, in agreement with far-infrared studies. Quasars contribute ~0.15 per cent of the cosmic star-formation rate density at z=0.5, rising to 0.4 per cent by z=2. The typical radio jet power is seen to increase with both increasing optical luminosity and black hole mass independently, but does not vary with redshift, suggesting intrinsic properties govern the production of the radio jets. We discuss the implications of these results for the triggering of quasar activity and the launching of jets.


INTRODUCTION
The fundamental physical mechanism that powers and defines active galactic nuclei (AGN) is the transfer of energy from the relativistically-deep potential well of a central supermassive black hole (SMBH) in a galaxy (Salpeter 1964;Zel'dovich 1964;Lynden-Bell 1969). The most luminous members of the AGN class are quasars, which can outshine their host galaxy on the order of hundreds or even thousands of times. Traditionally, the ratio of the radio flux density to the optical flux density of quasars has been used in the ★ E-mail: pnb@roe.ac.uk literature to divide quasars into two categories: radio-loud (RL) and radio-quiet (RQ; e.g. Kellermann et al. 1989). By this definition, of the order of 10% of quasars are observed to be RL, although this fraction is known to increase with both optical (Padovani 1993) and X-ray (della Ceca et al. 1994) luminosity. Nevertheless, even RQ quasars have been shown to emit weak radio emission when sufficiently deep radio data are available (e.g. Padovani 2016, and references therein).
Despite quasars being identified over half a century ago (Schmidt 1963), the question of whether RL and RQ quasars are two physically distinct populations has still not been answered. A dichotomy appears to exist in the sense that observations show an asymmetric distribution of flux ratios, with a long tail towards high radio lumi-nosities where only a small fraction of quasars lie. However, firm proof on whether this is due to a bimodal distribution or simply an asymmetric continuous radio luminosity distribution, is still elusive, as attempts to hunt for a bimodality in the radio loudness distribution of quasars have provided contradictory conclusions. Some authors have argued that the bimodality exists (e.g. Ivezić et al. 2002;White et al. 2007) while others claim there is no bimodality (e.g. Cirasuolo et al. 2003a;Cirasuolo et al. 2003b;Baloković et al. 2012). An inherent problem is the selection biases that exist in these studies due to the definition of radio-loudness using optical and radio information, and the use of flux-limited samples at both radio and optical wavelengths. Furthermore, the use of different observing bands can give varying results (Ivezić et al. 2002). Lacking the answer to such a fundamental question about the nature of RL and RQ quasars means we cannot obtain a complete understanding of the physical mechanisms at play. Moreover, this gap in our knowledge has significant consequences for theories of galaxy formation and evolution given the strong evidence for a close relationship between the growth of the central SMBH and the evolution of its surrounding host galaxy (e.g. see reviews by Fabian 2012;Heckman & Best 2014), such as the correlations observed between the mass of the central SMBH and properties of the host galaxy's bulge (e.g. Ferrarese & Merritt 2000;Gebhardt et al. 2000;Marconi & Hunt 2003;Häring & Rix 2004).
The source of the radio emission from 'RQ' quasars is also a debated issue (see Panessa et al. 2019, for a review). Star-formation (SF) produces free-free emission from HII regions as well as synchrotron radiation from electrons accelerated to relativistic speeds in supernova remnants (Condon 1992). As the hosts of quasars are often star-forming galaxies (see review by Heckman & Best 2014), the question then becomes whether SF in the host galaxy is sufficient to account for the observed radio emission from RQ quasars. Some studies have found that the radio emission from some RQ quasars could be explained solely by SF (e.g. Kimball et al. 2012;Condon et al. 2013). Others suggest that SF is not sufficient and, hence, that the majority of the radio emission from RQ quasars must come from the AGN (e.g. Zakamska et al. 2016;White et al. 2015White et al. , 2017, either in the form of small-scale jets (Cirasuolo et al. 2003a), or from other processes that produce weak radio emission such as AGNdriven winds or disk coronal activity (e.g. Laor & Behar 2008). In many cases, high resolution radio maps of RQ quasars have detected non-thermal radio core emission or emission coming from extended jet-like structures (e.g. Kukula et al. 1998;Blundell & Beasley 1998;Leipski et al. 2006;Klöckner et al. 2009;Maini et al. 2016;Herrera Ruiz et al. 2016;Jarvis et al. 2019;Hartley et al. 2019). It could be the case therefore that jets are a feature in all quasars but are often unresolved. This would not be surprising since such low luminosity radio jets are commonly seen in massive galaxies (although these radio-AGN are not optically classified as quasars; e.g. Heckman & Best 2014;Mingo et al. 2019;Baldi et al. 2021), with Sabater et al. (2019) finding evidence that locally they are essentially always present in the most massive galaxies. Mancuso et al. (2017) attempted to explain the abundances of SF-dominated and jet-dominated RQ quasars using a model of in-situ evolution whereby the origin of the dominant radio emission changes as the black hole grows.
Many studies have investigated the SF component of quasars, and in particular how the star-formation rate (SFR) relates to the optical or X-ray luminosity of the quasar, and how it evolves with redshift. The relation between SFR and quasar luminosity is particularly interesting because it effectively relates the growth rate of the SMBH to that of the galaxy around it, and hence the build-up of the black hole mass versus bulge mass relation. Netzer (2009) and Bonfield et al. (2011) both suggested a strong correlation between AGN lumi-nosity and star-formation rate using AGN from the Sloan Digital Sky Survey (SDSS; York et al. 2000). More recent studies of (often lower luminosity) quasars at higher redshifts find that the SFR increases strongly with increasing redshift out to at least ∼ 2, but that any trend of SFR with quasar luminosity is either weak (e.g. Harrison et al. 2012;Azadi et al. 2015;Stanley et al. 2017;Stemo et al. 2019) or insignificant (e.g. Rosario et al. 2012;Mullaney et al. 2012;Stanley et al. 2015) once the redshift effects are accounted for. However, the quantitative details remain widely debated. Furthermore, there is an exception to this result at the highest quasar luminosities, where most studies find that SFR does correlate with AGN luminosity (e.g. Shao et al. 2010;Rosario et al. 2012;Harris et al. 2016;Lanzuisi et al. 2017); this has been attributed to the triggering of these quasars by major mergers of galaxies.
The physical mechanisms that produce radio jets in RL (and at least some RQ) quasars are also still uncertain. Some proposed theoretical models are compelling but lack observational confirmation. Blandford & Znajek (1977) proposed that a spinning BH, threaded by magnetic field lines, can produce anti-parallel jets of energy, while Blandford & Payne (1982) proposed that outflows of matter and energy from a rotating disk of gas accreting on to a spinning black hole can form jets when certain conditions regarding the orientation of the magnetic field with respect to the disk are met. BHs can get spun-up by recent galaxy mergers and perhaps by accretion events (e.g. Dotti et al. 2013), andRawlings (2011) found that using this assumption they were able to explain the radio luminosity functions of both high-and low-excitation radio galaxies. However, no conclusive evidence has yet been found that suggests the BH spin is involved in generating radio jets in quasars, due to the extreme difficulty in measuring the spin of the BH (e.g. see the contradictory results for galactic-scale black holes from Steiner et al. 2013;Russell et al. 2013). Wilson & Colbert (1995) argued that BH spin must provide the fundamental distinction between RL and RQ quasars given that they found no dependence of radio loudness on other physical parameters, such as the mass of the BH and the accretion rate. However, while some later studies back up these findings (e.g. Woo & Urry 2002), many others do find strong dependencies of radio loudness on the BH mass (e.g. Laor 2000;Lacy et al. 2001;Dunlop et al. 2003;McLure & Dunlop 2004;Best et al. 2005). RL quasars are also found to be more highly clustered than RQ quasars (Retana-Montenegro & Röttgering 2017), suggesting that larger-scale environment may play a role. van Velzen & Falcke (2013) find a very tight relationship between jet and disk luminosities in RL quasars, and infer that if BH spin is a major factor in the jet power then the RL quasars must all have very similar spin values, with then a wide gap in spin to the RQ population; this would require a strong dichotomy in the quasar population. Sikora & Begelman (2013) argue that variations in the strength of the magnetic flux threading a spinning black hole may instead be the primary factor that controls the strength of radio jets, while several recent works (Klindt et al. 2019;Rosario et al. 2020;Fawcett et al. 2020) have found that red quasars have a factor ∼ 3 higher radio-detection fraction than blue quasars (see also Richards et al. 2003;White et al. 2007;Calistro Rivera et al. 2021), and argue that an evolutionary sequence may be occurring. Similarly, Morabito et al. (2019) found a higher radio-loud fraction in broad absorption line quasars and evidence for a link between the radio activity and an outflow phase.
Hence, there are several prominent questions still at large: (i) Are all quasars part of the same population?
(ii) Is SF sufficient to account for the radio emission observed from RQ quasars or are small-scale radio jets prevalent?
(iii) How does the SF rate depend on quasar properties (redshift, luminosity) and what does this tell us about the triggering of the quasar activity?
(iv) What are the physical mechanisms that influence the prevalence or strength of radio jets in quasars?
The LOw Frequency ARray (LOFAR; van Haarlem et al. 2013) Two-metre Sky Survey (LoTSS; Shimwell et al. 2017) is an ongoing radio survey of the northern sky in the frequency range 120-168 MHz, detecting an order of magnitude higher sky density of sources than any previous large-area radio survey (see Section 2.1 for more details). LoTSS is providing deep radio imaging of large samples of quasars. Furthermore, in the low-frequency radio regime, any extended radio structures (which are likely to have steeper radio spectral indices) are more prominent than at higher frequency, and therefore potential biasing effects due to Doppler boosting are much less pronounced than at GHz frequencies.
Recently, Gürkan et al. (2019) used the LoTSS data release 1 (DR1; Shimwell et al. 2019;Williams et al. 2019;Duncan et al. 2019) data over the Hobby-Eberly Telescope Dark Energy Experiment (HETDEX; Hill et al. 2008) Spring Field region and the LO-FAR Herschel-ATLAS North Galactic Pole survey (H-ATLAS NGP; Hardcastle et al. 2016) to examine the low-frequency radio properties of optically selected quasars from the SDSS Baryon Oscillation Spectroscopic Survey (BOSS; Dawson et al. 2013). The wide area and high sensitivity of the LoTSS DR1 allowed Gürkan et al. (2019) to determine radio luminosities, or place meaningful upper limits on these, for tens of thousands of quasars. They investigated how the radio loudness of these quasars depended on other galaxy and BH parameters such as BH mass, optical luminosity, radio luminosity, redshift, and Eddington ratio. Given their results, Gürkan et al. favour the scenario where AGN jets (of a wide range of powers) and star formation-related processes both contribute to the radio emission observed from quasars and that there is no RL/RQ dichotomy, but rather a smooth transition between the regimes where each of the two processes dominate.
We aim to build upon the results of Gürkan et al. by taking the simple approach of constructing and testing a numerical, two-component model of the radio flux densities of quasars. The model is a superposition of the two expected sources of radio emission from galaxies, the AGN (jets) and the SF, each modelled from physical prescriptions. The model implicitly assumes that no intrinsic bimodality exists but rather there is a smooth transition from a star-formation dominated to a jet-dominated regime as the radio jet power increases. Such an approach allows us to generate simulated samples through Monte Carlo realisations that can be compared to observed data. Quantifying the validity of the model and constraining its parameters will provide information relevant to questions (i) and (ii) above, while investigating how the model parameters change as a function of properties of the quasar such as redshift, optical luminosity (for which we use the absolute i-band magnitude as a proxy) and BH mass provides input into questions (iii) and (iv).
The outline of the paper is as follows. Section 2 outlines the data used during this research to build a large sample of quasars used to validate and constrain the model. The two-component model of the radio luminosity distribution of quasars is detailed in Section 3. Section 4 presents the results found, and the physical interpretations of these results are discussed in Section 5. Finally, a summary of our conclusions is given in Section 6. Cosmological parameters are taken to be (Ω , Ω Λ ) = (0.3, 0.7) and 0 = 70 km s −1 Mpc −1 .
LoTSS DR1 includes the radio images of the relevant region as well as a source catalogue (Shimwell et al. 2019), where the Python Blob Detector and Source Finder (PyBDSF; Mohan & Rafferty 2015) algorithm has been used to catalogue sources. However, the PyBDSF catalogue produced will contain a number of sources that were associated incorrectly. The reasons for this could be: the blending of physically distinct sources into a single catalogue entry; the separation of the components of extended sources into different PyBDSF catalogue entries; spurious emission or artefacts. To account for this, Williams et al. (2019) present a value-added catalogue in DR1 where significant effort (both statistical techniques and extensive visual analysis known as LOFAR Galaxy Zoo) has gone into ensuring, as much as possible, that the catalogue is a true representation of the radio sources in the relevant region. Furthermore, Williams et al. provided optical/infrared counterpart identifications (where detected) for all of the LoTSS sources, making use of optical and infrared data from the Panoramic Survey Telescope and Rapid Response System (Pan- STARRS Chambers et al. 2016) 3 survey and the Wide-field Infrared Survey Explorer (WISE; Wright et al. 2010). The detailed processes implemented to produce such a catalogue can be found in Williams et al. (2019).

Sloan Digital Sky Survey quasar sample
Optical data for quasars from the fourteenth data release (DR14Q; Myers et al. 2015;Pâris et al. 2018) of the SDSS were obtained, using the catalogue detailed in Pâris et al. (2018). The catalogue of over half a million quasars is described as a "superset" that is a compilation of all spectroscopically confirmed quasars from SDSS-I, II, III and IV (see also Richards et al. 2002;Ross et al. 2012, for earlier selection criteria). It was also necessary however, to obtain the SDSS-DR7 quasar catalogue (DR7Q) that consists of all spectroscopically confirmed quasars from SDSS-I/II (Schneider et al. 2010) in order to obtain the target selection flags of the 79,487 quasars that were not re-observed as part of SDSS-IV. The target selection flags allowed us to identify and remove quasars that were selected for observation solely based on their radio emission, in order to mitigate any selection bias (see Section 2.3).
Absolute i-band magnitudes were also obtained from the DR14Q catalogue. We convert the magnitude given in the catalogue (a magnitude K-corrected to = 2), to a magnitude K-corrected to = 0, ( = 0), assuming an optical spectral index of 0.5 (Richards et al. 2006). Galactic extinction corrections for the i-band were also applied. The extinction correction for each quasar, obtained from the Schlafly & Finkbeiner (2011) dust maps, was given in the DR14Q catalogue.

Our Sample
Our aim was to build a sample of quasars observed in the optical by the SDSS, which provides information on their redshift and optical luminosity, and to combine this with the radio properties of the quasars, such as their integrated flux density, from LoTSS. We started with the ∼ 51, 000 quasars within the LoTSS DR1 RA and Dec limits from the SDSS DR14Q catalogue described in Section 2.2. The following cuts were then applied in order to ensure a robust and unbiased analysis: (i) Sources brighter than = −40 were removed. These were likely misidentified during the automated process that generated the DR14Q catalogue.
(ii) The sample was limited out to a redshift of = 6, as beyond that the contamination from unreliable redshifts is high (in practice our analysis was restricted to even lower redshifts by sample size limitations).
(iii) Sources that lay outside the LoTSS coverage, or fell within gaps of the LoTSS mosaics, were removed, due to the lack of radio data.
(iv) Finally, a further 199 sources that had target selection flags in SDSS DR7Q or DR14Q that indicated they were included in the spectroscopic sample solely because of their radio emission were removed, to mitigate any possible bias (otherwise, in regions outside of the SDSS colour-space selections, radio-bright quasars would be preferentially included, potentially biasing the results).
The final sample had a size of 42,601 quasars. The distribution of these quasars in − space is shown in Fig. 1. As expected, the sample probes systematically higher optical luminosities at higher redshifts, but still has a good dynamic range in optical luminosity at each redshift.

Radio Flux Densities
To determine radio flux densities for our quasar sample, we first crossreferenced the coordinates of the two survey catalogues, LoTSS and SDSS. It should be noted here that the coordinates of the optical identification given in the LoTSS value-added catalogue (i.e. that of the cross-matched Pan-STARRS or WISE host galaxy) were used rather than the less accurate radio-derived coordinates from PyBDSF; the Pan-STARRS coordinates are aligned to the SDSS co-ordinates to typically much better than an arcsecond. Fig. 2 shows the number of quasars matched to a radio source (their nearest) as a function of the maximum cross-matching angular distance. To estimate the contamination one would observe for a given matching radius, we select random locations within the LoTSS DR1 coverage and measure the number of random matches as a function of radius. Fig. 2 also shows the number of matches out to each cross-matching radius after correcting for this random contamination. From this, the maximum cross-matching radius was chosen to be 1.5 . This resulted in just under 5750 direct matches with a predicted contamination of ∼ 0.1%. For these cross-matched sources, we extract integrated flux densities (and associated uncertainties) directly from the Williams et al. (2019) value-added catalogues. These integrated flux densities correctly incorporate any extended radio structures due to the source association process.
For the remaining optical quasars that were undetected in the LoTSS catalogue, a flux density was extracted directly from the LoTSS mosaics, using the image value in the mosaics at the coordinates of the quasars given in DR14Q. The uncertainty on this was extracted at the same location from the LoTSS rms noise maps. Although these sources are individually below the 5 catalogue S/N limit and so have measured flux densities dominated by the noise (with many having negative values), it is nonetheless expected that the genuine low significance emission will lead to considerable information being present in the exact flux density distribution (cf. Roseboom & Best 2014;Malefahlo et al. 2020, see also Fig. 3). An underlying assumption of extracting flux densities directly from the radio maps is that these faint sources must be compact compared to the 6 LoTSS beam (such that the peak flux density traces the integrated flux density); this is a reasonable approximation as any radio emission from these faint sources will primarily be star-formation on galaxy scales, radio cores and/or small-scale jets. It should be noted that this faint emission will not have been properly cleaned in the radio imaging step, which may lead to a slight under-estimate of the true flux densities.
For a small number (a few tens) of quasars, despite the lack of a LoTSS catalogue match, the radio flux density extracted from the mosaics had a flux density with greater than 5 significance. A selection of these sources were examined visually, and found to represent a mixture of cases such as sources PyBDSF had failed to detect, incorrect LoTSS IDs and quasars overlying the extended radio emission of other sources. In the former cases, the integrated flux density may be underestimated if the radio source is extended, while in the latter case it will be overestimated. As there were relatively few of these sources, and flux densities may be biased in either direction, no attempt was made to correct these errors.

Bolometric Luminosity and Accretion Rates
It is also useful at times during our analysis to relate the absolute iband magnitude 3 to a bolometric luminosity ( bol ) and to an estimate of the growth rate of the SMBH, BH . To do this, we first relate the absolute i-band magnitude to the absolute b-band magnitude using the relation determined empirically from a sample of 1046 quasars by Richards et al. (2006): − ( = 2) = 0.66 ± 0.31. This allows us to estimate the bolometric luminosity applying the empiricallyderived relation given in McLure  (1) Although there may be slight discrepancies between the b-band filters used to calibrate the above relations, the relations should give a sufficiently good estimate. The bolometric luminosities of the quasars in the sample range from about 10 38 to 10 40 W. The BH growth rate is then given by where is the efficiency at which the rest-mass of the accreting material is converted to energy; we assume a typical value of = 0.1.

Virial Black Hole Masses
To further our analysis of the radio loudness distribution of the sample beyond dependencies on optical luminosity and redshift, we have also obtained virial BH mass estimates. These have been obtained from catalogues provided by Shen et al. (2011) and Kozłowski (2017). These catalogues only include quasars up to the SDSS data release 12, and so do not include all of the quasars in our sample, but should represent a relatively unbiased subsample. The methods used to obtain the black hole mass estimates are briefly discussed here. Shen et al. (2011) present properties of quasars in the SDSS DR7Q catalogue based on spectral fits. We use the measurements of the virial BH mass based on the broad MgII and CIV emission lines, which have been calibrated as virial BH mass estimators. Several estimates are given for varying calibrations of the MgII line by Shen et al. (2011). Kozłowski (2017) make use of the broadband, extinction-corrected magnitudes from quasars in the SDSS DR12Q catalogue to derive monochromatic luminosities that can then be combined with the broad emission line widths of MgII and CIV to estimate BH masses.
The Shen et al. (2011) and Kozłowski (2017) catalogues were matched to the SDSS-LoTSS sample described in Section 2.3, using positional cross-match (all matches being found within ∼ 0.6 , with essentially no contamination). The two sets of black hole mass estimates are broadly comparable, and the precise choice of which to use has no qualitative effect on the results that we obtain. Where available, we use the estimates provided by Shen et al. (2011) and the weightedmean of the available estimates was taken. Otherwise, the estimates provided by Kozłowski (2017) were used, where the MgII estimate was prioritised over the CIV estimate because Kozłowski (2017) uncovers a bias with measurements from the CIV line. Of the 42,601 sources in the SDSS-LoTSS sample, 24,096 have estimates for the virial BH mass. The derived black hole masses typically range from 10 8.5 to 10 9.5 M (although at lower redshifts the lower-luminosity quasars have lower black hole masses, down to 10 7.5 M ).

A TWO-COMPONENT MODEL FOR RADIO EMISSION
We started from the simple approach of building a two-component model of the radio luminosity distribution of quasars. The model assumes that two sources (AGN and SF) contribute to the observed radio emission of every quasar. Thus, in our model, all quasars are assumed to display star-formation activity at some level (as would be expected, given the large gas content that is present in order to fuel the black hole accretion, and in line with many studies of quasar host galaxies; e.g. Shao et al. 2010;Floyd et al. 2013;Harris et al. 2016). In Section 3.1 we describe how we assign an SFR to each quasar, drawing from an inferred Gaussian distribution. In addition, in our model all quasars are assumed to possess radio jets. As we motivate in Section 3.2, the jet luminosity is allowed to vary in strength from the very powerful radio jets seen in the most radio-loud quasars down to the very weak small-scale radio jets that have been observed in high angular resolution, sensitive radio images of some radio-quiet quasars. The large range of possible jet powers is the primary factor that sets the overall radio luminosity of the system, and determines whether the AGN or the SF is the dominant source of radio emission.
We create simulated samples of quasars from our model through Monte Carlo realisations, summing the radio luminosity contributions from the SF and AGN components, converting these to a radio flux density, and adding noise (Section 3.3). We then compare these with the observed distribution of quasar flux densities in order to determine the validity of our model prescription for the sources of radio emission and to constrain the model parameters (as outlined in Section 3.4).
As it is expected that the strength of both the star formation and the jet component may be dependent on the optical luminosity of the quasar, and may vary with redshift, we carry out this comparison on subsamples of quasars produced by separating the main sample in − space by the grid lines shown in Fig. 1. We separately analyse each subsample in each grid square that hosts over 500 quasars. Analysing the distribution in the 2-dimensional opt − space is preferred to marginalising the distributions, as Jiang et al. (2007) showed that the strong correlation between opt and can lead to marginalised studies obtaining inaccurate results. Within each such grid square, the optical luminosity and redshift are reasonably constant (Δ = 1; Δ = 0.4), so fitting the same model parameters for the star formation and jet components for all quasars within each bin should produce robust results. By analysing the results from the different individual grid squares separately, we are able to recover detailed information as to how our model parameters for the SF and jet components vary as a function of cosmic time and optical luminosity.

SF Component
The radio luminosity function of star-forming galaxies is often modelled as a broken power-law, where the wide distribution of starformation rates arises from the large range of stellar masses of the star-forming galaxies, combined with the tight relation (scatter ∼0.2-0.35 dex) that star-forming galaxies show between their SFRs and their stellar masses (often called the star-forming main sequence; e.g. Noeske et al. 2007;Elbaz et al. 2007). Low star formation rates also arise from the quiescent galaxy population.
The host galaxies of powerful AGN are typically both massive (e.g. McLure et al. 1999;Best et al. 2005) and star-forming (e.g. Kauffmann et al. 2003). Investigations have found that the host galaxies of radiatively-efficient (quasar-like) AGN mostly lie on or above the SFR-mass relation (e.g. Mainieri et al. 2011;Heckman & Best 2014), avoiding the quiescent galaxy population 4 . Since they are all highmass star-forming galaxies, the SFRs of quasar host galaxies would therefore be expected to have a much narrower distribution than a full power law. Based on this, we model the radio emission of the SF component from each quasar (in a given bin of redshift and optical luminosity) as being drawn from a Gaussian distribution (in log space) with two free parameters: the mean, log( /[W Hz −1 ]) at the frequency of LOFAR (∼150 MHz), and the standard deviation, /dex. The normalisation of the star-forming main sequence is known to evolve strongly with redshift due to the higher availability of gas in the early Universe (e.g. Speagle et al. 2014): the mean star formation rates of massive galaxies increases by a factor ∼ 30 from = 0 to ∼ 2, in line with the evolution of the cosmic star formation rate density (e.g. see review by Madau & Dickinson 2014). Furthermore, as discussed earlier, many studies have investigated how the typical SFR of quasars varies with their AGN luminosity, finding different results. For this reason, the two Gaussian parameters (mean and standard deviation) are allowed to take different values in different bins of redshift and luminosity. In each Monte Carlo realisation, a random luminosity for the SF component of each quasar ( SF ) is drawn from this Gaussian distribution. The determination of the best-fitting values for and in the different redshift and optical luminosity bins then reveals how the SF component evolves across cosmic time and how it connects to the black hole accretion rate.
To yield physical information about the system, it is useful to relate the mean luminosity of the SF component to a SFR, Ψ. which agree to within 0.1 dex at the typical luminosities of the quasars in our sample ( 150 ∼ 10 23 − 10 24 W Hz −1 ). Here we use the latter relation, which has been calibrated using LOFAR data. To relate the width of the Gaussian, which is in log space, to a width in log Ψ space, Eqn. 4 gives Ψ = 0.93 .

AGN Component
The radio luminosity function of 'radio-loud' quasars (i.e. those with a dominant AGN component) is often modelled as a broken powerlaw (e.g. Dunlop & Peacock 1990;Kaiser & Best 2007;Kimball et al. 2012;Best et al. 2014). However, the radio luminosity function maps out the full source population, and it is likely that quasars of different optical luminosity may dominate different parts of this distribution. Furthermore, a broken power-law model requires two additional free parameters over a single power-law formalism, and the data do not justify these additional parameters: when attempting to model the AGN (jet) component in given redshift and radio luminosity bins as a broken power law, we found that we could not constrain the slope above the break luminosity and there were strong degeneracies between other parameters. Instead, a single power-law formalism was found to be sufficient to describe the AGN component of the model. To model the AGN component, we therefore draw a luminosity randomly from a single power-law distribution, with probability distribution function (PDF) where 0 is the normalisation and is the slope. The PDF was defined such that has units (Δ ) −1 . Defining as such means that the integral of the PDF must adhere to: where is the maximum radio luminosity obtained by any radio quasar, and is the minimum jet luminosity of the quasars. The choice of is not critical, as the integral of the model in the range → ∞ will be negligible provided that is sufficiently large. In practice, we set log( /[W Hz −1 ]) = 30, above which we do not expect to see any sources. The lower limit, , is then fixed, for a given normalisation and slope, by Equation 6.
To set the normalisation of the power-law, 0 , in to more physically intuitive units, we define the quantity as the fraction of the integral of the PDF at luminosities brighter than , where was set as log( /[W Hz −1 ]) = 26 (choosing a different value just leads to a scaling of , but does not affect the trends seen). The value of 5 Note that there is an error in the conversion from to SFR for a Chabrier IMF in footnote (b) of Table 3 of Brown et al. (2017), where the Salpeter to Chabrier IMF conversion factor appears to have been inversely applied; the values provided here use the correct conversion. was chosen to be high enough that jet emission will dominate and SF will be negligible, and thus variation of maps directly onto the variation in the fraction of high-power jet-dominated radio sources. The two free parameters in the model are now and .
It can be shown that, Asserting equation (6) now means can be computed:  (2019) investigate AGN fractions as a function of stellar mass, and find that the cumulative AGN fractions reach 100 per cent (and hence the luminosity function must turn over) at around 150MHz ∼ 10 22 W/Hz at the highest stellar masses, with the cut-off luminosity decreasing with decreasing stellar mass. The range of lower limits determined for the model, ∼ 10 20 − 10 22 W/Hz, agrees well with these observations, giving confidence that the model is producing sensible results.
To randomly sample from the single power-law, the inverse cumulative method was implemented. Briefly, this method involves building the cumulative distribution function (CDF): where is a uniformly distributed variate on (0,1), and then inverting the function to solve for . The luminosity variates of the AGN component can therefore be sampled from:

Total Simulated Flux Density
For each quasar, for a given set of values of the four free parameters in the model ( , , and ) a luminosity is randomly drawn (independently) for each of the components (i.e. a SF luminosity and an AGN luminosity). The two luminosities are then added together to give a total luminosity: Using the known redshift of the quasar, a simulated flux density can then be computed. Remembering that we need to apply a Kcorrection, assuming ∝ − , the flux density can be computed using: Figure 3. The observed flux density distribution for quasars with 2.0 < ≤ 2.4, −24 < ≤ −23 (solid histogram) and the best-fit simulated flux density distribution for these (dashed line). The simulated distribution is generated by the model using the set of model parameters that provide the best fit to the observed flux density distribution for this subsample (see Section 3.4). This subsample hosts 1127 quasars and for each quasar, 100,000 representations of the model have been returned (and averaged). In order to clearly display both the low flux densities (with some noise-dominated negative values) and the tail to very high flux densities, the distribution is plotted in two parts. The left panel shows the observed flux density distribution at low flux densities, using a linear flux scale and equally-sized bins; the vertical dotted line indicates the median flux density, showing that this is offset from zero. The right panel shows the distribution at higher flux densities on a logarithmic flux scale, binned using the Bayesian Blocks (Scargle et al. 2013) formalism (note that the full unbinned distribution is used in the KS test). For the right panel, the count, , has been divided by the width of the respective bin (therefore /mJy) since the bins have varying width. Over all flux densities, the model provides an excellent match to the observed distribution.
where 0 is the frequency of LOFAR observations, = 0.7 is the radio spectral index (e.g. Calistro Rivera et al. 2017) and is the luminosity distance. Gaussian noise is then added to simulate the observation, where the width of the Gaussian is taken to be the rms extracted from the LoTSS mosaic at the position of that quasar. This process is repeated for all quasars to produce a single representation of the simulated flux density distribution for the quasar population. We average across typically 1,000 Monte Carlo representations to then allow a direct comparison between the distribution of simulated and observed flux densities to be conducted.
An example of the comparison between the observed flux density distribution for a given subsample (as defined in Fig. 1) and the bestfit simulated flux density distribution for the same subsample can be seen in Fig. 3. The procedure to obtain the best-fit model for a given observed flux densitiy distribution is detailed in Section 3.4. It is clear that the model is able to provide a good match to the observed flux density distribution. It is also notable that the peak of both the observed and simulated distributions are offset from 0 mJy (this offset effectively constrains the SF component of the model) and that even well below the 5-sigma noise limit of the radio data (around 0.35 mJy) the observed flux density distribution is clearly non-Gaussian (with the tail of the distribution to high flux densities constraining the jet contribution); this highlights the valuable information available within the noise of the radio data, and which is well-fitted by the model. It is worth re-emphasising that our model assumes that no intrinsic bimodality exists in the radio flux distribution of quasars, as illustrated in Fig. 3.

Obtaining best-fit parameters
The two-sample Kolmogorov-Smirnov (KS) test is a non-parametric process to determine the probability that two sets of data have been drawn from the same PDF. The test can therefore be used to compare the simulated samples generated through the Monte Carlo realisations with the observed radio flux density distributions. This comparison will allow us to quantify the validity of our model as well as to constrain the values of the free parameters. Observed radio flux density distributions are taken from the subsamples described in Section 2.3 and one example of these is illustrated in Fig. 3.
We found from initial tests that a simple KS test on the full flux density distributions was not sensitive to the AGN component of the model, and that the degeneracy between the normalisation, , and the slope, , was not being broken during our fitting. The reason for this is the relatively few sources in the high luminosity tail (the KS statistic is determined by the maximum absolute offset between the model and observed flux density distributions, and a tiny fractional difference in shape of the bulk of the population can lead to a bigger absolute offset in the cumulative distribution than that of a significant fractional offset in the tail of AGN). This same issue prohibited a flux-binned chi-squared approach. Instead, to solve this problem, we split the flux density distributions at a threshold of thresh = 1mJy. This allows us to compare the relative number of sources above and below the threshold using a chi-squared test (based on Poissonian statistics; this ensures a broadly correct number of objects in the high luminosity tail), while additionally performing the KS test separately on the flux density distributions either side of the threshold; this approach gives additional weight to the AGN component (which dominates the above-threshold distribution). The value thresh = 1mJy was chosen as it is the flux density at which typically the observed flux density distributions transition from the large peaked component to the tail component of the population. In Appendix A we demonstrate that (to well within the errors) our results do not depend on the choice of thresh .
To optimise the parameters of the model, we extract the probabilities, > and < , returned from the KS tests above and below the flux threshold, respectively. We then combine these with the chisquared value for the Poissonian number test, 2 P , to derive the overall likelihood (L) which is then minimised. Here, where < model (or > model ) is the number of sources that the model predicts to lie below (or above) the threshold and < obs (or > obs ) is the number of sources observed below (or above) the threshold, When finding the best-fit model for each of our subsamples, in order to avoid the risk of minimisation routines getting stuck in local minima (especially for less-populated, and hence noisier, regions of parameter space), and given that the problem was computationally manageable, we implemented a "brute force" method: we computed the − ln L value at every point on a multi-dimensional parameter grid to find the global minimum and the uncertainties on this. In each dimension of parameter space, 30 grid points were defined, equallyspaced in the 4-dimensional parameter space: log( /[W Hz −1 ]), /dex, , log( ). An example of the results of applying such a minimisation routine can be seen in Fig. 4.
To ensure that the final best-fit model provides a valid explanation of the data, we perform a single KS test between the overall modelled flux density distribution and that of the observed quasars in each bin in luminosity-redshift space. In all cases, the returned probability is ≥ 0.25, indicating that our two-component model for the radio flux density distribution of quasars is valid across a wide range of optical luminosities and redshifts.

Variation of model parameters with and
Given that the model seems to accurately reproduce the observed data, we can now look at how the best-fit parameters change with optical luminosity and redshift. The best-fit parameters along with their errors for the subsamples are visualised in − space in Fig. 5. For clarity, the collapsed version of these results is given in Fig. 6, and a full table of best-fit parameters is provided in Appendix B (Table B1). The parameters and that characterise the Gaussian component of the model are expressed respectively as their SF equivalents, Ψ and Ψ , through application of the calibration given in Section 3.1.
One immediate result of note is the lack of discernible trends, with either redshift or optical luminosity, in the width of the SF component and the slope of the AGN component. These findings are perhaps expected. That the scatter in the relationship between SFR and quasar luminosity does not change much with quasar luminosity has been seen before in literature plots (e.g. Lanzuisi et al. 2017), and as this ratio is driven by gas distributions within the host galaxies it is not surprising that there is no strong redshift dependence either. The slopes of the radio luminosity function are also known to not evolve significantly with redshift (e.g. Dunlop & Peacock 1990;Smolčić et al. 2017, and references therein). Given the lack of observed trends of these parameters, we took the decision to fix and and conduct further analysis of the remaining two parameters. Fixing and allows us to increase the resolution with which we fit the remaining parameters as well as removing the possibility of artificial jumps in the best-fit parameters for and caused by random variations in and . We therefore reserve discussion of the trends of the mean SFR rate and the jet power normalisation to Section 4.2.

Fixing &
We fix the width of the Gaussian and the slope of the power-law with the values = 0.45 (corresponding to Ψ = 0.42 from Equation 4) and = 1.4; these are the (suitably-rounded) weighted-mean values from the different grid cell fits. We also increase the resolution of our parameter grid such that the number of grid points in each of the two remaining dimensions is now 40.
Best-fit parameter values are provided in Table B1. The collapsed distributions of the best-fit parameters with and are shown in Fig. 7, while the full − ln L parameter space is visualised in Fig. 8. In Fig. 7 we also relate to an estimate of the growth rate of the SMBH, BH , as described in Section 2.4.2, as it provides some interesting physical information; this is shown as an upper x-axis. Reassuringly, the results with fixed and are in agreement with the trends observed in Fig. 6 when these two parameters were not fixed, but now with lower noise.
The SFR, Ψ, of the host is seen to increase with increasing optical luminosity (or SMBH growth rate). We also find that the SFR of the host increases with increasing redshift out to ∼ 2 − 3 at which point, we observe the SFR beginning to turnover. In Appendix C we confirm that this redshift trend is genuine and is not caused by selection effects related to correlations between and within each bin.
We also find a strong dependence of the jet power normalisation with optical luminosity, in line with previous studies (e.g. Jiang et al. 2007): at the lowest optical luminosities, we observe values for the fraction of sources at high radio luminosities ( 150MHz > 10 26 W Hz −1 ), ≤ 0.5% while at the highest optical luminosities, we reach values of ∼ 6% (see Figures 7 and 8). No discernible trend of with redshift exists, suggesting that the principal mechanism that drives the distribution of jet powers must be an internal property of the system.

Black Hole Mass Dependence
We investigate the effect of the BH mass on the model parameters, for the case of fixed = 0.45, = 1.4, using only those sources that have a BH mass estimate, as described in Section 2.4.3. To investigate the effect of BH mass on the model parameters, we split each subsample in − space in two at the median BH mass of the given subsample, BH,med . The median BH mass of each subsample was chosen rather than a uniform value for the entire sample as the BH mass varies across − space, as can be seen from the first panel of Fig. 9. Across most of the parameter space, the difference in the median BH mass between the higher and lower BH mass bin is a factor of 2.5-3.5.
The results of including the BH mass information can be seen in Fig. 10; see also Table B2. The general trends of the parameters of the model are still in agreement with those found and discussed in Section 4.2 and so will not be repeated. Instead, we focus on how the inclusion of the virial BH mass estimate affects the trends. We find little variation of the SFR trend with optical luminosity between the high and low BH mass bins. Similarly, at < 2, the SFR seems to be independent of the BH mass. However, at redshifts above ≈ 2 (where the SFR flattens) in almost all luminosity-redshift bins the SFR in the higher BH mass bin is higher than that at lower black hole masses, by an average of 0.17 dex. Although the significance of this is low (typically 1-2 ) in each case, the consistent direction of the offset across almost all of the high-redshift bins suggests a genuine effect, although more data would be required to confirm this.
Regarding the fraction of sources at high radio luminosities, for quasars with < 2 there are indications that an increase in BH mass boosts the value of : although once again the differences are of relatively low significance (typically 1-2 ) within each individual luminosity-redshift bin, for all 13 luminosity-redshift bins at < 2 the value of derived for the higher black hole mass subsample is greater than or equal to that of the lower black hole mass subsample (giving a much higher combined significance). The average difference is a factor of ∼ 2 in (see Table B2). This result suggests that, at < 2, the jet power normalisation may depend on both BH and (or BH ), independently. At higher redshift, > 2, there is no evidence for a significant black hole mass dependence of . The results indicate that the width of the SF Gaussian distribution ( ) and the slope of the jet power distribution ( ) are broadly constant with both optical luminosity and redshift. The SFR (Ψ) increases strongly with optical quasar luminosity, and more weakly with redshift, while the jet-power normalisation ( ) increases with optical luminosity but is independent of redshift.

DISCUSSION & PHYSICAL INTERPRETATIONS
Before we discuss the physical interpretations of our results, it is worth quickly summarising what we have found.
• The radio emission of quasars can be explained by our simple model, which assumes two sources of radio emission contribute in all quasars: SF in the host galaxy and the AGN. The SF in the host galaxy is modelled as being drawn randomly from a Gaussian distribution of given centre and width. The AGN (jet) luminosity is drawn randomly from a power-law distribution of given normalisation and slope. The four parameters of the model have been constrained and their dependencies on redshift, optical luminosity and BH mass have been investigated. The upper-left panel shows that the SFR increases with increasing optical luminosity (i.e. BH growth rate); the dashed lines indicate ratios of Ψ/ BH of 10, 50 and 250. The upperright panel shows that the SFR of the quasar hosts increases with redshift out to redshift ∼ 2 and then flattens. The lower panels show that the jet power normalisation correlates strongly with optical luminosity but shows no consistent trend with redshift.
• We find that the width of the Gaussian SF component and the slope of the power-law AGN component do not vary significantly with either redshift or optical luminosity. This is perhaps expected, as explained in Section 4.1, justifying our approach to fix the two parameters during further analysis. This allowed us to increase the resolution and clarity with which we constrain the remaining two parameters.
• We observe that the SFR of quasar host galaxies increases with redshift out to ∼ 2 and then flattens in the range ∼ 2 − 3. Independently, the SFR also increases strongly with the quasar luminosity. We observe little dependence of the SFR on BH mass until we get to the highest redshifts, at which point we see some evidence that increasing the BH mass may correspond to a small increase in the SFR.
• The normalisation of the jet power (and hence the fraction of sources at high radio luminosities) shows very little trend with redshift. However, it increases with increasing optical luminosity and there are indications that it may also increase independently with the SMBH mass.

The ubiquity of radio jets
We find that our two-component prescription of the radio emission, SF and AGN, is valid in describing the observed radio flux density distribution of quasars across a wide range of redshifts and optical luminosities. As our model inherently assumes a wide and continuous distribution of radio jet powers, we therefore argue that the historical dichotomy of RQ and RL quasars can be simply explained by whether the radio emission of the quasar is dominated primarily by the SF of the host galaxy or the powerful large-scale jets.
Given the validity of our model, this implies that both jets and SF are contributing in all quasars, as argued by Gürkan et al. (2019). Since quasars typically reside in star-forming galaxies, we do expect a contribution at some level from SF. However, the ubiquity of quasar jets is not necessarily expected. It should be emphasized that such ubiquity is not conclusively demonstrated by our analysis: when SF dominates the radio emission, our flux density distribution comparison cannot distinguish between the scenario where some quasars have a very weak jet contribution that lies significantly below the SF contribution (the lower limit of the power-law distribution from which jet luminosities are drawn is typically in the range log( /[W Hz −1 ]) ≈ 20 − 22, whereas the mean of the SF Gaussian distribution is typically log( /[W Hz −1 ]) ≈ 22.5 − 24.5) and the scenario where some quasars have no jet contribution to the radio emission (e.g. if there is a threshold for jet production). However, our results do demonstrate that radio emission from quasar jets is contributing at least down to luminosities where the radio emission from SF becomes comparable, which is well into the traditionally 'radioquiet' quasar regime. This is consistent with the high detection rates of emission from small-scale low-luminosity jet-like structures in deep, high angular resolution radio observations of some radio-quiet quasars (e.g. Herrera Ruiz et al. 2016;Jarvis et al. 2019;Hartley et al. 2019), and with the presence of radio jets with luminosities down to at least 10 21 W Hz −1 in massive galaxies that do not host quasars (e.g. Sabater et al. 2019;Mingo et al. 2019;Baldi et al. 2021).
In Fig. 11 we present the fraction of quasars for which the jet contribution to the radio emission is larger than that of the starformation contribution in our model, as a function of redshift and quasar luminosity. The derived values of 10-20% agree well with the widely-adopted 'radio-loud fraction' for quasars (e.g. Kellermann et al. 1989). We find a clear trend for an increasing jet-dominance with increasing optical luminosity, from of order 10% at optical luminosities > −24 to ∼ 20% at < −25. This is in line with most literature results for the dependence of radio-loud fraction on quasar luminosity (e.g. Jiang et al. 2007). Overall we find no uniform trend for how the jet-dominated fraction varies with redshift. At high luminosities ( < −25) there is indication for a decreasing jetdominated fraction with increasing redshift, as suggested by Jiang et al. (2007), but at lower luminosities any redshift trend is weaker, absent or non-monotonic. This lack of redshift evolution is in line with recent studies out to the highest redshifts (e.g. Liu et al. 2020, and references therein).
It is interesting to compare these models in detail with traditional techniques of selecting radio-loud AGN. The traditional radioloudness parameter, = 5GHz / 4400Å , can be calculated for these quasars using the LoTSS data and an assumed spectral index of = 0.7; in all of the redshift vs absolute magnitude space studied, the LoTSS data are sensitive enough to probe down to below the usual = 10 cut-off value. The models do not predict the jet and star-formation components for each individual quasar, but rather predict the distribution of properties for a population of quasars in each luminosity-redshift grid cell. To allow a direct comparison, we therefore consider a single iteration of the model in each grid cell, rank the modelled sources by total radio luminosity, and crossmatch these against the observed data in the same grid cell, similarly ranked by radio luminosity. Figure 12 compares the observed value of = 5GHz / 4400Å for each quasar against the modelled AGN / SF for its rank-matched model source. It is important to note that the ranking process adopted is only viable for those sources detected  by LoTSS (with signal-to-noise above 2), and that therefore sources with upper limits on their radio luminosities are not shown on the plot. These sources would fill out the lower left region of the plot; the detection limits (and typical AGN / SF ratios) are different in each grid cell in redshift and optical luminosity, and this produces the apparent shape of the cut-off towards the lower-left (which is artificial).
As can be seen from Figure 12, the traditional radio-loud definition of = 5GHz / 4400Å > 10 cleanly selects a population of quasars for which the jet is by far the dominant source of radio emission (typically AGN / SF > 10); the cut-off value of = 10 is shown to be well-motivated as cuts at lower values of would begin to pick up a population of quasars whose radio emission is dominated by starformation (quasars with > 1 and AGN SF are objects where the star-formation is strong enough to give significant radio emission; typically these are lower-redshift, lower optical lumoinosity quasars). Figure 10. The variation of the SFR (Ψ; top) and the jet power normalisation ( ; bottom) with optical quasar luminosity (left) and redshift (right), split into bins of higher (solid lines) and lower (dashed lines) black hole mass, as outlined in the text. We observe little difference in the SFR between the two BH mass bins, except perhaps at > 2 where we see hints of higher SFRs in the higher BH mass bin. For the jet power normalisation, , we find this is typically a factor ∼ 2 higher in the higher BH mass bin, at least for < 2. Figure 11. The fraction of quasars in our sample for which the contribution to the radio emission from a jet component exceeds that from star formation, as a function of redshift and quasar luminosity. The 10-20% values derived match typically-quoted values for the radio-loud fraction. We find that the radioloud fraction is typically higher for higher quasar luminosities but shows no uniform trend with redshift. Figure 12. A comparison between the traditional 'radio-loud' definition based on = 5GHz / 4400Å (with a cut-off value at = 10) and the modelled AGN / SF ratios, for all quasars detected with / > 2 in LoTSS (sources with limits would fill out the lower-left region; see Sec. 5.1 for more details). The traditional cut-off is seen to select a clean sample of quasars which are highly jet-dominated; however, a number of sources that do not satisfy the traditional radio-loud cut also have their radio emission dominated by the radio jets.
On the other hand, it is also clear from Figure 12 that there is a significant population of sources with 5GHz / 4400Å < 10, that would hence be traditionally classified as radio-quiet, for which the jet is still the dominant source of the radio emission. For some of these 'radio-quiet' sources, the jet produces two orders of magnitude more radio emission than the star-formation.
The other definition of radio-loudness widely used in the literature is a selection purely on the basis of radio luminosity. Figure 13 shows the fraction of sources for which the jet is significantly the dominant source of radio emission ( AGN / SF > 5), as a function of radio luminosity, in different bins of redshift. In each redshift bin, the transition from less than 10% of sources satisfying this criteria, to over 90% satisfying it, happens over typically an order of magnitude range in radio luminosity. Furthermore, the radio luminosity at which 50% of sources satisfy the criteria increases with redshift; this is likely to be at least partially driven by the increase in quasar optical luminosity with redshift due to the sample selection effects. It is clear that a simple radio luminosity cut offers only a very crude manner of selecting jet-dominated sources.

The star-formation rates of quasar host galaxies
From analysis of the SF component of our model, we find that the typical SFR increases with increasing optical luminosity of the quasar. Intuitively, this appears reasonable: if there is more gas available in the galaxy then there will be both more for the BH to accrete (higher optical quasar luminosity) and more available to form stars, in line with the Kennicutt-Schmidt Law, SFR ∝ 1.4±0.15 gas (Kennicutt 1998). We also find that the SFR of the quasars increases with redshift out to at least ∼ 2. These trends of the SFR of the quasar host with redshift and optical luminosity are seen independently of each other.
The evolution of the SFR of the quasar hosts with redshift is in line with many results in the literature which find the same increase Figure 13. The fraction of radio sources at each redshift for which the jet is the dominant source of the radio emission ( AGN / SF > 5), as a function of radio luminosity. The plot is created by convolving the radio luminosities with a kernel density estimator of width 0.2 dex. The transition of the population from star-formation-dominated to jet-dominated happens over an order of magnitude in radio luminosity in each case, and the mid-point radio luminosity depends strongly on redshift (either directly, or indirectly due to the correlation between redshift and absolute magnitude). A radio-loudness classification based solely on radio luminosity would not produce a representative sample of jet-dominated sources.
(e.g. Bonfield et al. 2011;Rosario et al. 2012;Mullaney et al. 2012;Harrison et al. 2012). This is likely to be related to more gas being available at earlier times. However, the strong correlation of SFR with optical luminosity of the quasar appears, at first glance, to be at variance with recent studies that have concluded that, once redshift effects are accounted for, any correlation of SFR with quasar luminosity is either weak (e.g. Harrison et al. 2012;Azadi et al. 2015;Stanley et al. 2017;Stemo et al. 2019) or insignificant (e.g. Rosario et al. 2012;Mullaney et al. 2012;Stanley et al. 2015). However, these studies have been typically based on moderate luminosity AGN ( bol in the range 10 35 − 10 38 W), selected through deep Xray observations in relatively small fields. The AGN in the current study are rarer, higher-luminosity quasars ( bol ∼ 10 38 − 10 40 W): previous studies which have probed to these high AGN luminosities find, like our study, that SFR does correlate with AGN luminosity in this regime (e.g. Shao et al. 2010;Bonfield et al. 2011;Rosario et al. 2012;Harris et al. 2016;Dong & Wu 2016;Lanzuisi et al. 2017); at lower redshifts, this correlation has also been seen to extend to lower AGN luminosities (e.g. Rosario et al. 2012), in line with the strong correlation in SDSS found by Netzer (2009).
Following earlier investigations, we describe the dependence of the host galaxy SFR on luminosity and redshift (for these powerful quasars) as Using the best-fit SFRs derived for each subsample in redshiftluminosity space, we find constraints on and as shown in Fig. 14.
As is evident from the upper left panel of Fig. 7, which shows that the SFRbol relation is not strictly linear, the fitting function does not provide a particularly good fit to the data across all parameter space (reduced chi-squared ∼ 7, dominated by the lower redshift lowest luminosity points). However, the fits work well at ≥ 1 and allow a like-for-like comparison of our results against previous studies. Fig. 14 includes the results obtained by Serjeant & Hatziminaoglou (2009) and by Bonfield et al. (2011). In both cases these previous studies used far-infrared estimates of the SFR, on the assumption that the far-IR luminosity is dominated by the SFR contribution. As can be seen, our results agree very well with those of Serjeant & Hatziminaoglou, and our value of = 1.61 ± 0.05 agrees with that of Bonfield et al.; their is marginally discrepent, probably due to the relatively small area of sky included in their study. A more recent Herschel-based study by Dong & Wu (2016) (also shown on Fig. 14) found = 0.46 ± 0.03, in very good agreement with our derived value of = 0.47 ± 0.01. It is also notable the scatter around the relation found by Dong & Wu is a few tenths of a dex, consistent with the value of Ψ = 0.42 dex found in our model for the width of the distribution of star formation rates at given luminosity and redshift. Overall, the close agreement between our radio-derived SF properties of the quasar hosts and those derived from far-IR data gives confidence that the weak radio emission does indeed arise from star-formation, and is not dominated by other processes such as disk coronal activity (see Panessa et al. 2019), which could also be correlated with the quasar luminosity (e.g. Laor & Behar 2008). The dependence of the SFR on quasar luminosity at a given redshift can be combined with the quasar luminosity function to investigate the overall contribution of quasar host galaxies to the cosmic star-formation rate density. We fit a linear function to the SFRbol relation in each redshift range, and combine this with the quasar luminosity function derived by Ross et al. (2013). Specifically, we use Ross et al.'s recommended double power-law model with pure luminosity evolution over the redshift range 0.3 < < 2.2 and luminosity and density evolution over the range 2.2 < < 3.5 (with parameters derived from the Stripe 82 data). We integrate the resultant luminosity function down to quasars with luminosity 0.01 * at each redshift. The results are shown in the upper panel of Fig. 15, with statistical errors determined by combining the uncertainties in our fitted SFR- Figure 15. Top: the contribution of quasar host galaxies to the cosmic starformation rate density, derived by combining the quasar luminosity function of Ross et al. (2013) with the relations between SFR and derived in Fig. 7. Bottom: the fractional contribution of quasars to the total cosmic starformation rate density. In both panels, the grey-shaded region represents the uncertainty.
bol relation with the uncertainty in the faint-end slope derived by Ross et al.: these parameters dominate the statistical error budget. Systematic errors may arise from the choice of integration limit, or from a flattening of the SFRbol relation at lower luminosities (e.g., see Lanzuisi et al. 2017), but these are more likely to shift the whole distribution vertically than to change its shape.
The SFR density in quasar host galaxies is seen to increase with increasing redshift from the current epoch back to ∼ 2, where it flattens and then declines to higher redshifts. This mirrors the overall cosmic SFR density (e.g. Madau & Dickinson 2014). To compare these, the lower panel of Fig. 15 shows the fractional contribution of quasar host galaxies to the cosmic SFR density, derived by dividing the SFR density in quasars by the functional fit to the total cosmic SFR density provided by Madau & Dickinson (2014). Madau & Dickinson integrate the cosmic SFR density down to 0.03 * SFR at each redshift, which is broadly comparable to our integration limit for the quasars, but we include an additional 0.1 dex in quadrature in the uncertainty to account for such systematic differences. Fig. 15 shows that quasar host galaxies account for approximately 0.15% of all cosmic star formation at ∼ 0.5, rising to 0.4% at ∼ 2 and then flattening towards higher redshifts. We discuss a possible explanation for these results in the following subsection.

Merger triggering of powerful quasar activity
Analysis of moderate luminosity AGN samples has indicated that they lie close to the SFR versus stellar mass 'main sequence' observed for galaxies, with the increase in their SFR with redshift simply mirroring the evolution of the star-forming main sequence (see review by Heckman & Best 2014). This, together with the lack of evidence for a higher fraction of galaxy mergers or interactions in these AGN compared to a control sample (e.g. Kocevski et al. 2012), indicates that the AGN activity in these objects is triggered by internal secular processes (see also Smethurst et al. 2019, and references therein). This is consistent with the observation that the host galaxies of moderate luminosity AGN often possess 'pseudo-bulges' (e.g. Capetti & Balmaverde 2006); these are rapidly rotating bulges, with power-law profiles and disky isophotes, which are believed to form through secular processes such as bars and disk instabilities in spiral galaxies (see Kormendy & Kennicutt 2004). Diamond-Stanic & Rieke (2012) argue that in such AGN, the SFR on nuclear scales (radius < 1 kpc) correlates well with the black hole accretion rate, as both trace the very central gas densities of the system, but the extended star-formation rates (which is all that can be measured in high-redshift systems) do not correlate strongly; indeed any residual correlation between global SFR and AGN luminosity may be associated with both properties having a mutual dependence on galaxy mass (e.g. Stemo et al. 2019).
In contrast, at the highest quasar luminosities ( bol ∼ > 10 38 W), and especially towards lower redshifts, the observed SFRs of the quasar hosts lie significantly above the star-forming main sequence, indicating that these objects are associated with starburst activity. It has been widely argued that major galaxy mergers are the most likely origin of both the starburst and the associated powerful quasar (e.g. Sanders et al. 1988;Hopkins et al. 2006); these objects are typically hosted by very massive ellipticals, regardless of whether they are radio-loud or radio-quiet (e.g. Dunlop et al. 2003;Pagani et al. 2003). Some observations also find direct evidence of a high merger fraction at these high luminosities (e.g. Treister et al. 2012;Goulding et al. 2018); this is particularly the case for the luminous reddened quasars (Glikman et al. 2015), with some authors suggesting that these are an evolutionary phase during which the dust associated with the merger starburst is blown out of the galaxy, before unreddened quasars are observed (e.g. Calistro Rivera et al. 2021). It must be noted, however, that other authors have argued that the role of major mergers is sub-dominant (Hewlett et al. 2017) or unimportant (Marian et al. 2019) even at the highest luminosities, and therefore the requirement for major mergers to trigger the most luminous quasars remains controversial. If these objects are triggered by major mergers though, then the available gas mass and the dynamical time of the system would influence both the SFR and the black hole accretion rate, leading naturally to a correlation between these two properties, as we observe in our study. Lamastra et al. (2013) argue that starburst activity becomes increasingly important (relative to quiescent star-formation) at earlier cosmic times: they find that the fraction of the total cosmic SFR density that is associated with starbursts increases by a factor of 4 between ≈ 0.1 and ≈ 5. Similarly, Martin et al. (2017) find a factor-of-two increase from = 0 to ∼ 2 in the fraction of cosmic SF that is associated with merging systems. From a theoretical standpoint, the simulations of Hopkins et al. (2010) also suggest that merger-induced starbursts contribute 1-5% of all star formation at ∼ 0, rising to 4-10% at > 1 and flattening at higher redshift. In Fig. 15 we found that the fraction of star formation in powerful quasar host galaxies increased by a factor 2-3 between = 0.5 and ∼ 2 and then flattened; this mirrors the trends found for starbursts and mergers, as would be expected for the picture where powerful quasars are triggered by galaxy mergers.
Finally, in Section 4.3 we considered the influence of black hole mass on the results that we find. We showed in Fig. 10 that for the bulk of the population at ∼ < 2, the SFR -AGN luminosity relation is the same in both higher and lower black hole mass bins. This is to be expected if these systems are driven by major mergers, since the black hole itself does not play a major role in determining how quickly gas will be funnelled down on to it. At > 2, however, if the possible small increase in the SFR in higher black hole mass systems compared to lower black hole masses is indeed real, then this could be understood if, at those redshifts, we are reaching a regime where the gas fractions in massive galaxies are high enough that the SFRs and optical luminosities that we observe can be achieved in galaxies lying on the star-forming main sequence; in such galaxies, a higher black hole mass is likely to be correlated with a higher stellar mass and hence a higher SFR.

Quasar lifetimes
Numerically, we can convert the quasar optical luminosity to an estimated growth rate of the central SMBH, BH (as in Section 2.4.2). Although the bolometric correction from optical luminosity is more uncertain than that from e.g. X-rays, relying on an empiricallyderived correlation, it is sufficiently accurate for our purposes.
In Figure 7, we find that the ratio of the star-formation rate in the host galaxy to the growth rate of the SMBH is typically Ψ/ BH ∼ 40 (varying from ∼ 20 − 100 across different redshift and luminosity bins). This ratio is an order of magnitude lower than the ratio of bulge mass to BH mass in present-day bulges ( bulge / BH ≈ 500 − 700; e.g. Marconi & Hunt 2003;Häring & Rix 2004), suggesting that the phase of quasar activity must be an order of magnitude shorter than the duration of the star-formation activity of the galaxy. This result is also consistent (for a merger-triggered scenario) with the ratio between the fraction of the cosmic SFR density in quasar host galaxies (∼ 0.4% at z∼ 2; Fig. 15) and that in merger-induced starbursts (∼ 5% at ∼ 2; Hopkins et al. 2010). In major mergers, the peak starburst activity is understood to last for a few tens of Myr (e.g. Genzel et al. 1998;Bernloehr 1993;Mihos & Hernquist 1994); this then implies quasar lifetimes consistent with the lower end of the range typically suggested by observations; ≈ 10 6 − 10 8 yr (Martini 2004). Furthermore, it is also interesting to note that we observe lower values of Ψ/ BH at the highest accretion rates, suggesting that the lifetimes of the most luminous quasars may be shorter than those of lower luminosities. This is reasonable, since the higher luminosity quasars will be consuming their gas supply more quickly.
There is an important caveat to these conclusions: it may well be the case that the peak of quasar activity does not correspond to the peak of star-formation activity, particularly if these quasars are part of a merger-driven evolutionary sequence from ultra-luminous infrared galaxies (starbursts) to quasars (Sanders et al. 1988). Wild et al. (2010) studied the growth of BHs in galactic bulges in which strong bursts of SF have recently occurred. They find that the black hole growth peaks around 300 Myr after the burst in star-formation. They propose that BH growth has been driven primarily by slow stellar ejecta from intermediate mass stars (cf. Norman & Scoville 1988), and that at earlier times the black hole growth is suppressed by supernovae feedback. These effects would mean that the currently estimated value of Ψ may under-estimate the star formation that occurs during the whole burst, and would then permit longer quasar lifetimes. The quasar luminosity may also vary over the quasar lifetime. However, any such variations are unlikely to affect the qualitative conclusion that quasar lifetimes are shorter than the period of star-formation activity, and that the most luminous quasars have the shortest lifetimes.

The powering of quasar jets
No consistent discernible trend with redshift is observed for the normalisation of the jet power, . We do see hints of a decrease in with increasing redshift for the highest optical luminosities (see Fig. 7), in agreement with some previous studies (e.g. Jiang et al. 2007;Baloković et al. 2012;Kratzer & Richards 2015); this is where the radio sources tend to be the most luminous and extended, and might therefore be due to the increasing importance of inverse Compton scattering losses in these sources towards higher redshifts, as suggested by Gürkan et al. (2019). Overall, however, the lack of strong redshift dependence suggests that the physical properties that govern the power of the radio jet are probably local properties of the system. In agreement with the aforementioned studies, we do find strong evidence of an increase in the fraction of sources at high radio luminosities with increasing optical luminosity (or black hole accretion rate): we find that increases by just over an order of magnitude as the optical luminosity brightens by 4 magnitudes, corresponding to ∝ bol with ∼ 0.65. This is broadly in line with the exponent of 0.85 suggested by White et al. (2007) from analyses with the much shallower radio data from the Faint Images of the Radio Sky at Twenty centimetres survey (FIRST; Becker et al. 1995). We emphasize this optical luminosity dependence is a scaling of a full power-law distribution of jet powers to typically higher powers: at all optical luminosities a large range of jet powers is seen.
We also find weak evidence that the fraction of sources at high radio luminosities may increase with an increase in the mass of the central SMBH, at least out to = 2. Dependencies of the radio luminosity, or RL fraction, on the BH mass have been argued by previous studies (e.g. Laor 2000;Lacy et al. 2001;Dunlop et al. 2003;McLure & Jarvis 2004;Best et al. 2005). It is possible that any such dependence of jet power on black hole mass could be due to a residual correlation with other properties, such as stellar mass of the host. Sabater et al. (2019) found (albeit for radiatively-inefficient AGN but the argument is the same) that although the fraction of radio AGN increases with increasing BH mass, once they disentangled the correlation between BH and stellar mass, they found that the fraction of radio AGN was mainly driven by the stellar mass. As we do not have information about stellar masses for our sample, it is difficult to draw any direct conclusions as to which property most directly drives any increase in jet power.
Our results are unable to provide direct evidence for the physical mechanisms that produce quasar jets, however, they do allow us to speculate. First, as discussed in Section 5.1, our model suggests that radio jets are ubiquitous in quasars -or at very least they exist down well into the traditional 'radio-quiet' regime where the jet luminosity is comparable to the starburst luminosity (or other sources of radio emission, cf. Panessa et al. 2019). This result is consistent with the detection of jet-like structures in deep high-resolution observations of such quasars (e.g. Hartley et al. 2019). This suggests that the jet-launching mechanism operates in all quasars, but with different powering efficiency. Second, as discussed in Section 5.2, we interpret that the majority of the (high luminosity) quasars in our sample are triggered by galaxy mergers. Third, we observe trends of the jet power normalisation with optical luminosity (or SMBH growth rate) and BH mass (or perhaps stellar mass).
One popular explanation for the varying power of radio jets is a dependence on black hole spin (Blandford & Payne 1982). However, our observations suggest a very wide range of jet powers, which would require a correspondingly wide range of black hole spin parameters: numerical models find it hard to produce this (Volonteri et al. 2013). Furthermore, if the majority of our quasars (of high and low jet powers) are indeed triggered by mergers, then since the orbits of two colliding BHs would give rise to a significant amount of angular momentum in the resulting BH, it would be surprising to find such a large population of radio-quiet quasars. It therefore seems unlikely that variations in black hole spin can be the main factor influencing radio loudness.
An alternative hypothesis was put forward by Tchekhovskoy et al. (2011) and developed by  and Sikora & Begelman (2013): that the main parameter driving the wide range of jet production efficiencies is the magnetic flux threading a spinning black hole. Sikora & Begelman argued that periods of hot accretion are efficient at depositing magnetic flux close to the black hole. Therefore, if a period of (cold-accretion) quasar activity had been preceded by a period of hot accretion, as might happen for example when a giant elliptical galaxy undergoes a merger with a disk galaxy, then very powerful radio jets would be likely to result. We have argued that the quasars in our sample are predominantly triggered by major mergers, and hence in all cases these should have the spinning black hole required. These mergers will have involved a mix of progenitor galaxies (disk-elliptical and disk-disk mergers) and so the black holes may be threaded by a wide variety of magnetic fluxes, producing a significant range of jet efficiencies. Furthermore, as well as high power radio jets produced in these 'magnetically-choked accretion flows', Sikora & Begelman (2013) predict that low power or intermittent jets in quasars will arise from fluctuating magnetic fields arising in the corona above a thin accretion disk, or in a hot inner region of an accretion flow. Combined, these could give rise to the continuous distribution of jet powers that our model adopts. Finally, for these accretion flows  show that (for fixed other parameters) the jet power increases with the accretion rate (optical luminosity), broadly in line with the correlation that we observe. Thus, our observations can all be well explained by this model.

CONCLUSIONS
We present a model of the radio luminosity distribution of quasars that assumes that the radio emission of every quasar is a superposition of two components: active galactic nuclei (jets) and star-formation. We compare Monte Carlo simulated samples to our sample of ∼ 42, 000 quasars from the Sloan Digital Sky Survey quasar catalogue fourteenth data release with the radio emission measured by the LOFAR as part of the first data release of the LoTSS.
We find that our two-component model is valid in describing the observed radio emission across a wide range of redshifts and optical quasar luminosities. We therefore argue that an intrinsic bimodality in the radio loudness distribution of quasars does not exist; instead, 'radio-loud' quasars are simply the luminous tail of a continuous jet power distribution. Our analysis cannot prove whether or not quasar jets are ubiquitous, but our results do suggest that the radio emission from jets is contributing down at least to the level that radio emission from SF becomes comparable. Our model naturally leads to the expectation that some radio-quiet quasars will have their radio emission dominated by small-scale jets, and others by star-formation, in line with observations.
Given the validity of our model, we investigate how the parameters of our model depend on redshift, optical luminosity (which we relate to the SMBH growth rate) and BH mass. The width of the Gaussian SF component and slope of the power-law AGN component are found not to vary significantly so we fix the two parameters in our model. We find a strong correlation between the mean SFR and the optical quasar luminosity, with ∝ 0.47±0.01 bol . These results are in line with other recent studies that probe the high quasar luminosities characateristic of our sample, and with far-infrared determinations. Unlike at lower quasar luminosities (where the AGN activity is believed to be triggered by secular processes and the host galaxies lie close to the star-forming main sequence, leading to little correlation between star-formation and BH acrretion rate), these high luminosity quasars are understood to be triggered by massive galaxy mergers, where the gas fraction and dynamical time of the system will influence both star formation and black hole accretion, leading naturally to the observed correlation. The ratio of the black hole growth rate to star formation rate are observed to be an order of magnitude higher in these quasars than the current black hole mass -bulge mass ratio, implying that quasar activity must be an intermittent phase.
We also investigate the cosmic star-formation history of quasar host galaxies. We see an increase in the SFR of the quasar hosts out to ∼ 2, which then flattens in the range ∼ 2 − 3. The integrated starformation rate density in quasar host galaxies contributes roughly 0.15% of the total cosmic star-formation rate density at ∼ 0.5, increasing to around 0.4% at ∼ 2 and then flattening. This trend mirrors that of the importance of merger-induced starbursts to cosmic star formation. We observe little dependence of the SFR on BH mass until we get to the highest redshifts, at which point we see weak evidence that increasing the BH mass may correspond to a small increase in the SFR. These highest redshift quasars lie closer to the star-forming main sequence and some may be secularly triggered, in which case this observation would be naturally explained by a correlation of both parameters with the stellar mass of the host.
The normalisation of the jet power distribution is shown to have little dependence on redshift, suggesting that the physical properties responsible for producing powerful radio jets are local to the system. We do observe an increase in the fraction of sources at high radio luminosities with increasing optical luminosity (or BH growth rate) and indications of an increase with BH mass, in line with previous studies. Although our results do not allow a definitive answer to be reached on the physical mechanisms that produce radio jets, by considering the possible interpretations of our results, we conclude that the model which can best explain our combination of results is the one of Sikora & Begelman (2013), where the magnetic flux threading the black hole is the primary factor influencing jet production efficiency. This model is able to naturally produce the very wide range of radio loudness required by observations, while also giving a jet power to optical luminosity correlation. We use our model to investigate the effectiveness of different literature definitions of 'radio loudness' for quasars. We find that the traditional radio-loudness selection based on the ratio of radio-to-optical luminosities, = 5GHz / 4400Å > 10, cleanly selects a sample of jet-dominated sources ( AGN / SF > 1), but does so in a substantially incomplete manner: many quasars classed as 'radio-quiet' by this criteria have the majority of their radio emission associated with the jet, which can be up to two orders of magnitude brighter than that from star formation. We find that definitions of radio-loudness based solely on radio luminosity perform relatively poorly.
The potential of the Monte Carlo approach that we have adopted, particularly if subsequently adapted to use a Bayesian framework, is vast for this research. Wider areas of the LOFAR survey (Shimwell et al., in prep.), and new, much lower noise data in the LOFAR Deep Fields (Tasse et al. 2020;Sabater et al. 2020) can be seamlessly added to the existing sample of quasars, to constrain the parameters of the model further. With the larger samples, additional parameters can be investigated, such as the quasar colour recently studied by Klindt et al. (2019) and Rosario et al. (2020), to help disentangle evolutionary effects. Further work to disentangle the dependencies on other relevant properties of the system such as stellar mass and the Eddington ratio, would also be informative.

APPENDIX A: THRESHOLD TESTS
Since the definition of the threshold flux density is arbitrary, it is important to test that our results do not depend on our choice of thresh . To do this, we computed the best-fit values of the model parameters for given subsamples of − space for a range of threshold flux densities. The results for one such subsample are shown in Fig. A1; similar results are found for other subsamples. As can be seen from Fig. A1, the best-fit values for the parameters of the model are stable, within the error bars, almost right down to a threshold flux density that lies deep in the main quasar population at thresh = 0.15mJy. Therefore, a threshold flux density of thresh = 1mJy should be suitable for our analyses.

APPENDIX B: TABLES OF BEST-FIT VALUES FOR MODEL PARAMETERS
In this appendix, we provide tables of best-fit values for the model parameters when fitted to the full sample (Table B1), and when split by black hole mass (Table B2).

APPENDIX C: M I − Z BINNING TESTS
It is important to test that the observed SFR trend with redshift is not due to selection effects related to the binning in − space, coupled with the strong dependence of the SFR on . Each bin has a width of Δ = 1 and so, for a given slice of (e.g. −24 < < −23 across all ), the correlation between and imprinted by selection effects means that higher redshift bins within the slice may contain quasars which are typically more luminous than those in lower redshift bins. In turn, this could give rise to an apparent trend between SFR and , arising solely due to binning biases.
To test this, the distribution of the mean absolute i-band magnitude, ,mean , in each slice of was plotted as a function of redshift, offset by the central absolute i-band magnitude of the slice. The Table B1. The variation of the best-fit values of the model parameters as a function of redshift and the optical luminosity of the quasar, for the full sample. Columns 1 and 2 indicate the redshift and optical luminosity bin under consideration. Column 3 gives the black hole accretion rate corresponding to the optical luminosity (see Section 2.4.2). Columns 4-13 give the best-fit values for the model parameters, and their errors, as plotted in Figure 6: specifically the mean luminosity ( ) associated with the star-forming component, the corresponding star-formation rate (Ψ), the width of the Gaussian distribution ( Ψ = 0.93 ), the power-law slope of the AGN component ( ) and its normalisation term ( ). Δ values indicate the relevant uncertainties. Columns 14-17 give the best-fit values when and are fixed (Section 4.2), as displayed in Figure 7 Table B2. The variation of the best fit star-formation rate (Ψ) and jet power normalisation parameter ( ) as a function of redshift and the optical luminosity of the quasar, split by black hole mass. Columns 1 and 2 indicate the redshift and optical luminosity bin under consideration. Columns 3-6 give the best-fit values for the two model parameters, and their errors, for black hole masses below the median value in each bin, while columns 7-10 give the equivalent values for black hole masses above the median value. The data in this table are plotted in Figure 10.  Figure A1. The four panels show how each of the best-fit values for the model parameters depends on the threshold flux density used to split the subsample in order to implement the method discussed in Section 3.4. The subsample shown is the same one as shown in Figure 3. The plot demonstrates that the results obtained are robust against the choice of flux threshold to well within the errors; a similar conclusion is found for other subsamples. results of this test can be seen in Fig. C1. We do observe a relatively small increase in ,mean , particularly at lower luminosities and especially towards the upper end of each of their respective redshift ranges. However, from Figure 7, such changes in are only 0.2-0.3 magnitudes at most, which correspond to a difference of only ≈ 0.1 dex in SFR; this is much smaller that the variations seen in SFR with redshift. With this bias being relatively small, and with evidence of the SFR in the range −26 < < 25 flattening off despite the bias at high , we can infer the derived results such as in Fig. 7 are not significantly affected by the bias. This paper has been typeset from a T E X/L A T E X file prepared by the author.