A novel Bayesian approach for decomposing the radio emission of quasars: I. Modelling the radio excess in red quasars

Studies show that both radio jets from the active galactic nuclei (AGN) and the star formation (SF) activity in quasar host galaxies contribute to the quasar radio emission; yet their relative contributions across the population remain unclear. Here, we present an improved parametric model that allows us to statistically separate the SF and AGN components in observed quasar radio flux density distributions, and investigate how their relative contributions evolve with AGN bolometric luminosity ($L_\mathrm{bol}$) and redshift ($z$) using a fully Bayesian method. Based on the newest data from LOFAR Two-Metre Sky Survey Data Release 2, our model gives robust fitting results out to $z\sim4$, showing a quasar host galaxy SFR evolution that increases with bolometric luminosity and with redshift out to $z\sim4$. This differs from the global cosmic SFR density, perhaps due to the importance of galaxy mergers. The prevalence of radio AGN emissions increases with quasar luminosity, but has little dependence on redshift. Furthermore, our new methodology and large sample size allow us to subdivide our dataset to investigate the role of other parameters. Specifically, in this paper, we explore quasar colour and demonstrate that the radio excess in red quasars is due to an enhancement in AGN-related emission, since the host galaxy SF contribution to the total radio emission is independent of quasar colour. We also find evidence that this radio enhancement occurs mostly in quasars with weak or intermediate radio power.


INTRODUCTION
Quasars (QSOs) are known to have a profound impact on the evolution of their host galaxies, both through radiation-driven processes and through jet feedback (e.g.see reviews by Fabian 2012;Heckman & Best 2014).These physical processes and their impact on the host galaxies can be traced across multiple wavelengths, therefore helping us build knowledge towards the evolutionary paths of galaxies (e.g.Kormendy & Ho 2013).
Radio observations play a particularly important role in analysing the impact of active galactic nuclei (AGNs) on galaxy evolution, since they are a major tracer of AGN jets (Heckman & Best 2014;Hardcastle & Croston 2020) , in which relativistic electrons create synchrotron radiation that is detectable in radio bands.AGN jets have a highly collimated structure originating from an optically thick launching region also known as the jet base (Blandford & Königl 1979;Reynolds 1982), and give rise to radio lobes with a ★ E-mail: bohan.yue@ed.ac.uk (BY) steep continuum profile due to highly relativistic electrons.The nonthermal processes from relativistic electrons give the radio cores a high brightness temperature (Blundell & Beasley 1998), and radio emission from a jet origin is often highly polarized.These properties can be used to identify radio jets even in the lowest power regime, given sufficient angular resolution (e.g.see the resolved detection of NGC 4151 in Pedlar et al. 1993;Williams et al. 2017;Mundell et al. 2003).
Radio emission can also arise from regions of star formation (SF; e.g.Condon 1992), which could be associated with the AGN activity.Radio emission from SF is characterised by its steep spectrum at low frequencies (with a spectral index1 of  ∼ 0.7), which is caused by the acceleration of electrons in supernova remnants.
Separating SF and AGN features in radio-faint quasars is particularly difficult, as the hosts of most AGNs are star-forming galaxies (Heckman & Best 2014).The most popular approach where suitable multiwavelength optical / IR data are available is through spectral energy distribution (SED) fitting (e.g.Calistro Rivera et al. 2017;Delvecchio et al. 2017;Whittam et al. 2022;Best et al. 2023); however, the number of sources with reliable SF-AGN identification is limited by the scope of optical/IR surveys and the number of available photometric bands.Recently, without having to utilize the multiwavelength data, Morabito et al. (2022) applied very-long baseline interferometry (VLBI) techniques to calculate the radio brightness temperature and separate SF and AGN contributions in spatiallyunresolved quasars, but similar tasks are still difficult to accomplish for sources without VLBI-resolution observations.Our lack of knowledge about the ongoing physical processes driving the radio emission in AGNs and their host galaxies has a profound impact in addressing some of the most important questions regarding the quasar properties and the impact of the quasars on their host galaxies.
Historically, quasars have been separated into two categoriesradio-quiet (RQ) quasars and radio-loud (RL) quasars -based on their radio-loudness R, which is typically defined as the ratio between fluxes in the optical and radio bands:  =  (4400 Å)/  (6cm) (Kellermann et al. 1989).Results from large radio surveys indicated an apparent dichotomy between RL and RQ quasars, thanks to an asymmetric distribution of radio flux densities with a long tail towards the radio-bright end (due to powerful jet activities), combining with a peak at low flux density regime (due to host galaxy star formation or small-scale radio emissions from the AGN).However, debates are still open on whether star formation in host galaxies can provide sufficient radio emission as observed in RQ quasar samples (for a thorough review on the topic, see Panessa et al. 2019): some studies found that the radio emission in RQ quasars can be explained by SF alone (e.g.Kimball et al. 2011;Condon et al. 2013;Bonzini et al. 2013), while others suggest that the majority of the emission needs to come from AGN activity (e.g.Zakamska et al. 2016;White et al. 2015White et al. , 2017)), in the form of small-scale jets, AGN winds (Mullaney et al. 2013;Zakamska et al. 2016;Morabito et al. 2019;Petley et al. 2022), or accretion disk coronae (Laor & Behar 2008;Chen et al. 2023).In support of the latter argument which links RQ quasar radio emission with AGN activities, weak radio jets have recently been identified within several spatially resolved RQ quasars (e.g.Leipski et al. 2006;Herrera Ruiz et al. 2016;Jarvis et al. 2019).These studies suggest that the weak jet activity in RQ AGNs might be merely a scaled-down version of RL jets, and the only difference lies in the powering efficiency of acceleration on sub-parsec scales, which therefore would indicate that radio jets may contribute to the radio emission in all quasars, whether they are low-luminosity small-scaled or radio-loud extended.
The LOw-Frequency ARray (LOFAR; van Haarlem et al. 2013) is a state-of-art radio telescope observing at 120-168 MHz with its highband antennae.Thanks to LOFAR's wide field of view combined with its overall sensitivity, including sensitivity to low surface brightness emission, we can now detect radio sources at an order of magnitude higher sky density than any previous large-area radio surveys and obtain deep radio images for a large sample of quasars with optical counterparts (see Section 2.1 for details).
Motivated by the high sensitivity of LOFAR observations, Macfarlane et al. (2021, hereafter M21) studied in detail the radio flux density distribution of quasars from the Sloan Digital Sky Survey (SDSS; e.g.Pâris et al. 2018).They adopted the conclusion from Gürkan et al. (2019) that every source hosts a contribution from both jet activity in the AGN and star formation in the host galaxy, and proposed a two-component model that characterises the radio emission in the quasar population, where both the SF and AGN components were modelled from physical descriptions.The contributions to the overall quasar radio emission from host galaxy SF and AGN jets can thus be studied independently with the two-component model.While investigating the variation of radio emission from SF and jet components with bolometric luminosity, redshift, and black hole (BH) mass, their two-component model was able to provide a good fit to the data across all parameter space, thus strongly indicating the lack of an RL/RQ dichotomy.
As more evidence point toward a continuous distribution of quasar jet power rather than a dichotomy (e.g.Cirasuolo et al. 2003a,b;Baloković et al. 2012), which factors affect the powering efficiency of radio jets has thus become a more interesting topic for discussion.Some studies have argued that jet strengths are related to the bolometric luminosity of the AGN; earlier studies show a higher jet fraction in low-luminosity RQ quasars (e.g.Ho & Ulvestad 2001;Blundell & Rawlings 2001;Ulvestad & Ho 2001), while in the M21 model the fraction of sources with powerful jets (denoted as  ) does increase with  bol sublinearly (  ∝  0.65 bol ).Others have found strong dependencies between jet power and other parameters such as black hole mass (e.g.Laor 2000;Lacy et al. 2001;McLure & Dunlop 2004;Best et al. 2005).Some recent works (including Retana-Montenegro & Röttgering 2017) have found higher angular clustering in RL quasars, suggesting that larger-scale environment might also be a factor.Morabito et al. (2019) found a lower RL fraction in broadabsorption line quasars (BALQSOs) which links the radio-loudness to an outflow phase.Based on their two-component model that separates the host galaxy contribution from the observed quasar radio emission, M21 has revealed positive correlations between typical jet power and optical luminosity or black hole mass -but not with redshift.Therefore, they suggested the production of radio jets is more likely to be governed by intrinsic properties.
Recently, it has been shown that a small population of unusually red quasars (rQSOs; e.g.Richards et al. 2003) shows a significant excess in radio emission compared to a control sample of QSOs with blue or average colours (e.g.Glikman et al. 2007;Urrutia et al. 2009;Glikman et al. 2012;Klindt et al. 2019;Calistro Rivera et al. 2021) .Several theories have been proposed to explain the nature of this rQSO population.Netzer (2015, and references therein) argued that from the standard AGN model's point of view, the red quasars are simply typical blue quasars with their core area (the accretion disc and broad line region) partially obscured by the dusty torus.However, more recent studies have connected red quasars with other phenomena including flatter bolometric luminosity functions (e.g.Banerji et al. 2015), occurrence of major mergers (e.g.Urrutia et al. 2008;Glikman et al. 2015) and higher incidence of strong AGN outflows (e.g.Urrutia et al. 2009;Banerji et al. 2012).These phenomena cannot be explained by the torus obscuration (see Calistro Rivera et al. 2021).More specifically, based on the higher prevalence of radio activity in the SDSS-selected red quasar sample, Klindt et al. (2019) concluded that red quasars are a fundamentally different population from the typical blue quasars (see also Rosario et al. 2020;Fawcett et al. 2020;Calistro Rivera et al. 2021;Fawcett et al. 2022;Calistro Rivera et al. 2023).This evidence supports the quasar evolution model proposed by Sanders et al. (1988), where the dichotomy between the red and blue quasars arises during an evolutionary phase that connects dust-rich star formation and AGN activity through gas feedback between AGN and its host galaxy.According to the model, within the red quasar phase, the wind/outflow of the central black hole gradually drives away the obscuring dust generated during a period of fast star formation activity (a starburst phase, perhaps driven by merger activity).Eventually, the outflow shuts down SF and reveals the unobscured central black hole, which appears as one of the typical blue quasars (see also Hopkins et al. 2006Hopkins et al. , 2008;;Farrah et al. 2012;Glikman et al. 2012).
Unfortunately, these theories lack direct evidence to validate themselves.Klindt et al. (2019) and Fawcett et al. (2020) found that the excess in radio emission of red quasars is mostly seen in compact and radio-faint systems, which often lie around the traditional threshold between radio-quiet and radio-loud sources.Rosario et al. (2020) used more recent radio data from the LOFAR Two-metre Sky Survey data release 1 (LoTSS DR1; Shimwell et al. 2019) and found the modelled SF contribution to the total radio emission shows little difference between red quasars and blue quasars, thus concluding that the reddening is likely linked to AGN activities within the system.This argument is further supported by Rosario et al. (2021) where they used high resolution e-MERLIN data to show that the radio emission in red quasars is more extended on the most compact scales, indicating a greater AGN contribution to the radio emission of red quasars compared to blue quasars.Fawcett et al. (2022) used data from the X-shooter spectrograph and argued that a dusty environment can fully explain most of the colour differences between red and blue quasars, while the radio excess is more likely to be connected with jet interactions in a higher-opacity interstellar medium (ISM)/circumnuclear environment rather than accretion disk activities or outflows.However, all of these studies were unable to separate the contribution from SF and AGN activities due to the shortage of sufficiently high quality radio data.
In this work, we adopt the assumptions from the two-component model in M21 and propose an improved parametric model using a Bayesian approach that can let us study radio emission from the host galaxy SF and AGN jet activity independently within a wider parameter space while obtaining robust results on the possible correlations.Increased quasar samples in the LoTSS DR2 (Shimwell et al. 2022) and SDSS DR16Q (Lyke et al. 2020) catalogues allow us to further investigate the contribution to the quasar radio emission by any physical processes, including those associated with the quasar colour.As a result, we can finally provide a quantified view of the leading powering mechanism of RQ quasars and the nature of the radio excess in red quasars.
This paper is structured as follows: Section 2 describes the data we used to build a quasar model with multiband measurements; Section 3 explains the parametric model we proposed to characterise quasar radio emission and the validation of our model; Section 4 presents our improved result on the dependence of quasar radio emission on optical luminosity and redshift; Section 5 shows our result and discussion on the origin of radio excess in red quasars; finally, a summary of our conclusions can be found in Section 6.Throughout this work, we assume a ΛCDM cosmology with parameter values published in the WMAP9 result (Hinshaw et al. 2013).

LOFAR Two-metre Sky Survey (LoTSS)
LoTSS (Shimwell et al. 2017(Shimwell et al. , 2019(Shimwell et al. , 2022) ) is the LOFAR HBA (highband antenna) wide-field imaging survey that aims to cover the entire Northern sky in the 120-168 MHz radio band, with a target sensitivity of ∼ 100Jy beam −1 RMS, an angular resolution of 6 ′′ and a positional accuracy of < 0.2 ′′ .To date, the most complete catalogue is the LoTSS DR2 catalogue published by Shimwell et al. (2022), spanning over 5,720 deg 2 of sky area and reaching a median sensitivity of 83Jy beam −1 .The LoTSS DR2 survey has ∼ 10 times more coverage compared to the previous data release, DR1 (Shimwell et al. 2019) 1995), LoTSS DR2 reaches ∼ 10 times better sensitivity, assuming a spectral slope of  = 0.7.In total, the LoTSS DR2 catalogue contains ∼ 4, 400, 000 radio-detected sources.The sky coverage (compared with LoTSS DR1 and the SDSS DR16 quasar catalogues) is shown in Fig. 1.Most of the sky area in LoTSS DR2 is also covered by SDSS DR16 quasar catalogue (described in the section below), providing an ideal combination for multi-wavelength studies.
Radio sources in the LoTSS DR2 catalogue were extracted from LoTSS images using the Python Blob Detector and Source Finder (PyBDSF; Mohan & Rafferty 2015), identifying sources with peak radio flux densities above the 5 limit of the LoTSS DR2 images.While PyBDSF detects regions of radio emission, they are not always correctly grouped into physical sources.To prevent wrong association of the radio sources by the PyBDSF classification, both statistical techniques and extensive visual inspection (using LOFAR Galaxy Zoo) have been used to ensure the radio catalogue represents the actual distribution of the radio sources (Hardcastle et al. 2023).In addition to the updated radio catalogue, Hardcastle et al. present the cross-matching of LoTSS DR2 sources with optical-infrared counterparts from WISE (Wright et al. 2010) and DESI Legacy Imaging (Dey et al. 2019) surveys using a method similar to that described in Williams et al. (2019) and Kondapally et al. (2021).

Sloan Digital Sky Survey (SDSS) quasar sample
Optical data for the sample of quasars is drawn from SDSS DR16 quasar catalogue (DR16Q; Lyke et al. 2020).The DR16Q catalogue is created based on observations from SDSS-I/II/III and IV epochs, reduced using the final eBOSS SDSS reduction pipeline (v5_13_0), and spectrally confirmed using the criteria provided in Table 1 of Pâris et al. (2018).It also includes previously detected sources from the DR14Q, DR12Q and DR7Q catalogues (Ross et al. 2012;Schneider et al. 2010).For quasars with redshift between 0.8 <  < 2.2, the quasar classification is done via the decision tree presented in Dawson et al. (2016); for quasars with redshift  > 2.2, the classification uses Ly forest measurements in Myers et al. (2015) instead.Our quasar optical properties, including i-band magnitude and ( − ) colours, show no signs of variation between different selection methods at  = 2.2.The total number of detected quasars is > 480, 000 and > 239, 000 for the two categories, respectively.Additional optical properties of quasars from the SDSS data, including absolute i-band magnitude2 are supplemented by survey data including GALEX, UKIDSS, WISE, FIRST, 2MASS, XMM-Newton and Gaia through cross-matching between catalogues.

Building a LoTSS-SDSS quasar sample
We aim to build our quasar sample set by extracting LoTSS radio flux density measurements for parent SDSS DR16Q quasar samples; we therefore cross-match the extracted radio data with SDSS optical data based on the quasar sky positions.Before creating our cross-matched catalogue, we applied several extra restrictions to the parent SDSS dataset; these are the same cuts that were applied to the LoTSS DR1 sample in M21: (i) Sources with absolute i-band magnitude brighter than -40 were discarded, since they were likely artefacts from the SDSS pipeline.
(ii) Sources with redshift  > 4 were discarded due to the small number statistics and larger possibility of quasar misidentification and/or erroneous redshift measurement.
(iii) As shown on Fig. 1, only SDSS sources that fell in the 13h field of LoTSS DR2 were included.Sources in the 0h field have a systematically larger uncertainty in radio flux density compared to the 13h field, and are therefore discarded to ensure the conformity of our sample.Other SDSS quasars were excluded if they fall outside the LoTSS coverage, either because of being outside the target field or being in the gaps between the LoTSS mosaics.
(iv) Finally, another 199 sources in the SDSS DR7Q and DR14Q catalogues were removed; they were selected in the SDSS catalogue only because of their radio emission, and removing them would prevent selection bias toward radio-loud sources outside the SDSS colour-selection region.
The final sample consists of 361,123 quasars, with a sample size nearly 10 times larger than the 42,601 quasars used in M21.The samples are characterized by their distribution in the i-band magnitude (  )-redshift () plane, as shown in Fig. 2. While the sample is naturally biased towards brighter magnitudes at higher redshifts, it still maintains a good dynamical range across the entire parameter space.

Radio Flux Densities from LoTSS
To obtain LoTSS radio flux density measurements for the selected quasars, we adopted different strategies as stated below: (i) We cross-matched their sky positions in the SDSS catalogue with those in the LoTSS catalogue.This will include quasars from our parent sample with radio flux densities above the LoTSS DR2 5 detection limit.Based on the test result from M21, we have set the matching radius to 1.5 ′′ for sources above 5 detection in LoTSS, in order to balance between the selection rate and random contamination.It is worth noting that the coordinates from the crossmatched optical catalogue of Hardcastle et al. (2023) were used whenever possible for the LoTSS positions, since these coordinates  are more accurate than the flux-weighted radio positions obtained from the PyBDSF code.We have successfully cross-matched 47,902 quasars from this step, all of which are included in our catalogue.
(ii) For quasars with radio flux densities below the LoTSS 5 detection limit and which thus do not have a match in the LoTSS DR2 catalogue, we performed a forced photometry on LoTSS mosaics using the method described in Gloudemans et al. (2021).This step is motivated by Roseboom & Best (2014) and M21 where considerable information has been extracted from sources with detected radio emission dominated by noise.The extracted radio flux density is assumed to be the absolute peak pixel value within the 3px × 3px box (where pixels in the LoTSS mosaic are 1.5 arcsec in size, compared to a beam size of ≈6 arcsec) surrounding the corresponding position of the source in the SDSS catalogue; the flux density uncertainty is determined by the standard deviation of pixel values in a cutout region of 100px×100px surrounding the central pixel.By performing this extraction, we assume the sources undetected by LoTSS are all compact, spatially-unresolved sources; this assumption is reasonable since these radio-quiet sources are likely to be tracing star-formation on galaxy scales, compact radio cores and/or small scale jets -all of them being relatively small-scaled compared to the 6 ′′ beam size of LoTSS survey.In fact, only 814 out of 361,123 targets have multiple radio components in the LoTSS DR2 images according to the final optical cross-matched LoTSS DR2 catalogue (Hardcastle et al. 2023); they are also among the most radio-bright sources in our catalogue, therefore having minimal effect on our results.We have extracted radio flux densities of the remaining 313,221 quasars using this process, bringing the total sample size to 361,123.

Bolometric Luminosities
We follow the approach in M21 and use the absolute -band magnitude (  ) as a proxy for bolometric luminosity of the quasar, and hence the black hole accretion rate.
To convert   into bolometric luminosity, it can first be converted into the absolute -band magnitude (  ) using the empirical relationship in Richards et al. (2006):   −   ( = 2) = 0.66 ± 0.31.The bolometric luminosity is then estimated from the absolute band magnitude based on the relation in McLure & Dunlop (2004):   = −2.66(±0.05)log [ bol /W] + 79.36(±1.98).We note that the usage of the absolute -band magnitude as a proxy for absolute magnitude does not take into account possible effects of dust obscuration for the quasars.For most quasars this is expected to be small, but this may be more significant in red quasars.We discuss this further in Section 5.

Colours
We obtained the quasar ( −) colours for samples with  ≤ 2.5 from the observed SDSS − and −band magnitudes.We then used these colours to estimate the dust extinction for quasars across the entire redshift range, and hence the  ( −) colour excess, which we used to describe the colour-dependent properties of our quasar sample.To be more specific, we compared the colour against the redshiftdependent modal colour of ( − ) and ( − ) bands published in Table 1 of Hopkins et al. (2004) which therefore accounts for Kcorrection effects, and used the colour difference between the modal colour and observed ( − ) colour to obtain the reddening (relative to the average quasar) in the ( − ) bands.Following the work in Glikman et al. (2022) using WISE-2MASS selected quasars, we used an SMC 3 -like extinction law (Pei 1992) to convert the reddening into the relative colour excess Δ ( − ) of the quasar samples.Note that while these  ( −) values (resembling the difference from the average colour) will scatter around zero with apparently unphysical negative values indicating less reddening than the average quasar, our interest is in the tail towards the red end, where the choice to use Δ ( −) allows us to more easily compare across different redshift ranges.

PARAMETRIC MODEL OF OBSERVED RADIO FLUX DENSITY DISTRIBUTION
To construct a statistical model that translates the observed radio flux density distribution to a parametric likelihood, we adopt the framework presented in Roseboom & Best (2014) to allow the statistical likelihood of each individual source to be stacked together and modelled as an ensemble.This method makes full use of every source available and allows the overall probability to be constrained consistently.We will briefly summarize the Roseboom & Best (2014) approach below and then show how it can be adapted into our work.
The radio data consists of a set of sources with sky positions  and observed flux densities   .Instead of using a traditional approach to estimate the model parameters (e.g.M21), we instead use a bottomup approach by assuming the underlying flux densities to be   , in which case the pixel intensity can be modelled as: where   is a series of random statistical errors drawn from a Gaussian distribution  ( = 0,   ), and   represents the different uncertainties on the flux density measurements for each sample.If a certain model  can be constructed to predict the probability () 3 Small Magellanic Cloud of a source having underlying flux density , then the probability of observing a pixel intensity  for that source is: by convolving equation ( 1) over the entire parameter space for , folding in the Gaussian noise distribution.The stacked probability of a population with observed flux density distribution   and RMS flux density error  , can then be expressed as: Applying Bayes theorem, we can derive the likelihood of the model M based on the current distribution   : thus, the best-fitting parameters of model M can be found by maximizing the stacked likelihood in equation ( 3) under any given prior to the model (()).

Two-component model for radio emission
To construct the parametric expression of the probability of a source having true radio flux density  (()), we base our approach on the two-component model proposed in M21.Each quasar in the sample is assumed to have contributions to the radio flux density from both star-forming activity and AGN (jet) activity, while the total sample set is binned into grid cells on the   −  plane, each of which is modelled separately (see Figure 2).Here, we will first introduce the two model components individually, before explaining how we will combine them into a parametric expression of ().

SF Component
The host galaxies for radiative-mode AGNs are known to be massive and star-forming (e.g.Kauffmann et al. 2003;Best et al. 2005).M21 modelled the radio emission arising from the host galaxy starforming activities within each grid cell with a log-Gaussian probability distribution centred at a certain 150 MHz radio luminosity, log(  /[W Hz −1 ]) (tracing typical SFR Ψ within this population), and a scatter,   /dex.The conversion between SFR (Ψ) and radio luminosity (  ) is provided by Smith et al. (2021), being calibrated with LoTSS data and a Chabrier ( 2003) initial mass function: log with scatters within 0.3 dex in our typical SF galaxy radio luminosity range   ∼ 10 23−24 W Hz −1 .The scatter for SFR is propagated from the scatter in luminosity:  Ψ = 0.96  .Note that while there is some evidence for a mass-dependent SFR-radio luminosity relation (e.g.Algera et al. 2020;Delvecchio et al. 2021;Smith et al. 2021), we use the mass-independent form of the correlation in Smith et al. (2021) since we do not have robust measurements for the host galaxy stellar mass.Given the tight correlation between SFR and stellar mass at a fixed redshift (e.g.Speagle et al. 2014, and references therein) and that the estimated SFRs for the quasar sample are relatively high (see Section 4), we can infer relatively high stellar masses (> 10 10  ⊙ ) within a given bin.Therefore, the weak mass dependency in the Smith et al. (2021) relation is unlikely to have a significant impact on the result when compared to the intrinsic scatter in the  150 -SFR correlation.
Since the physical properties are expected to be similar for sources sharing a grid cell on the   −  plane, we expect a tight correlation between radio luminosity and star formation rate if either SFR correlates with gas accretion rate, or if host galaxies lie on the SF main sequence.This provides a broad justification for the log-Gaussian model we use here, and results from M21 further support the feasibility of our model.Note that values of   and   may vary between different grid squares: the normalisation of the star-forming main sequence has a strong dependence on redshift, due to more gas being available at the peak of the cosmic star formation rate density (e.g.Madau & Dickinson 2014); many studies have also found variance between typical host galaxy SFRs with AGN luminosities (e.g.Shao et al. 2010;Bonfield et al. 2011;Rosario et al. 2012;Dong & Wu 2016).
The probability distribution function (PDF) for the star formation component thus becomes (where L is the radio luminosity of the source):

AGN Component
Since the radio luminosity of radio-loud AGN is dominated by their jet component, we can extrapolate the luminosity function of radioloud AGNs to lower luminosities to reflect the distribution of jet luminosities for both radio-loud sources and radio-quiet sources, assuming they share the same jet mechanism, only with varying power efficiency.Note that while there are several other proposed scenarios for AGN-related radio emission in RQ AGNs, we only consider the AGN jets here for simplicity.Impacts from other mechanisms are negligible in the radio-bright end of our model, but might affect the radio-faint end; we will further discuss such impacts in Section 5. M21 used a single power-law distribution to model the emission from radio-loud AGNs (and thus the jet component emission), which results in a probability distribution function such that: where  is the power-law slope,  is the normalisation parameter which can be translated into the radio-loud fraction  defined in M21 (see the discussion in Section 3.1.3),and  min jet is the lower limit of the jet luminosity, required for normalisation (see below).
While the luminosity function for radio-loud AGNs is often modelled as a broken power-law distribution (e.g.Dunlop & Peacock 1990), this luminosity function traces the entire population of radio AGNs; therefore, quasars with different optical luminosities may occupy different parts of the function.As a result, if we are only sampling a part of the entire population in a given grid cell, the distribution may not follow the same pattern as the integrated population.M21 found that the lack of high-luminosity sources in a number of grid cells caused strong degeneracies between parameters, and the slope above the break luminosity could not be constrained.Here we adopt the M21 conclusion and use a single power-law function to model the jet component instead.We demonstrate later that this provides a very good fit to the data.
A practical issue with a power-law probability distribution function is that it grows monotonically when moving towards the faint end.
To tackle this problem, M21 set a lower luminosity limit ( min jet in Equation 7) for all quasars (below which the probability dropped to zero) such that the integral of the jet probability distribution comes to unity.They found that this lower luminosity limit is typically 10 19 ∼ 10 20 W Hz −1 , which is in line with the expectations from previous models (e.g.Mauch & Sadler 2007;Cattaneo & Best 2009;Sabater et al. 2019).

Total Simulated Radio Emission
While M21 randomly drew an  SF and  jet from the corresponding PDF and summed them in Monte Carlo simulations to get the overall luminosity distribution, in this work we need to derive a PDF for the individual quasar luminosity that we can convert to the probability function () in Equation 2. We take a different approach of rescaling and summing the SF and jet component PDFs (rather than summing luminosities drawn from the PDFs), as outlined below.
Under the two-component assumption, when the jet is very weak, the quasar luminosity will be dominated by SF emission; thus its distribution will follow the distribution function of the SF component.However, the growth of the jet distribution function at the faint end would instead let the jet component dominate at low flux densities if the two PDFs are simply summed, giving a clearly incorrect combination.We hereby make a simplifying assumption of letting  min jet =   in Equation 7; we then rescale and add the two probability distributions (Equation 6 and 7) together to form the probability density distribution for the total underlying radio luminosity ().
Under this assumption, when the jet luminosity falls below   , () is proportional to the PDF of the SF component: () ∝  SF (), since the SF component already dominates the entire luminosity distribution.On the other hand, when  jet ≫   the jet component becomes dominant in the distribution function and the overall PDF becomes the scaled version of  jet ().
To achieve this in practice, for any given jet distribution, we let  be the fraction of the sources that have  jet >   .The value of the AGN-dominated source fraction  is given by: Since the jet component dominates the luminosity distribution of these sources, we use () =  •  * jet () to describe their luminosity function, where  * jet () =  jet () for  jet >   and  * jet () = 0 for  jet <   .For the remaining (1 − ) fraction of the sources, the total luminosity is dominated by the SF contribution.Therefore, for these sources, we have () = (1 − ) SF .We thus define the luminosity function (or overall PDF) in our model as: For  ≤   we therefore have: and for  >   we have: .
We can thus convert between  and  using: Combining this scaling relation with Equation ( 8) gives the following conversion between the jet normalisation parameter  used in this work and the parameter  in M21: Having established the PDF for radio luminosities, (), the probability distribution of radio flux densities () in Equation (2) (assuming fixed redshift within each grid) becomes: where  () is the luminosity distance at redshift , and  () is the K-correction4 applied at .Note that here () is the probability of a particular quasar having an underlying radio flux density .
We refer the reader to Table 1 for a brief summary of the model parameters presented above.

Fitting with a Parametric Approach
With the expression of (), we can finally calculate the stacked likelihood of a certain parameter set {  ,   , , } under a radio flux distribution   and an RMS error   , using Equations ( 2) to (4).We then use a Monte Carlo Markov Chain (MCMC)-based algorithm (emcee; Foreman-Mackey et al. 2013) to determine the best-fit parameters based on the marginalised median of the probability density distribution defined in Equation (4).
Figure 3 gives an overview of our parametric model described above, for one well-populated grid cell.The left panel shows the flux density distribution in log space, which explicitly shows the twocomponent model discussed in Section 3.1; the right panel shows the model in linear space -reaching negative values -to better show the effects of the noise.The orange line and green line represent the radio flux density PDF of the log-Gaussian SF component (Equation 5) and the single power-law jet component (Equation 7) respectively.The pink dashed line represents the combined PDF of the two-component model (Equation 9), while the red dashed-dotted line shows the final PDF after convolving with the observational error (Equation 14).We have drawn from the quasar parent sample a grid cell containing 20,937 sources with optical magnitude −24 <   < −23 and redshift 0.8 <  < 1.2, and binned them in radio flux density (shown in blue crosses); the best-fit PDF agrees well with the distribution of the observed sources.There is a small mismatch at  just below   (the blue data points above the brown dashed-dotted line), perhaps caused by our method of combining the PDFs assuming no jet contribution here, while there might still be a weak amount from jet just below   , leading to an underestimation in the constructed PDF.Other well-populated grid cells share a similar pattern, despite being located at different places in the parameter space.
To test the limits of viability and robustness of our model within the parameter space explored, we have created mock catalogues with different numbers of sources (ranging from 250 to 10,000) sampled from the luminosity function in Equation 11 with different input parameters, and fitted our model with the mock catalogues.We then compare the fitted values of the model parameters against the input to test whether our fitting approach can retrieve the actual parameter values.Details on the validation results using mock data can be found in Appendix A1.
The test results show that the lower limit of   to which our approach can probe depends on the number of quasars in the bin -having more quasars allows us to probe fainter   .The lower limit of   corresponds to a flux density around 5 times the stacked noise level.Assuming an average noise level of LoTSS and a typical SFR of the SDSS quasars, our approach accurately retrieves input parameter values for grid cells containing at least 1,000 sources.To comply with the test results, we consider only grid cells with more than 1,000 sources (10,000 sources for analysis regarding other physical parameters, since we subdivide our samples into 10 bins) in the following analysis; this is indicated by the colour-coding of grid squares in Figure 2.
The mock test results indicate a degeneracy between   and  in our fitting; this is most likely due to the relative scarcity of sources in the radio-loud regime, which makes the best-fit result of  depend more on the sample distribution in the region where both the SF and jet make a significant contribution to radio emission.On the other hand, assuming a single power-law distribution in the radioloud quasars provides us a simple way to estimate  without running MCMC fits.As the host galaxy SF contribution becomes insignificant in the radio-bright tail of the distribution, the power-law slope of the distribution at high radio luminosities provides a good estimation for  in our model.We can then use the fitted values of different power-law slopes across different bins as strong informed priors to the  values in our model, and thus resolve the degeneracy between   and  in the radio-bright end.
The results of the analysis on the  prior are provided in Appendix A2, where Figure A3 shows the radio-bright end of the quasar radio luminosity function for all grid bins investigated in this work (blue line).The lower luminosity limits for analysis of the radiobright luminosity functions are selected using the following criteria: (i) not lower than the 10 radio flux density limit and (ii) at least 1 dex above the estimated   to ensure negligible SF contribution, so that the distributions only show a single power-law feature in the log  − log  space (values listed in Table A1).We used a single power-law model (() ∝  − ) to fit the radio-bright distributions, as shown by the red dotted lines in each grid cell.The single power- Note that the actual probability density of our Gaussian flux density error peaks close to zero in linear units of flux density, with the difference between the two panels being that the bin size in the left panel (with logarithmic bins) becomes smaller towards the radio-faint end, hence the number count within each bin will drop with the bin sizes.The red dashed-dotted line shows the PDF after convolving with the observational error.We have drawn sources from a grid cell in our catalogue (−24 <   < −23, 0.8 <  < 1.2) and binned them in radio flux density to get the observed number count within each bin (shown in blue crosses, with the uncertainty shown in thick blue lines along the y-axis).The observed distribution agrees well with the proposed model.
law model provides a good fit across the entire parameter space, which further justifies our model for the jet components, and therefore supports the ubiquity of jets in our quasar samples.The fitted value of  is shown in Figure A3 for every grid cell, while Figure 4 gives an overview of the  values across some of the most populated grids.These values show little dependence on the redshift or optical magnitude, which agrees with the conclusion in M21 using the chi-square fitting.
As a result, within each grid, the informed prior of  is defined as a Gaussian distribution: where  0 is fixed at 1.5 at all grid bins due to little evolution with redshift or optical magnitude, and   is set to 0.05 to give  a tight range.
For the rest of the parameters, we have adopted simple box priors (19 < log   < 30, 0.05 <   < 0.5, −5 < log  < −0.25) to maximise the information obtained from the MCMC fitting while ruling out the unrealistic results.The final expression for the likelihood of the model is therefore given by Equation 4assuming () = ().

DEPENDENCIES WITH 𝑀 𝑖 AND z
In this section, we present the results related to the best-fit parameters from our proposed model.To better compare our results with the previous ones from M21, we use the scaling relation in equation 13 to convert the jet power normalisation parameter  obtained in this work to the parameter  presented in M21.In Section 3.1 we also proposed a more physically motivated parameter -the jet probability scaling factor () -defined as the fraction of quasars with radio emission brighter than the average host galaxy contribution, on which our analysis will focus mainly.
Figure 5 shows the variation in our fitted parameters across different subsamples binned by their absolute -band magnitudes and redshifts, while we present the detailed values of best-fit parameters in Table A1.Compared to previous results in M21, our method extended the investigated redshift range to  ≈ 3.8 and the optical luminosity range up to a magnitude deeper, thanks to the wider sky area coverage of LoTSS DR2 and the improved parametric fitting method, which requires fewer samples for a good fit.The fitted values for the star formation rate of the host galaxies (Ψ) exhibit a much higher precision when compared to the similar results using chi-square-based approach in M21 -another major improvement made possible by our Bayesian fitting method.However, the values of Ψ in our work deviate significantly more from the best-fit values in M21 than the errors, in both faint-  and high- sources.This is believed to be largely due to the differences between our model and the M21 model.Firstly, we have calculated a higher prior for  ( = 1.5) compared to the M21 model ( = 1.4) due to the difference in our fitting strategy (see M21 for details), which leads to an increased fraction of jetdominated sources in the radio-intermediate regime.As a result, the scatters of the SF component () that we derive is notably smaller compared to M21 since fewer SF-dominated sources are required in the radio-intermediate quasar population, which would therefore lead to higher fitted values of peak SFR (Ψ) to keep a consistent jet normalisation in the radio-loud regime. 5Secondly, the way that we combine the PDFs for the SF and AGN components is subtly different, which predominantly influences sources around   ; where   corresponds to a flux density below the noise limit, this can impact the results.Thirdly, we convolved our theoretical PDF with different radio flux uncertainties of each individual sources (see Equation 2), while the M21 model used the average radio flux uncertainty for the full sample; although our method maps the actual physical scenario more accurately, it can also lead to deviations from the M21 result.However, whilst all of these factors affect our comparison with the of M21 values, variations between quasar sub-populations within our model are very robustly determined.
The absolute -band magnitude (  ) is a good tracer of the bolometric luminosity, and is hence associated with the BH accretion rate of a quasar (see Section 2.4.2).Note that although there is a degeneracy between BH mass and BH accretion rate in their correlations with the bolometric luminosity, in this paper we will focus on the BH accretion rate side and leave the investigation of possible correlations with BH mass for future work.In this study, we have observed a positive correlation with   / bol in both host galaxy SFR and AGN activity level (characterised by ).By assuming SFR ∝   bol under fixed redshifts, we have  = 0.26 ± 0.02, which is in agreement with the value found in Bonfield et al. (2011) where they used far-IR Herschel measurements instead as a SF tracer.We have also compared our result with the  SF,1.5GHz -  correlation presented in White et al. (2017), where they obtained 1.5 GHz radio measurements of 70 radio-quiet quasars at  ∼ 1 using the Karl G. Jansky Very Large Array (JVLA).Their best-fit relationship, with  SF,1.5GHz converted to a 150 MHz radio SFR using the radio spectral slope defined previously, is plotted as the grey dashed line in Figure 5, and is in agreement with our SFR- bol correlation at 0.8 <  < 1.2 (albeit under a different assumed best-fit function).
Our results also agree with the results in M21 when comparing the fitted values of log  .Adopting their definition of  ∝   bol , we find  = 0.67 ± 0.03, while M21 gives  ∼ 0.65.This positive correlation is in line with our current knowledge on galaxy evolution and AGN properties (e.g.Jiang et al. 2007), as the gas ensemble within the halo fuels both host galaxy star-forming activity (Kennicutt 1998) and the accretion activity around black holes.Therefore, higher quasar optical luminosities would indicate a more abundant gas reservoir in the galaxy haloes, which is tied to more intense star formation in host galaxies (e.g.Koss et al. 2021, and references therein).Meanwhile, we emphasize that despite the increase in the fraction of high jet power sources at higher optical luminosities, quasars spanning the full range of the distribution of jet powers are seen across the whole optical luminosity space -the increase in  is the result of scaling the full power-law distribution to higher  values in the overall probability distribution.
At fixed   , our modelled host galaxy SFR also shows a strong increase with redshift out to  = 2.5, which is in line with previous studies on quasar host galaxies (e.g.Bonfield et al. 2011;Rosario et al. 2012;Mullaney et al. 2012;Harrison et al. 2012).Note that while the host galaxy SFR increases with both quasar bolometric luminosity and redshift, the effect of Malmquist bias is minimal in our fitting process.This is because our fitting approach uses quasar samples within a small range of bolometric luminosity and redshift, therefore our dependencies on luminosity are derived for bins at fixed redshifts and vice versa.While the M21 correlation showed hints of a turnover beyond  = 2.5, there were not enough redshift coverage and too few   bins to properly characterise it; alternatively, our result shows that despite a flattening in the correlation, the SFR continues to increase with redshift out to  ∼ 4. The reason behind such an increase might be connected to the high prevalence of merger-induced starburst activities at earlier cosmic times, which is often linked to the triggering mechanism of powerful quasars (e.g.Sanders et al. 1988;Hopkins et al. 2006).Lamastra et al. (2013) predicted that the percentage of burst-dominated star forming galaxies increases with redshift, from ≤ 0.5% at  ∼ 0.1 to ∼ 20% at  ∼ 5, due to mergerinduced starbursts.The SFRs of quasar host galaxies that we derive at higher redshifts lie beyond the predicted values of the star-forming main sequence, hence showing an association with starburst galaxies, and thus an enhancement at higher redshifts relative to the global cosmic SFR density.Duncan et al. (2019) further posed observational constraints on the merger histories up to  ∼ 6, indicating a higher merger fraction at  > 3, which is consistent with the cosmic epochs when our SFH of AGN host galaxies differentiates from the cosmic SFR density in star-forming galaxies.
On the other hand, the jet activity level shows little sign of evolution with redshift (see middle panel of Figure 5).There is a small upturn in the lowest redshift bins in the  plots, but this is a result of integrating the power-law jet distribution to a fainter limit at low redshifts due to the lower value of   , rather than a true evolutionary trend: the lack of a trend of either  (Figure 8) or  (Figure 4 and A3) with redshift implies that the distribution of jet powers of the quasars is largely unchanged.While the AGN activity and host galaxy star formation are both tied to the gas reservoir within the galactic halo, the lack of redshift dependency suggests that the former is more likely to be    2021) (  ), jet probability scaling factor from this work (), and their variation with absolute -band magnitude (  , top row) and redshift (, bottom row).The solid lines show the result of this work, and we compare our result with previous results in M21 (dashed lines).In line with M21, we observe a significant dependence of both SFR (Ψ) and jet activity (  or ) with   (a proxy for bolometric luminosity and accretion rate), while only the SFR evolves strongly as a function of redshift for fixed   (see the text for further discussion).Our SFR-  correlation at 0.8 <  < 1.2 also shows good agreement with the best-fit relationship in White et al. (2017) where they studied the radio emission for 70 radio-quiet quasars at  ∼ 1 (grey line).dictated by local activities around the galactic core rather than the time evolution of the gas reservoir.
It is not possible to provide direct evidence to the mechanism that powers the AGN radio emission within the scope of this work; however, if we assume this emission arise from radio jets, since the level of jet activity links to the jet powering efficiency of the investigated quasar population, we can still speculate on the physical processes behind AGN jet production upon knowing the evolution pattern of parameter .Some studies (e.g.Blandford & Payne 1982;Woo & Urry 2002) back up the claim that black hole spin is related to the efficiency of jet production (Wilson & Colbert 1995;van Velzen & Falcke 2013), with more rapidly-spinning black holes producing more powerful radio jets.In this picture a RL/RQ dichotomy exists, since the black hole spins tend to be either high or low depending on their accretion history.On the other hand, the excellent fits to the data for our model, in which there is no dichotomy between RL and RQ quasars but rather a continuous distribution of jet powers from very high (∼ 10 30 W Hz −1 ) down to sub-dominant compared to the radio emission from SF, is inconsistent with that theory: the wide range of jet power distribution requires a corresponding wide range of black hole spin parameters, which is hard to find in simulations (e.g.Volonteri et al. 2013).Instead, our result tends to favour the alternate theory presented in Sikora & Begelman (2013) that the jet power is controlled by the magnetic flux threading a spinning black hole.Their model is able to create radio jets spanning a wide power range from a variety of magnetic black hole accretion flows, including "magnetically-choked accretion flows" for high-luminosity jets and magnetic field fluctuation in the coronae above thinner accretion disks for low and intermediate power jets.They therefore predict an increase in jet power with increasing black hole accretion rate; both predictions are in line with our results.Note that the discussions above assume a moderate range of black hole masses, while more massive black holes can also produce a stronger magnetic field that leads to increased jet power.We aim to break the degeneracy between black hole accretion rate and black hole mass in our future work.
As for the two remaining parameters, we are able to obtain fitting results on   through the Bayesian scheme and  through examining the radio-loud sources.Our results for   showed little sign of variation with either   or  (see Table A1), while a similar situation with  has already been discussed earlier and is shown in Figure 4.Both results agree with the analysis in M21.

LINK BETWEEN DUST ATTENUATION AND RADIO EMISSION
Thanks to the LoTSS DR2 data, we now have additional sample statistics available within the most populated regions of the parameter space.With these extra measurements, we can use the method proposed in this work to study the variance of host galaxy SFR and jet power distribution separately with any measured physical quantity for the quasars, by further splitting the parameter space in our fitting.
In particular, in order to definitively identify the origin of the extra radio emission associated with red QSOs discussed in Section 1, we split our samples by optical colour.
In this work, instead of defining a subsample of red quasars based on a redshift-dependent percentile cut-off in observed  −  colour (Klindt et al. 2019) or Δ ( − ) colour excess (deviation from the average  ( − ) value; Glikman et al. 2022), we investigated the colour-related variations across the full range of colour excess; we then discuss our results in the context of previous similar studies on the properties of red/blue quasars.This allows us to speculate into the continuous evolutionary trend of the relevant parameters, rather than only comparing the differences between two populations.As discussed in Section 2.4, we focus on  ( − ) instead of observed  −  colour, as we assume that the redder colour in rQSOs is due to dust extinction only, which has been verified in Fawcett et al. (2022).By using  (−), we can therefore better compare between different bins without trends between different   −  bins introducing any biases.
Figure 6 shows the cumulative distribution of the redshiftcorrected Δ (−) (i.e. the difference of  (−) from the average value) in the most densely populated redshift bin ( = 1.4) for the parameter space explored in this section.The grey dotted line and dashed line compare the selection criteria for rQSOs in Klindt et al. (2019, upper 10% of colour distribution) to the criteria in Glikman et al. (2022, Δ ( − ) > 0.25) .Both selections broadly agree with each other across the entire luminosity range, suggesting our sample shows a similar intrinsic colour distribution with the previous works.
To define a continuous trend of parameter variation, we binned each grid in the   −  space into 10 sub-grids based on the colour percentile.Note that we only picked grids with more than 10,000 sources in this process so that each sub-grid contains at least 1,000 sources, which is consistent with the required minimum source number to obtain good fits stated in Appendix A1.
We characterise the samples in each sub-grid by their radio luminosity ( 150 ), -band magnitude (  ) and redshift ().The distributions of these physical properties across different colour sub-   Such uniformity enables comparison of fitted values across different colour percentile samples.The distributions of  150 , on the other hand, show an extended radio-bright wing in the reddest samples, while the locations of the star-forming peak (as defined in Figure 3) remain unchanged across different colour percentiles.A detailed analysis is presented in Figure 8.Such distributions are seen in all   - grids used for quasar colour analysis.
grids are presented in Figure 7, using the quasar samples from the −24 <   < −23, 0.8 <  < 1.2 grid.The similarity of both   and  distributions across the colour sub-grids indicates there are no systematic biases when studying the evolution of SF and AGN activities with quasar colour percentiles.On the other hand, the distributions of  150 differ in the reddest populations (top, i.e. 30% colour percentile) from the rest, in terms of an extended wing on the radio-bright end.This trend is seen across all grids inspected in our study (see Figure A4).The difference in the distribution of radio powers can thus be accurately modelled by our Bayesian approach.Figure 8 shows the relative variation of host galaxy SF activity (characterised by the mean radio SFR Ψ) and AGN activity (characterised by the jet normalisation parameter  ) with the  ( − ) colour percentile.Here, the term 'relative variation' is defined as the ratio between the fitted value within each colour sub-sample (Ψ,  ) and the fitted value using the entire population within the parent grid cell prior to splitting by colour (Ψ 0 ,  0 ).
Since our model has the ability to separate the contribution from SF and AGN in the radio flux density distribution, we can demonstrate that the AGN activity is the main driver behind the reddening of quasars as it strongly correlates with the optical colour (right panel); quasars with the reddest colour also have the highest jet fraction.The host galaxy SF, on the other hand, shows a much weaker correlation with optical colour (left panel).Comparing the correlations across different bins (grey lines), we see little sign of redshift evolution, which suggests that rQSOs are not associated with the star forming activity of their host galaxies, but more likely with the AGN evolutionary stages.The difference in the redness of quasars is therefore related to the dust content surrounding the AGN rather than the dustiness of the host galaxy.Note that our colour-split samples only cover three   −  grids beyond the cosmic star formation peak at  ∼ 2; therefore deeper radio observations are needed to draw a complete census across cosmic times.Using additional multiwavelength data in the smaller LoTSS Deep Field DR1 (Tasse et al. 2021;Sabater et al. 2021) coverage instead, Calistro Rivera et al. (2023) developed an independent method to separate AGN and SF components in quasar SEDs and came to a similar conclusion that radio emission in rQSOs comes almost exclusively from AGN.Furthermore, as the host galaxy SF also correlates with the absolute i-band magnitude (Figure 5), this weak trend between SF and optical colour may as well be the side effect of extinction.
Extinction will result in the measured -band magnitude being fainter than it should be due to the dust, especially for quasars in the redder bins.Since correction based on the extinction law will again increase the measurement of   , the corrected values for  bol will increase in higher bins; as a result, the expected SFR will also be higher due to the correlation between Ψ and   in Figure 5, which may produce an artificial trend shown by the orange dashed line in Figure 8.To test this, we deduced Δ log Ψ = Δ  from the fitted value, where  is the slope inferred from the SFR −   correlation in Section 4, and Δ  is the correction in the measured -band magnitude based on the extinction law and the given  ( − ).The weak correlation observed between fitted Ψ and Δ (−) is removed after the correction (solid line), confirming that this correlation is merely an artefact of the dust reddening in the photometric bands.This result shows that the SF contribution to the radio emission remains unchanged between red and blue QSOs of the same bolometric luminosity.We have taken a similar approach to the log  − Δ ( −) relation using a slope inferred from the  −   correlation in Section 4, and the corrected relation is shown as an orange solid line in the right panel.The corrected result still shows a strong positive correlation between AGN activity and optical reddening, and we therefore rule out the possibility that this correlation arises from the side effect of dust extinction.
Recently, Fawcett et al. (2022) argued that radio excess in rQSOs occurs primarily in quasars with lower radio powers, likely due to the interaction between weak/intermediate jets and the opaque interstellar medium/circumnuclear environment in rQSOs.Our model allows us to speculate about the details of this process by exploring whether  is changing.Figure 9 shows the variation in fitted  values across different colour percentiles, using the fitting method described in Section A2.Our results reveal a universal increase in fitted  values for rQSOs, which indicates that more radio-faint quasars have been enhanced than radio-bright quasars in the reddest bins, thus causing the rise in the slope of flux density distribution in the radio-bright end.We therefore confirm previous literature results (Klindt et al. 2019;Rosario et al. 2020;Fawcett et al. 2020Fawcett et al. , 2022;;Calistro Rivera et al. 2023), and conclude that the change in evolutionary phase that caused the reddening of rQSOs takes place mostly in radio-quiet and radio-intermediate quasars.
It is also worth noting that our results do not rule out the contribution of wind shocks in the radio enhancement of red quasars, as investigated in Petley et al. (in prep) based on a sample of broad absorption line quasars (BALQSOs).While our model can effectively rule out the host galaxy contribution in the radio enhancement, the AGN component in our model consists of any radio emission generated in BH activity that does not follow a Gaussian distribution centred at the average radio SFR, including jet and wind activities.The steepening of  could also be indicative of faint AGN radio emission other than weak jet (e.g.wind shock) being enhanced in the red quasars.If the radio luminosity function of that other component becomes reasonably well understood in the future, the methodology developed in this work can easily be expanded to explicitly include a third component (e.g.merger history, wind shocks, etc.).

CONCLUSIONS
In this work, we adopted the assumption made in M21 and developed a fully Bayesian two-component model (see Figure 3) that statistically disentangles the host galaxy SF and AGN jet contributions to the quasar radio emission, using the quasar radio flux density distribution as input.We assumed a log-Gaussian distribution for the SF flux densities with two parameters   (mean radio SF luminosity) and   (Gaussian scatter of radio SF luminosity), and a single powerlaw distribution for the jet flux densities with an additional two parameters  (or ; the jet power normalisation factor) and  (powerlaw slope).To investigate the variation of best-fit values of our model parameters with a number of factors including bolometric luminosity ( bol ), redshift () and optical reddening (characterised by the colour excess  ( −), we binned our samples into grid bins based on their location in the parameter space of   ,  and colour, before fitting the quasar radio flux density distributions within each grid bin against our model.
After various analyses, we have reached the following conclusions: (i) Our model confirms the main results of M21 but with a higher accuracy and a wider coverage of the parameter space, probing sources at a magnitude deeper and with a redshift out to  ∼ 4. We obtained good fits using our two-component model across the entire parameter space, which argues against an RL/RQ dichotomy.
(ii) The host galaxy SFR (Ψ) and jet activity () share a positive correlation with   , which can be associated with the black hole accretion rate.While the host galaxy SFR also positively correlates with redshift, the AGN jet activity experiences little evolution across cosmic time (Figure 5).Our result supports the speculation in M21 that the overall gas fraction and dynamical time within the halo will affect both star formation and black hole accretion; the fact the host galaxy SFR continues to grow when  > 2.5 goes beyond the M21 result, as they lacked sufficient source to probe beyond the redshift limit, and is different from the cosmic star formation rate density measured in star-forming galaxies (e.g.Madau & Dickinson 2014).This trend can be explained by the high prevalence of galaxy mergers in earlier cosmic times, which are likely to be responsible for the triggering mechanism of both powerful quasars and starburst activities.
(iii) The host galaxy SFR is found to have little impact on the excessive radio emission in red quasars (rQSOs), while the jet activity (SF component)   f (AGN component) plays a significant role as it is found to be increasing at all colours redder than the average quasar colour (Figure 8).This provides direct evidence to the evolutionary model that explains the radio enhancement in rQSOs, and rules out the possibilities of SF influence.The slope of the radio-bright end of the flux density distribution shows signs of positive correlation with the quasar colour (Figure 9), which indicates rQSOs with weak or intermediate jet activities are more likely to experience a radio enhancement.
This work shows the flexibility of our model that will enable us to investigate the role of host galaxy and AGN separately in terms of the correlation between quasar radio emission and various other physical parameters, including the black hole mass, clustering environment, AGN winds and outflows, etc.We plan to extend our method to examine the correlations stated above using the newest LoTSS and SDSS survey data.Future spectroscopic surveys including WEAVE-LOFAR (Smith et al. 2016) will provide more robust measurements to the continuum and emission line properties of LoTSS-detected quasars, thus giving a deeper insight into the details of the accretion processes, outflow, circumnuclear environments and dust composition of quasars.Other surveys such as DESI (DESI Collaboration et al. 2016) and WEAVE-QSO (Pieri et al. 2016) will also further expand the samples of optically selected quasars.Combining the techniques and new observational data, we are hopeful of understanding the sources of radio emission within quasars and their relative contributions.Note that here the relative difference is defined by Equation A1.We adopted a simple box prior to   and a sharp Gaussian prior to  centred on  = 1.5, as described in Section 3.2.We thus investigate the dependency on the input values of the remaining two parameters:   and  (fraction of jet-dominated quasars).For   , the bestfit result reaches maximum difference in the bottom left corner where both SF and jet components are weak, therefore making it hard to separate the two components.As the jet component grows stronger (or the SF luminosity increases), the SF contribution becomes easier to separate from the entire distribution.The pattern remains valid for all sample numbers, while the lower limit of the parameter space with good fits corresponds to 5 times the stacked noise level.Lower source counts lead to higher stacked noise levels, thus limiting the suitable parameter space.

𝛿
Figure A1 and A2 shows the fitting quality (characterised by  log   and  log  ) of   and  under different model parameters (characterised by the log   -log  parameter space) and different number of fitted sources.We targeted these two parameters for investigation because they are pivotal to the model construction (see Section 3.1), and have important roles in explaining the physical process of quasars (  translates to the mean luminosity of host galaxy star formation and  to the fraction of jet-dominated quasars).For the rest of the parameters, we adopted the priors described in Section 3.2 for consistency.
Both figures show a minor trend between the fitting quality and the number of sources included.As the number of sources decreased from 10,000 to 1000, the lower limit of intrinsic   required for a good fit increased, in a trend corresponding to the increase of stacked noise level of simulated sources.The poorest fitting quality of   occurs at the faint end of both SF and jet emission, and gradually improves with increasing intrinsic SF and jet luminositynote that there is no significant difference in the fitting quality once the SF luminosity reaches 10 times the stacked noise level, although higher jet fraction leads to slightly worse results in   under the same intrinsic SF luminosity.The fitting quality of  , on the other hand, correlates more with the intrinsic jet fraction compared to the situation in   .While the limit of intrinsic SF is still governed by the stacked noise level, an additional limit to intrinsic jet power is found in the parameter space, as the fitting quality for  worsened in the high-SF, high-jet power regime.This is mostly due to the confusion between the bright end of the SF distribution and the faint end of the jet distribution, when the increased mean SF luminosity and jet power normalisation moved these two populations into the same space in the radio flux density distribution.The fitting quality within the parameter space for good fits also slightly deteriorated for  , when compared to   .When the source count falls down to 250, the fitting quality continues to deteriorate, so that we cannot determine a clear boundary within our parameter space for guaranteed good fits.Therefore, we chose 1,000 sources as the minimum number of quasars required for a good fit.We selected the target grid bins for analysis in the   −  plane based on this criterion.For colour-dependent studies, the minimum number is set to 10,000 sources since we are binning each grid bin into 10 sub-grids based on the colour percentile, so that each sub-grid contains at least 1,000 sources.

A2 Priors
Following the discussion in Section 3.2, we present some detailed results used to determine the prior on  in order to eliminate the degeneracy between  and   / Ψ in our fitting.Figure A3 shows the selection of the bright-end quasar radio luminosity function fits used to determine the  prior.The determined values of , together with their uncertainties, are given in the final two columns of Table A1.A1.The variation of best-fit parameter values (for  we show the best-fit values obtained in Section A2) and uncertainties within the   - plane for the full sample, together with the lower luminosity cut-off for determining the radio-loud luminosity function defined in Appendix A2 ( min ).We also included the black hole accretion rate (log  BH ) for reference, which is calculated from Equation 2

A3 Radio flux density distribution for red quasars
In Figure A4, we present the actual radio luminosity distributions of different colour percentile quasars across all   −  grids used for the study of colour dependence in Section 5 (note that sources with negative flux densities are excluded when plotting with radio luminosities).The radio luminosities are binned into the same bins defined for the top panel in Figure 7, and are plotted in the colour-code defined for different colour percentiles as shown in the bottom-right legend.The red dashed line shows the radio luminosity distribution for the entire population within the grid.Here we would like to remind the readers that a fair amount of sources lie under the 2 flux density uncertainty limit; therefore the distributions presented in Figure A4 cannot be used to deduce the exact levels of SF and AGN activity.Readers should still refer to Figure 8 for the model-fitted values of   and  in different coloursplit samples.Across all   −  grids, the radio flux density distributions for rQSOs (defined as the top 10% in colour excess) show an extended radiobright wing while sharing similar star-forming peaks when compared to the distributions for the entire population.This is in accordance with the upper panel in Figure 7    show an extended wing on the radio-bright end while sharing similar star-formation peaks when compared to the entire population.Note that the shaded areas reflect sources with flux densities below the 2-sigma LoTSS DR2 uncertainty limit and therefore cannot be used to visually determine the exact positions of the star-formation peak or the relative jet power unless fitted with our two-component model.

Figure 2 .
Figure 2. Distribution of sources used in this work in the absolute i-band magnitude (  )-redshift (z) plane.The samples are later separated into subsets in order to study the dependency of best-fit model parameter on quasar optical properties, using the grid lines in this plot.The grid with solid lines shows the parameter space explored in Macfarlane et al. (2021) with LoTSS DR1 data, while the dashed grids show the new parameter space explored in this work.Grids with orange lines indicate the parameter space used to investigate the colour dependency in this work, i.e. grid cells with more than 10,000 quasars.Thanks to the new data from LoTSS DR2, we are able to investigate sources with fainter radio emissions and higher redshifts.

Figure 3 .
Figure3.Separate and combined probability distribution functions (PDF) of radio flux densities (), illustrating different components in our model.The left panel shows the flux density distribution on a logarithmic scale which offers the clearest representation of the model, while the right panel shows the flux density distribution on a linear scale, which is important to show the negative values (due to noise) and is the form used in our fitting algorithm.In each panel, the orange line and green line represent the SF component from the host galaxy activity (log Gaussian) and the jet component from the AGN activity (single power-law) respectively.The pink dotted line represents the combined PDF of the two-component radio flux density distribution model.The vertical grey dashed line indicates the lower radio luminosity limit that we used to fit  value in this   - grid (see Section 3.2 for details).The grey solid line compares the Gaussian error introduced by flux measurement (using average flux error in the extraction from DR2 mosaics) with the modelled distribution and the actual distribution.Note that the actual probability density of our Gaussian flux density error peaks close to zero in linear units of flux density, with the difference between the two panels being that the bin size in the left panel (with logarithmic bins) becomes smaller towards the radio-faint end, hence the number count within each bin will drop with the bin sizes.The red dashed-dotted line shows the PDF after convolving with the observational error.We have drawn sources from a grid cell in our catalogue (−24 <   < −23, 0.8 <  < 1.2) and binned them in radio flux density to get the observed number count within each bin (shown in blue crosses, with the uncertainty shown in thick blue lines along the y-axis).The observed distribution agrees well with the proposed model.

Figure 4 .
Figure 4.The  value inferred from the power-law slope at the bright end of the radio luminosity distribution, plotted for the most populated grids.The values of  show little change with redshift (left panel) or absolute -band magnitude (right panel), agreeing with the conclusion in M21.Therefore, we fixed the value of  to  ≈ 1.5 when fitting the general sample (not binned with colour), using a tight Gaussian prior.

Figure 5 .
Figure5.The best-fit values of the mean host galaxy SFR (Ψ, as determined from   using Eq.5), jet power normalisation fromMacfarlane et al. (2021) (  ), jet probability scaling factor from this work (), and their variation with absolute -band magnitude (  , top row) and redshift (, bottom row).The solid lines show the result of this work, and we compare our result with previous results in M21 (dashed lines).In line with M21, we observe a significant dependence of both SFR (Ψ) and jet activity (  or ) with   (a proxy for bolometric luminosity and accretion rate), while only the SFR evolves strongly as a function of redshift for fixed   (see the text for further discussion).Our SFR-  correlation at 0.8 <  < 1.2 also shows good agreement with the best-fit relationship inWhite et al. (2017) where they studied the radio emission for 70 radio-quiet quasars at  ∼ 1 (grey line).

Figure 6 .
Figure 6.The cumulative distribution of Δ ( −  ) (deviation in  ( −  ) from the average value) in the most populated redshift range (1.2 <  < 1.6).The   and  values in the legend are the mid-points of the corresponding   −  bins.The observed  −  colour has been corrected to the dustindependent  (− ) colour excess using the correction method presented in Section 2.4.The grey dotted lines show the definition of red QSOs in Glikman et al. (2022) (Δ ( −  ) > 0.25) and Klindt et al. (2019) (upper 10% of the colour distribution).The two criteria intersect close to the cumulative colour distribution curve, suggesting that our sample shows a similar intrinsic colour distribution with both works.

Figure 7 .
Figure7.Distributions of sample properties within a representative   - grid, characterised by radio luminosity ( 150 ), -band magnitude (  ) and redshift (), while separated by 10-percentile colour sub-grids.The distributions of   and  show no significant variation across the sub-grids, indicating a uniform sample property among different colour percentiles.Such uniformity enables comparison of fitted values across different colour percentile samples.The distributions of  150 , on the other hand, show an extended radio-bright wing in the reddest samples, while the locations of the star-forming peak (as defined in Figure3) remain unchanged across different colour percentiles.A detailed analysis is presented in Figure8.Such distributions are seen in all   - grids used for quasar colour analysis.

Figure 8 .Figure 9 .
Figure8.The relative variation of model parameters (Ψ/Ψ 0 ,  /  0 ) plotted against the deviation in quasar colour excess from the average value (Δ ( −  )), where Ψ,  are the best-fit parameter values within each colour percentile bin, and Ψ 0 ,  0 are the best-fit value using the entire   - grid.Within each   - grid, the quasars are binned into 10 bins based on the number percentile of  ( −  ) colour.Grey lines indicate the trend in each grid cell, whereas the orange dashed line shows the average variance of all investigated grid cells.The left panel shows the variation in the mean host galaxy SFR (Ψ), which appears to be a slight increase in Ψ for redder colours (orange dashed line) but shows no variation with Δ (  −  ) once the effect of dust extinction biasing the   values has been removed (orange solid line; see text).The right panel displays the variation in radio jet normalisation (  ), while the dashed and solid lines represent the variation before and after correcting for dust extinction effect, respectively.As the quasar colour becomes redder, we see a significant increase in the relative jet intensity, showing that the cause of the radio excess in red quasars is indeed due to increased AGN (jet) activity rather than the host galaxy SF.
Figure A1.Relative difference between input and best-fit values of   (mean host galaxy radio SF luminosity).Note that here the relative difference is defined by Equation A1.We adopted a simple box prior to   and a sharp Gaussian prior to  centred on  = 1.5, as described in Section 3.2.We thus investigate the dependency on the input values of the remaining two parameters:   and  (fraction of jet-dominated quasars).For   , the bestfit result reaches maximum difference in the bottom left corner where both SF and jet components are weak, therefore making it hard to separate the two components.As the jet component grows stronger (or the SF luminosity increases), the SF contribution becomes easier to separate from the entire distribution.The pattern remains valid for all sample numbers, while the lower limit of the parameter space with good fits corresponds to 5 times the stacked noise level.Lower source counts lead to higher stacked noise levels, thus limiting the suitable parameter space.

Figure A2 .
Figure A2.Relative difference between input and best-fit values of  (fraction of jet-dominated quasars).While the pattern and the limit for suitable parameter space remain similar to FigureA1, the best-fit values of  is more sensitive to input, and the uncertainties on the derived values for  increase in the regions of parameter space where the input jet-dominated quasar fraction is high or the SFR is low.
and the best-fit parameter values presented in Figure8

Figure A4 .
Figure A4.The radio luminosity distributions separated by 10-percentile colour sub-grids (solid lines) compared with the radio luminosity distribution of the population (red dashed lines), plotted for every   −  grid investigated in Section 5. Higher percentiles indicate higher values in Δ ( −  ) colour excess, hence redder quasars compared to the entire population.Across all   −  grids, the distributions for rQSOs (top 10% of the Δ ( −  ) colour excess) show an extended wing on the radio-bright end while sharing similar star-formation peaks when compared to the entire population.Note that the shaded areas reflect sources with flux densities below the 2-sigma LoTSS DR2 uncertainty limit and therefore cannot be used to visually determine the exact positions of the star-formation peak or the relative jet power unless fitted with our two-component model.
, which covers 424 deg 2 .When compared to the Faint Images of the Radio Sky at Twenty Centimetres survey (FIRST; Becker et al.
LoTSS DR2 contains ∼4,400,000 radio sources -a ten times increase from the DR1 catalogue, which is ideal for multi-wavelength studies.

Table 1 .
A list of our model parameters and their definitions.All luminosities refer to the extracted radio luminosities from LoTSS catalogue.  /[W Hz −1 ] Average SF luminosity of the population   /dex Gaussian scatter in the SF luminosity distribution  Jet power-law slope as in  jet () =  −  Jet normalisation factor as in  jet () =  −  Fraction of sources with  jet >   (AGN-dominated) It is also worth noting the relationships between the jet normalisation parameter  or the AGN-dominated fraction  used in this work and the similar parameter  defined by M21.The parameter  in M21 is defined as the fraction of sources with jet luminosity brighter than   , where log(  /W Hz −1 ) = 26.Therefore, we have  = ∫ ∞    − .Our definition of , on the other hand, gives in M21.