The intrinsic X-ray luminosity distribution of an optically-selected SDSS quasar population

In active galactic nuclei, the relationship between UV and X-ray luminosity is well studied (often characterised by $\alpha_\text{ox}$) but often with heterogeneous samples. We have parametrized the intrinsic distribution of X-ray luminosity, $L_\text{X}$, for the optically-selected sample of SDSS quasars in the Stripe 82 and XXL fields across redshifts 0.5-3.5. We make use of the available XMM observations and a custom pipeline to produce Bayesian sensitivity curves that are used to derive the intrinsic X-ray distribution in a hierarchical Bayesian framework. We find that the X-ray luminosity distribution is well described by a Gaussian function in ${\log_{10}}L_\text{X}$ space with a mean that is dependent on the monochromatic 2500A UV luminosity, $L_{2500}$. We also observe some redshift dependence of the distribution. The mean of the $L_\text{X}$ distribution increases with redshift while the width decreases. This weak but significant redshift dependence leads to $L_{2500}$-$L_\text{X}$ and $L_{2500}$-$\alpha_\text{ox}$ relations that evolve with redshift, and we produce a redshift- and $L_{2500}$-dependent $\alpha_\text{ox}$ equation. Neither black hole mass nor Eddington ratio appear to be potential drivers of the redshift evolution.


INTRODUCTION
The energetic processes associated with the fuelling of Active Galactic Nuclei (AGN) produce radiation across the electromagnetic spectrum.An optically thick accretion disc is expected to emit thermally resulting in a blackbody across the optical/UV (Shakura & Sunyaev 1973).Meanwhile, the bulk of the X-ray emission is thought to be produced by inverse Compton scattering of accretion disc photons accelerated to X-ray energies in some form of corona following a power-law spectrum.The geometry of the corona is unclear.Various models exist to describe this corona: from a lamp-post geometry where the corona illuminates the disc from its position above the black hole (Fabian et al. 2017), to a slab corona that sandwiches the disc (Haardt & Maraschi 1991).X-ray polarimetry with the Imaging X-ray Polarimetry Explorer (IXPE; Weisskopf et al. 2022) has begun to constrain the geometry of the corona in a handful of AGN (MCG-05-23-16 [Marinucci et al. 2022;Tagliacozzo et al. 2023], IC 4329A [Ingram et al. 2023], NGC 4151 [Gianolli et al. 2023]).The presence of a 'soft excess' (i.e., X-ray emission ≲1 keV exceeding what would be expected from an extrapolated power-law) suggests an additional component to the X-ray production and is often attributed to an inner warm disc (Petrucci et al. 2018(Petrucci et al. , 2020)).
The relationship between the X-ray and UV luminosity has been known for some decades (Avni & Tananbaum 1982, 1986) and parametrized as  X ∝   UV with  ∼ 0.6.The relationship is generally considered to be tight (Lusso & Risaliti 2017;Bisogni et al. 2021), although some scatter is observed, motivating models where the processes involved in producing the X-ray and UV emission are ★ E-mail: amy.rankine@ed.ac.uk (ALR) dependent on some common parameter of the AGN (e.g., accretion rate, black hole mass; Lusso & Risaliti 2017;Kubota & Done 2018).There is little evidence of evolution with redshift (e.g., Vignali et al. 2003;Steffen et al. 2006;Just et al. 2007;Green et al. 2009;Lusso & Risaliti 2017;Timlin III et al. 2021); however, see Shen et al. (2006) and Kelly et al. (2007) who do see some redshift-dependence of the relation.In fact, the lack of significant redshift evolution and the general tightness of the relation has lead to claims that the correlation between the X-ray and UV luminosities (or more precisely the X-ray and UV fluxes) can be used to infer cosmological parameters (Salvestrini et al. 2019;Lusso et al. 2020).
The spectral index of a power-law between the UV luminosity and X-ray luminosity, specifically the monochromatic 2500 Å and 2 keV luminosities (first introduced by Tananbaum et al. 1979), is denoted  ox and is often used to parametrize the  2500 - X relationship.Jin et al. (2023) have recently shown that  2500 is appropriate as a tracer of the accretion disc emission and that it is sufficient as a single parameter to describe the optical/UV emission over a broader wavelength range such that  2500 is suitable for determining the relation of the UV to X-ray emission.The non-flat relationship between  ox and  2500 shows that the X-ray luminosity increases less than monotonically as UV luminosity increases suggesting that the spectral energy distribution (SED) becomes more disc-dominated.The physical driver of this relation is unclear and revealing the true relation between  2500 and  X free from selection effects will be a step towards understanding the physical mechanism(s) that govern the relations.
The observational  2500 - X and  2500 - ox relations are plagued by selection effects due to both the choice of the parent quasar sample and the limitations of the available X-ray data.Many studies

DATA
We use X-ray data from XMM and UV/optical data from SDSS both taken in the XXL and Stripe 82 fields.In brief, we are using the optically-selected quasars from SDSS DR16 (Lyke et al. 2020) at redshifts 0.5 <  < 3.5 across the two regions and have re-reduced the XMM data using the xmmpype custom pipeline outlined in Georgakakis & Nandra (2011).We crossmatch the XMM sources with the SDSS sources using Nway (Salvato et al. 2018).Detailed descriptions of each dataset are provided below; however, some readers may wish to peruse Table 1 and move onto Section 3.

Optical/UV
We use the SDSS DR16 quasar catalogue (Lyke et al. 2020) to create the optically-selected sample of AGN within the XXL and S82 fields.We filter the DR16 quasar catalogue with the multi-order coverage maps (MOCs) of the XXL and S82 XMM observations (see Section 2.2) in Aladin to select only the objects within SDSS that fall within the footprints of XXL and S82.The quasar catalogue is further limited to the optically-selected quasars which we define as the CORE sample from the BOSS and eBOSS targets (Myers et al. 2015).The CORE sample is produced by selecting objects with the following of SDSS's Bitmasks activated: bit 40 (QSO_CORE_MAIN) of mask BOSS_TARGET1, bit 10 (QSO_EBOSS_CORE) of EBOSS_TARGET0, and bit 40 (QSO1_EBOSS_CORE) of EBOSS_TARGET1.With this selection, we aim to only include quasars that were selected and targeted based on their optical properties.In doing so, we avoid biasing our results by including, for example, the X-ray selected quasars in the XXL field which were observed as part of the large SDSS ancillary programme led by A. Georgakakis.
The optical sample contains both X-ray detected and undetected objects (see Section 2.2) with a total of 2292 quasars.Note that our selection does not remove quasars with broad absorption lines in their spectra (BAL quasars;Weymann et al. 1991) or radio-loud (i.e., jetted) quasars (Kellermann et al. 1989).This choice was made in order to assess the X-ray properties of the truly optically-selected quasar population; removing them would introduce additional selection biases.However, we acknowledge that BAL quasars tend to be X-ray weak compared to non-BAL quasars (e.g., Gibson et al. 2009;Luo et al. 2014) and the X-ray emission of radio-loud quasars can be dominated by jets (Shang et al. 2011;Zhou & Gu 2021).BAL quasars can only be identified at  > 1.5 at which redshifts any potential C iv absorption systems are within the observed wavelength window of SDSS, meaning that we are only able to identify BAL quasars in 54 % of our sample.There are only 58 quasars identified as BAL quasars in our optical sample, of which 5 are X-ray detected.Radio-loud quasars number only 32 in the sample (Lyke et al. 2020, based on having a match to FIRST in the SDSS DR16 quasar catalogue;), with 17 of these detected in the X-ray.It is unknown if there are additional radio-loud quasars in our sample that would be detected with deeper data, further justifying our decision to not apply a radio cut on our sample.

Optical/UV properties
The SDSS spectra are reconstructed using the ICA technique outlined in Rankine et al. (2020) which essentially provides high S/N versions of the spectra over the restframe wavelength range 1260-3000 Å from which the continuum luminosity at restframe 2500 Å can be measured,  2500 .The left-hand panel of Fig. 1 contains an example spectrum and reconstruction. 2500 is estimated by calculating the median flux in a 10 Å window centred on 2500 Å and converting to a luminosity.Where available, redshifts from Rankine et al. (2020) which are based on an independent component analysis (ICA) of the optical spectra are used, otherwise, the redshifts reported in Lyke et al. (2020) are used.The differences between the two redshift samples are of order 300 km s −1 with only a few as different as ∼1000 km s −1 .The updated redshifts from Rankine et al. (2020) will not significantly affect the calculations of luminosities; however, they will produce more accurate black hole mass measurements, particularly C iv-based masses due to the correction derived from the 'blueshift' of the emission line (see Section 6).All in all, the changes are minimal.
We correct the luminosities for Galactic dust extinction with the dustmaps Python module (Green 2018) and the dust map of Schlegel et al. (1998) updated by Schlafly & Finkbeiner (2011) in tandem with the extinction module (Barbary 2016) and the reddening curve of Fitzpatrick (1999), producing median  ( − ) = 0.02 for the XXL sample and 0.03 for S82 and S82X.2500 Å is redshifted out of the BOSS spectrograph at  ≳ 3.2 ( ≳ 2.7 for the SDSS spectrograph) which would ordinarily prevent the measurement of the 2500 Å monochromatic luminosity of quasars above this redshift.However, reconstructing the spectra with the ICA technique which utilises the spectral information, including emission lines and the continuum shape, across the rest of the available spectrum above 1260 Å allows the 2500 Å luminosity to be estimated reliably.We checked the accuracy of extrapolating the reconstructions with a sample of quasar spectra in which 2500 Å was present but only included the wavelength range 1260-2200 Å in the fitting and found good agreement with the reconstructions that used the full available wavelength range between 1260-3000 Å. See an example of this extrapolation in the right-hand panel of Fig. 1.Only 46 quasars of our optically-selected sample require extrapolation of the reconstructions.
Uncertainties on  2500 are calculated by propagating the errors on the weights of the ICA spectral components produced during the reconstruction process.The median errors on  2500 for the subset of objects without restframe 2500 Å in their spectra and  2500 was extrapolated from the reconstructions are ∼0.04 dex compared to ∼0.02 dex for the subset with restframe 2500 Å which reflects the indirect measurement of  2500 .Errors from the spectrum reconstructions will be much less than those from the spectrophotometry; however, we do not propagate the  2500 errors further, since the main source of uncertainty is the X-ray luminosities, and so do not make an attempt to quantify them here.
The left panel of Fig. 2 shows the distribution of  2500 versus redshift for the X-ray detected and undetected quasars.The CORE Stripe 82 and Stripe 82X samples contain very few quasars above  ∼ 2.2 compared to the XXL sample due to the differing SDSS selection between SDSS II and SDSS III/IV with all of the CORE Stripe and Stripe 82X quasars originating from SDSS II.

X-ray
We start from the 294 XMM pointings in the North field of XXL (Pierre et al. 2016).XMM-XXL North covers ∼25 deg 2 with an exposure time of 10 ks per XMM pointing.
Stripe 82 is an equatorial region of sky covering ∼300 deg 2 which has been repeatedly observed with SDSS.Approximately 28 deg 2 of Stripe 82 has been observed with XMM.This combines the 198 pointings from the Stripe 82X survey (S82X) at ∼5 ks per XMM pointing and 33 additional archival pointings (S82; 7-66 ks per pointing) extracted from the XMM archive (LaMassa et al. 2013;LaMassa et al. 2016).

Reduction
We use the xmmpype XMM pipeline, which is based on the methods and techniques described in Georgakakis & Nandra (2011).In brief, the pipeline creates images in the different energy bands, sources are detected and astrometric corrections are applied before X-ray fluxes are estimated and any optical counterparts to the X-ray sources are identified.One advantage of employing the pipeline is the greater accuracy of the sensitivity curves which are generated with a robust and well-quantified Bayesian approach (following the methods of Georgakakis et al. 2008).The sensitivity curves allow for an accurate characterisation of the selection function of a sample using analytic relations instead of cumbersome and computationally expensive simulations and can naturally account for non-detected sources.
In particular, at faint fluxes the Bayesian sensitivity curves correctly account for the effects of Poisson statistics on the X-ray detection and photometry in the low-counts regime and the impact of Eddington bias. Figure 3 contains the area curves for the S82, S82X, and XXL fields in the full band.In general, at a given flux, the XXL sample is most sensitive, followed by the S82 archival pointings and finally the S82X survey.The nature of our investigations means that correcting for the X-ray detection probability is necessary and will be most significant at faint fluxes.Additionally, the pipeline coadds overlapping XMM observations to increase the X-ray depth.It is also designed for large-area serendipitous X-ray surveys which greatly facilitates the post-processing of the various products in the case of surveys that extend over large sky areas.We limit our sources to those detected in the full band (0.5-10 keV) where a detection is defined by a "false detection probability"  false < 4 × 10 −6 , where  false is the probability of the observed counts (or higher) being produced purely by a fluctuation of the background.Column 1 of Table 1 lists the number of X-ray point sources resulting from the reduction of the XMM pointings, totalling 14 493 sources.Comparing to the S82 reductions of LaMassa et al. (2016), we find 5529 X-ray sources in the combined S82 regions, whilst LaMassa et al. ( 2016) produced a catalogue of 4668 sources with XMM detections in the full band.We find that the log -log  relations of Georgakakis et al. (2008), the ExSeSS catalogue (Delaney et al. 2023), and the CDWFS (Masini et al. 2020) are in good agreement with those of our sample (see Fig. 4) providing confidence in the source detection and sensitivity maps of the xmmpype reductions (see Appendix A for comparisons in the hard and soft bands).

Crossmatching
We perform an initial search for possible optical counterparts in SDSS DR16 (Ahumada et al. 2020) with xmatch (Pineau et al. 2020) and a search radius of 40 arcsec around each X-ray source which yields an optical catalogue of 1 484 651 sources.We use Nway to match the X-ray observations to this catalogue with a 20 arcsec maximum radius.X-ray RA and Dec positional uncertainties were generated during the reduction with median uncertainties of ∼1.5 arcsec.We supply constant 0.1 arcsec positional uncertainties for the optical catalogue.We supply Nway with the total sky area of the reduced XMM observations -calculated from the multi-order coverage maps (MOCs) generated by xmmpype -and estimate the sky area of the input optical catalogue by creating a MOC with Aladin (Bonnarel et al. 2000) and a radius around each X-ray observation of 40 arcsec, producing a total area of 14.23 deg 2 once overlaps Median luminosity errors are presented in the legend.The 1-D redshift and luminosity distributions for the three samples are plotted above and to the right, respectively, of their corresponding axes.The  2500 panel contains both X-ray detected and undetected quasars, whereas the  X panel contains only the X-ray detected subsample.
between 40-arcsec regions have been accounted for.Nway produces a matched catalogue containing all possible matches for each X-ray source and corresponding probabilities. any is the probability that an X-ray source has a true counterpart in the provided catalogue and   is the probability that a particular match is the true counterpart.
As such, a combination of  any and   and limits on each can be invoked to produce a final catalogue of robust optical counterparts of the X-ray sources.Nway calculates the average source density on the sky from the provided sky areas which leads to a scaling of the counterpart probabilities,  any and   .Combining the XXL and Stripe 82 fields leads to an average sky density across the two fields which will affect the relative counterpart probabilities for sources in different fields.However, in our use case of Nway we do not use any absolute  any or   thresholds to determine the final matches; instead we are only ever comparing  any and   values between different objects across small physical scales; i.e., optical sources that are potential matches to the same X-ray source such that they are within the same region.As such, the scaling of the probabilities does not affect the final matching.
We include magnitude priors in the crossmatching to preferentially select counterparts with optical magnitudes that match the magnitude distribution of quasars which are less likely to be spurious alignments and more likely to be the true counterparts to the X-ray sources.We perform the matching with -band information from SDSS with priors pre-determined based on the magnitudes of the optical quasars in the optical catalogue compared to the non-quasar objects.
We make the final match selection by prioritising counterparts that are classed as AGN which we define as either having spCl (the spectroscopic class) as 'QSO' or if the object is found in the SDSS DR16 quasar catalogue (Lyke et al. 2020) having performed a simple 1-arcsec crossmatch between the DR16 and DR16Q catalogues.To implement the AGN-prioritisation, we take the possible matches from Nway, and inspect the match with the highest product of  any and   that is also an AGN.The product avoids multiple X-ray sources having the same AGN optical counterpart and gives priority to the X-ray source with the highest probability of having a counterpart in this optical catalogue ( any ).Only using   would result in 195 X-ray sources having an optical match already associated with another X-ray source.Only if the   for this optical AGN is > 0.01  of the original best match is the AGN selected as the counterpart.We perform a false-positive calibration by offsetting the X-ray positions and running Nway with this mock X-ray catalogue (and corresponding mock optical catalogue obtained with xmatch and a 40 arcsec search radius).The AGN number density on the sky is low such that for our AGN-prioritisation scheme a  any threshold of zero is sufficient to maintain a false-positive fraction <1 %.
Ultimately, we obtain the highest completeness when including the -band quasar-based magnitude prior; however, the QSOprioritisation scheme leads to only a few X-ray matches changing depending on the prior used.We end up with 26 % of all our X-ray sources having an optical counterpart that is spectroscopically identified as an AGN.Given that we are starting from an optically-selected subsample of the SDSS DR16 quasar catalogue, we limit the sample to the X-ray sources that have optical counterparts identified as AGN based on their inclusion in the DR16 quasar catalogue.Our final sample thus contains 2292 optically-selected AGN, 770 (34 %) of which are X-ray detected (see Table 1).

X-ray properties
X-ray flux measurements for the full 0.5-10 keV band are calculated during the reduction with Galactic absorption taken into account (estimated from the H i maps of the LAB survey; Kalberla et al. 2005) but assume a photon index of Γ = 1.4.We are specifically selecting X-ray sources associated with (broad-line) quasars and so expect them to have unabsorbed X-ray spectra.We check this using the hardness ratios, defined as where  and  are source counts in the hard (2-10 keV) and soft bands (0.5-2 keV) normalised by exposure, and confirm that they are, on average, consistent with Γ = 1.9 (see Fig. 5).We convert the flux measurements to a Γ = 1.9 using conversion factors from webpimms based on H i column densities of 2 × 10 20 and 3 × 10 20 cm −2 for The median errors are plotted in the bottom right.The horizontal lines mark the average hardness ratios for Γ = 1.4,1.6, and 1.9 assuming appropriate H i column densities of 2 × 10 20 and 3 × 10 20 cm −2 and PN and MOS detectors.
the XXL and S82 fields, respectively.We apply a K-correction when calculating rest-frame luminosities that also assumes a photon index of Γ = 1.9.The X-ray luminosity distribution with redshift is plotted in the right-hand panel of Fig. 2. The sensitivity of the X-ray data is apparent in the lower bound on  X with redshift.The S82 and XXL samples are similar in their  X - distributions; however, XXL extends to higher redshifts due to the SDSS selection.

MEASUREMENTS OF THE INTRINSIC DISTRIBUTION OF 𝑳 X AS A FUNCTION OF 𝑳 2500 AND REDSHIFT
We aim to arrive at a model that describes the distribution of X-ray luminosity as a function of UV luminosity and redshift.In Fig. 6 we plot the distribution of X-ray luminosity for our X-ray detected quasar sample in bins of  2500 and  (solid colour histograms).In what follows we will make use of the Bayesian sensitivity curves provided by xmmpype in order to account for the X-ray undetected quasar population and derive the underlying  X distribution function.

Completeness-corrected distribution at a given 𝑳 X
We attempt to account for the undetected X-ray sources in each ( X ,  2500 , ) bin by calculating the probability of a source having an X-ray luminosity  X given its  2500 and : The numerator is the number of X-ray detected quasars in each ( X ,  2500 , ) bin.The denominator takes into account the probability that quasar  with redshift   would be detected if it had an X-ray luminosity corresponding to the centre of the  X bin and is summed over all X-ray detected and undetected quasars in that ( 2500 ,  bin).
Δlog 10  X is the width of the  X bin.The corrected counts in a given ( X ,  2500 , ) bin can then be calculated by the following: . (3) In the limit where (det| X , ) = 1 for all quasars in a given  X bin (i.e. the X-ray data are sufficiently deep that any quasar with that  X should be detected),  corr =  det and thus corresponds to the "uncorrected" (solid) histograms in Figure 6.Thus, as expected, at the highest  X the corrected (open histograms/error bars) and uncorrected (solid histograms) estimates are consistent.
The corrected counts are plotted in Fig. 6 as coloured outlined histograms and Poisson errors are generated based on the  det in each bin and applying Gehrels' method for small number statistics (Gehrels 1986).As expected, the correction is larger at low X-ray luminosities.The corrected  X distribution is perhaps Gaussian with the centre, , and width,  potentially varying with  2500 and ; however, this model can only correct bins with  det > 0 and significant binning is required.Additionally, while narrower bins leads to higher resolution, the uncertainties on  X are comparable to the  X bin-width.In the following section we move on to using Maximum Likelihood Estimation (MLE) to arrive at a fully unbinned approach to determining the  X distribution as a function of  2500 and .

Maximum likelihood fitting
The observed and corrected distributions in Fig. 6 suggest that log 10  X is normally distributed for quasars of a given  2500 and : with mean  and width  both of which may depend on  2500 and/or .In this section, we will attempt to fit the X-ray luminosity distribution function from equation 4 via maximum likelihood estimation (MLE) and will investigate the requirement for  2500 -and -dependence.The log-likelihood (which we derive in Appendix B) is given by ln where the first and second terms account for the X-ray detected and undetected quasar samples, respectively.The black histograms are a random sample drawn from the assumed Gaussian distribution of  X with parameters determined by the maximum likelihood estimation with the best-fitting model (vii) (see Section 3.2).From both the binned corrected counts and the MLE results, it is clear that the X-ray detected sample is skewed towards the high  X sources.In the majority of the  2500 and  bins the stacked  X from the MLE results agrees with the stacked data (black and coloured vertical arrows with 1- error bars).Only the ( 2500 , ,  X ) bins populated with X-ray detected quasars can be corrected via the binning method outlined in Section 3.1, providing motivation for the MLE detailed in Section 3.2).
number of counts from a source with  X which is determined via The energy conversion factor, ECF, exposure,  exp , background counts,   , and total counts,   are specific to each X-ray detection and calculated during the reduction.The encircled energy fraction, EEF, is 70 % for our adopted aperture and is based on the point spread function at 2 keV, the average energy (weighted by the response) of the full band.The luminosity distance,   (  ), and Kcorrection,  corr (  ), are calculated for the redshift   of the quasar.Since the X-ray luminosities (fluxes more precisely) and errors are calculated from a Poisson distribution by xmmpype, (  | exp ) will be maximal when the integration variable  X equals the estimated X-ray luminosity of quasar ,  X  .We note that the maximum of (  | exp ) corresponds to our nominal best estimate of the X-ray luminosity,  X  , for a given detected quasar.The X-ray undetected term, similarly to the detected term, depends on the probability of undetected quasar  having an X-ray luminosity  X given its UV luminosity  2500  and redshift   , ( X | 2500  ,   , ).This probability is multiplied by the probability of quasar  remaining undetected if it were to have X-ray luminosity  X : where (det| X ,   ) is calculated from the area curves (Fig. 3) via The integration limits are set as log 10  X = 40, ∞.In practise the upper limit is set by the maximum X-ray flux probed by the sensitivity curves (10 −10 erg s −1 cm −2 ).

Distribution of 𝑳 X in fixed 𝑳 2500 and redshift bins
We first aim to determine if and how the X-ray luminosity distribution changes as a function of  2500 and .We divide X-ray detected and undetected quasar samples between equally spaced redshift bins: 0.5 <  < 1.5, 1.5 <  < 2.5, and 2.5 <  < 3.5.We also split the samples across six  2500 bins such that there are approximately equal numbers of quasars in each bin.For each (  ,  2500  ) bin we fit for  and  by maximising the log-likelihood in equation 5 with the Python package emcee (Foreman-Mackey et al. 2013).The best-fitting parameters are presented in Fig. 7 as circles. is clearly dependent on  2500 with the mean  X increasing with increasing  2500 across all redshift bins.On the other hand, there is little evidence for a  2500dependent  at any redshift but it is possible that  decreases as redshift increases suggesting that the distribution of  X is narrower at greater redshifts.In light of these correlations, we remove the bins contain data and so some bins are missing from the analysis.The lines and shaded regions correspond to the  and  parameters obtained when modelling a linear dependence of  2500 on  and no dependence for  (model (ii)).Significant trends with  2500 exist for  in all redshift bins and the results from models (i) and (ii) match well.Within the errors, models (i) and (ii) agree with a  2500 -independent  but a clear decreasing  with increasing redshift.
2500 binning in the next section and model the relationship between  2500 and  (and ) as linear.

𝑳 2500 -dependent distribution of 𝑳 X in fixed redshift bins
When binning by  2500 the model parameter  (i.e. the average of the log  X distribution) appears to increase as  2500 increases.We model this dependence on  2500 for  and  as linear with log 10  2500 : We perform MLE on the -binned data to constrain   ,   ,   , and   (model (ii)), thus removing the need to bin our quasar sample according to  2500 .It is not clear that  varies with  2500 , thus we repeat the MLE for the model where  does not depend on  2500 , formally   = 0 therefore  =   (model (iii)).
To compare the different models (with different numbers of free parameters) we will use the Akaike Information Criterion defined as with  dim the number of free parameters and L the maximum of the likelihood function (Equation 5).The AIC penalises models with a large number of parameters and models with lower AICs are considered to better represent the data.To calculate the AICs of the models with binning, we treat the model as a piecewise function such that the maximum log-likelihood is the sum of the maximum log-likelihood over all  bins (and  2500 bins for model (i)) and  dim is the total number of parameters across all bins.Model (iii) is formally a better fit with a lower AIC than model (ii) (see Table 2) and so we plot the results of model (iii) in Fig. 7 as the straight lines and shaded regions.Within the errors the linear dependence on  2500 agrees with the  2500 -binned fits of model (i) (circles; Section 3.3).Across redshift bins, the intercept of the  relation and  in general change.At a given  2500 ,  increases as redshift increases and  decreases.This can be seen more clearly in Fig. 8.The gradient of the - 2500 relation appears to be relatively constant with redshift; however,   is clearly increasing as redshift increases and  is decreasing.We thus move on to model the dependence of these parameters on redshift to arrive at a fully unbinned MLE in Section 3.5.

Continuous model of the redshift evolution
In this section we arrive at a selection of models to describe the whole data sample in a continuous manner instead of discrete  2500 or  bins.We model any possible redshift evolution of  and  via a linear dependence of   ,   , and  on .We perform the MLE with various models with different  dependencies, explicitly: (iv) no redshift evolution, with parameters   ,   , and .This is the equivalent of model (iii) in the limit of one redshift bin.
(v) no redshift evolution, with parameters   ,   ,   , and   .This is the equivalent of model (ii) but assuming a single, broad redshift bin.
(vi) only redshift evolution of , with gradient and intercept parameters for   () and   (), and a constant .
(vii) redshift evolution of only   and  with gradient and intercept parameters for   () and () and a constant   .
From the AIC values in Table 2, the model which best represents the data is model (vii) which allows for redshift evolution of   and  parametrized by, We also make use of nested sampling via the MultiNest algorithm (Feroz & Hobson 2008;Feroz et al. 2009Feroz et al. , 2019) ) and its Python implementation PyMultiNest (Buchner et al. 2014) for model comparison assessed via the Bayesian evidence which also prefers model (vii) (see the 7th column of Table 2).In fact, there is greater evidence for the models with redshift evolution (vi-viii) compared to those without (iv, v).The grey lines and shaded regions in Fig. 8 are the best-fitting parameters for model (vii) (also listed in Table 3).The -binned   values from Section 3.4 are systematically higher than the continuous redshift modelling in this section.This is due to the distribution of objects within the relatively broad redshift bins and the intrinsic redshift evolution of  within such bins.The higher  2500 and  X sources within a given redshift bin are preferentially Increasing the number of redshift bins by a factor of two removes this systematic bias but also reduces the number of objects in each bin and thus leads to greater statistical uncertainties in the parameters.

UNDERLYING 𝑳 2500 -𝑳 X DISTRIBUTION
With model (vii) in hand, for a given , as expected the peak of the intrinsic  X distribution increases as  2500 increases.Perhaps not as obvious is that for a given  2500 , the intrinsic  X distribution shifts to higher  X as redshift increases (  increases since the gradient of   () is found to be positive) and also narrows ( decreases since the gradient of () is negative).We compare the underlying distribution of  X of our opticallyselected quasar sample to the observed data and original binned corrections in Fig. 6.For each object in our sample with a given  2500 and redshift, detected or otherwise, we draw 100 samples from the distribution function (equation 4) with the best-fit parameters listed in Table 3, effectively creating a mock sample of  X measurements if there were no limitations in X-ray depth.In Fig. 6 we then normalise to the number of objects in each  2500 and redshift bin (black histograms).Unlike the binned corrections, we can infer the source counts in  X bins with zero observed sources.In the majority of  2500 and  bins the binned corrections and the MLE corrections agree.However, in some panels (e.g., second row, second to last column) the binned corrections are significantly lower than the MLE distribution which we believe to be the combination of using the Bayesian sensitivity curves which appropriately account for Eddington bias but are not suitable for the crude binned corrections carried out in Section 3.1 in bins where the majority of the sample is around where redshift evolution is modelled as a linear dependence of  on   ,   , and .The grey squares represent the higher resolution redshift binning used to check the cause of the systematically higher binned points in   .
the flux limit.As noted previously, in all bins, the detected sources are only probing the high  X tail of the distribution.
As mentioned in Section 1, there is a well-known correlation between the optical and X-ray luminosities of AGN with  X ∝   UV and  ∼ 0.6.In this work, we have found that the X-ray luminosity distribution is a function of  2500 and redshift, with the peak of the distribution given by ( 2500 , ).In the left panel of Fig. 9 we plot the peak of the  X distribution for a constant  = 1, 2, 3 where an increase in redshift produces a higher  X for a given  2500 .
In order to check our results, we produce a stacked value of  X from the X-ray counts extracted at the positions of all of our quasars.We do so by calculating individual X-ray luminosities for each optical source and then produce a mean  X .This produces an  X value (black square in the left panel of Fig. 9) that is higher than the centre of the contours due to the mode of a log-normal distribution (which the  X distribution is) being different from its mean.As a sanity check, calculating the average  X in the same way from the mock data used to produce the contours results in a higher  X value than the contours would suggest (red cross); however, it is consistent with the stacked  X from the data.We do the same in the middle panel of Fig. 9 in bins of  2500 and  to compare to the relations.Again, we find that the  X values from the stacked data are systematically higher than the relations and, although suffering from small-number statistics with this relatively high-resolution binning, the stacked mock data is in agreement.In fact, the stacked data in the highest redshift bin (red squares) appear to agree too well with the relations; however, this is due to the distribution of redshifts within this redshift bin.If instead we were to plot the relations for the mean redshift within each redshift bin, the red line would shift lower and the squares would be offset.The redshift dependence on the width of the  X distribution is also having an effect: at low redshifts where the distribution is wider, the discrepancy between the stacked data and the relation is greater.This in turn reduces the redshift dependence in the stacked data.
In order to compare more directly to the literature, we take the mock sample (red contours in the left-hand panel of Fig. 9) and fit a straight line and we obtain a  ≃ 0.62 which is in agreement with the literature (right-hand panel of Fig. 9).It is not obvious that this best-fit line is in agreement with the contours; however, the median log 10  X in  2500 bins is consistent with the  ≃ 0.62 relation (blue points and error bars in the right panel of Fig. 9).

CORRECTED 𝜶 OX
The spectral slope between the X-ray and optical,  ox , is often used as a means of describing the relationship between the X-ray and UV luminosities, and is calculated as follows, where  2 keV is the monochromatic X-ray luminosity at 2 keV and corresponding frequency  2keV . 2500Å is the frequency equivalent to 2500 Å.We calculate  ox for our X-ray detected sample by converting the full band  0.5−10 keV into  2 keV and plot these values against  2500 in Fig. 10 colour-coded by redshift.The flux-limited nature of the parent sample is clear here in that the highest  2500 quasars are only found at high redshifts; however, the completeness curves generated from the sensitivity curves in Fig. 3 reveal that the completeness of our optically-selected sample drops significantly as X-ray luminosity decreases ( ox decreases) across the full range of  2500 probed by our quasar sample.
In what follows, we make use of the derived underlying X-ray luminosity distribution as a function of both  2500 and redshift to produce a corrected  2500 - ox relation.We calculate  ox for the mock sample in Section 4 with equation 13 and produce the red contours in Fig. 10.The true underlying  ox distribution suggests that we are missing the  2500 -moderate,  X -faint population which any  ox relation should account for.
In order to check our results, we produce a stacked value of  ox .We take the stacked  X value from Section 4 for all of the quasars, log this value and convert to  ox with the mean log 10  2500 of our data.This produces an  ox value (black square in Fig. 10) that is higher than the centre of the contours due to the mode of a lognormal distribution (which  ox is) being different from its mean.Calculating the average  ox in the same way from the mock sample used to produce the contours results in a higher  ox value than the contours would suggest; however, it is consistent with the stacked  ox from the data.We caution that simple linear stacked measurements to infer relations with broad, log-normal shapes will not correspond to the peak (mode) of the distribution but be biased high as we have found here.
The relationship between the peak of the  ox distribution and the  2500 and redshift is given by where  = −0.264+0.013 −0.013 ,  = 0.159 +0.013 −0.013 and  = 6.095 +0.400 −0.395 .In short, the relation is derived by converting the peak of the full-band  X (0.5-10 keV), given by model (vii) with the parameter values from Table 3, to the 2 keV monochromatic luminosity and substituting this in equation 13.The full derivation is presented in Appendix D. Thus far we have not considered whether the parameters of our model are independent; however, parameters  and  are correlated since both are functions of   and so the uncertainty on  ox will include, at the very least, the covariance of  and .We consign the equation for Δ ox and its derivation to Appendix D but note here that we assume the parameters of our model in Table 3 are independent (but see Appendix D and Fig. D1).We provide the posterior distributions of the parameters as supplementary data.
In Fig. 11 we plot  ox from equation 14 for a constant  = 1, 2, 3 where an increase in redshift produces a vertical shift towards lessnegative  ox values (upwards on the plot).Our model that describes how the peak of the intrinsic distribution of  ox depends on  2500 at different redshifts (accounting for X-ray sensitivity limits and the underlying redshift evolution of the relation between  2500 and  X over this redshift range) produces significantly steeper relations (solid lines in Fig. 11) than most prior estimates that use X-ray upper-limits and often combine a wide redshift range (e.g., Just et al. 2007;Nanni et al. 2017;Timlin III et al. 2021, as shown by the dashed lines in Fig 11,right).For comparison to the literature, we use our mock sample that corrects for the X-ray incompleteness (but not the uneven sampling of the quasar samples in terms of  2500 and ) and fit a linear relation between  2500 and all of our mock  ox values.Fitting the mock sample with a single linear relation produces an  ox with a flatter slope that is in better agreement with the literature but has a lower normalisation.The black line is lower in normalisation for two reasons: i) it accounts for X-ray fainter sources that tend to be below the sensitivity limits, and ii) it tracks the quasar sample that is dominated by lower redshift ( ≲ 2) sources, which we find to have lower  ox (at a given  2500 ).
We consider the effect of a soft excess which is thought to be important at restframe energies below 1 keV (Halpern 1984;Arnaud et al. 1985;Done et al. 2012) Thus it can only be observed in our lowest redshift sources (0.75 keV at  = 0.5).None the less, we conducted simulations with xspec (Arnaud 1996), finding that at  = 0.5 we could be overestimating the 2 keV luminosities by only 10 % which is insignificant when compared to the redshift dependence suggested by our results.Additionally, we see no evidence for a soft excess in the hardness ratios of our sample (Fig. 5).
Although we are aiming to produce an  2500 - ox relation corrected only for observational selection effects and are not considering any intrinsic absorption, we investigate if the observed redshiftdependence can be explained by intrinsic absorption of the X-ray emission.We draw intrinsic log  H values uniformly between 20 and 22, and Γ from a normal distribution centred on 1.9 with a standard deviation of 0.2.We produce X-ray spectra with these pa-rameters using xspec and calculate the observed 2 keV luminosity from the observed 0.5-10 keV luminosity.At  = 0.5, we could be underestimating the 2 keV luminosity by 20 %, 8 % at  = 2, and 5 % at  = 3.5.Although there is a systematic underestimation that is correlated with redshift, the discrepancies are again insignificant compared to the redshift-dependence we observe.Additionally, rerunning the maximum likelihood analysis with the sample of Peca et al. (2023), who performed spectral analysis of the Stripe 82 sample to account for intrinsic absorption in their calculations of X-ray luminosities, produces still a redshift-dependent relation. 1Hard-band X-ray emission will be less affected by absorption than the full band and so we re-run our whole analysis (including Nway matching) on the X-ray sources detected in the hard band and find that model (vii) is still the most successful model in explaining the data.Models (vi) and (viii) are also given a viable joint-second place, but importantly, both of these models involve redshift-dependence.The source numbers are smaller in the hard band, thus we choose not to use this for our main analysis.See Appendix C for the summary statistics using the hard band.The SDSS spectra of our quasars also show little evidence of any intrinsic extinction; thus any changes in intrinsic reddening across our redshift range are negligible and cannot explain the observed redshift evolution of the  2500 - X relation.

DISCUSSION
Using a sophisticated Bayesian framework, we have shown that the intrinsic distribution of X-ray luminosities of the SDSS quasar sample evolves with redshift, shifting toward higher  X at a given  2500 and with decreasing scatter in the distribution as redshift increases.Our finding is in disagreement with a number of prior works that do not find any evolution in this relation, albeit for distinct samples 1 While their source extraction method differs from ours resulting in a different and smaller sample this comparison provides at least a first order test of the effect of intrinsic absorption.
Figure 11.The same  2500 - ox space as Fig. 10 with the X-ray detected sample as black circles and mock sample described by the contours.Left: The solid blue and red lines are calculated via equation 14 with constant  = 1, 2, 3 and the shaded areas are the 0.5- widths of the relation.Right: The dashed lines are relations from the literature.The orange line is produced by fitting a straight line to the mock sample which is consistent with the median  ox in  2500 bins (blue points and error bars).and without applying the sophisticated analysis techniques that we present here (see Section 1).However, Kelly et al. (2007) also find evolution of the  ox - relation in a sample of radio-quiet quasars across  = 0.1-4.7 with  ox increasing as redshift increases (with  ox depending linearly on the age of the Universe).Additionally, Shen et al. (2006) also perform a similar maximum likelihood analysis using soft X-ray detection from RASS in the SDSS DR3 and found a redshift dependent  2500 - X relation, albeit weaker than found here. 2n a quick glance, our redshift-dependent  2500 - X relation would suggest that this relation cannot be used as a cosmological probe as Lusso & Risaliti (2017) suggest, for example.However, we want to stress that our results are applicable only to the overall optically-selected quasar population.Only with carefully chosen sub-samples of quasars can it be possible to use the  2500 - X for cosmological purposes (Salvestrini et al. 2019;Bisogni et al. 2021).Regardless of whether quasars can be used in this way to test cosmological models (see Khadka et al. 2023, who suggest the answer is uncertain), our aim is to eventually use this relation to determine the underlying causes of the link between X-ray and UV emission in the larger quasar population.
As mentioned above, the purpose of this work is to derive the intrinsic distribution of  X as a function of  2500 and redshift that applies to the optically-selected SDSS quasar sample, specifically.While our finding of redshift evolution of the  2500 - X relation is at odds with the consensus (see Section 1) we do not dwell on complex comparisons as our results apply to a specific (but well-defined) sample.However, one advantage of our work is that we have carefully considered the X-ray sensitivity limitations, without so doing would result in a different answer for the  2500 - X relation (or  ox ).These relations are important for understanding the balance of energetic output coming from the corona versus the accretion disk, and in order to compare to physical models (e.g., Kubota & Done 2018) one must account for X-ray sensitivity limitations of the sample.
While we do not aim to come up with a detailed physical model to explain the observed redshift evolution, it is informative to look at the black hole properties of the optically-selected SDSS quasar sample across redshift and compare to the trends observed in Kubota & Done (2018).We estimate the black hole masses (BHM) and Eddington ratios ( Edd ) of our quasars using the ICA-based spectrum reconstructions (see Fig. 1 and Section 2.1), calculating BHMs from the full width at half maximum (FWHM) of the C iv1550 and Mg ii2800 emission lines, redshift-permitting, and the 1350 Å and 3000 Å monochromatic luminosities.For the C iv-derived BHMs we apply the relation of Coatman et al. (2017) which accounts for the non-virial component and subsequent asymmetry of the C iv emission line.The Mg ii BHMs are estimated with the Vestergaard & Osmer (2009) relation.We apply bolometric corrections of 5.15 and 3.81 to the monochromatic 1350 Å and 3000 Å luminosities, respectively, to estimate the bolometric luminosities (Shen et al. 2011).With the BHMs and bolometric luminosities in hand, the  Edd are calculated.We calculate the median BHM and  Edd in bins of redshift along with the standard error on the median values and the standard deviation of the distributions using either the C iv-and Mg ii-derived values or both values at redshifts where both lines are within the spectral window (see Fig. 12).We do not focus on the absolute values of the quantities but instead focus on the general trends of BHM and  Edd with redshift, observing that neither BHM Intrinsic X-ray luminosity distribution 13 or  Edd show any significant trend with redshift across the majority of the  2500 bins.
The spectral energy distribution (SED) model of Kubota & Done (2018) predicts that the X-ray luminosity scales linearly with BHM, after fixing  X = 0.02 Edd motivated by the SED fits of a handful of AGN (but applied to a larger sample by Mitchell et al. 2023).In fact, Mitchell et al. (2023) explicitly show that  ox (and therefore the  2500 - X relation) has a dependence on both the BHM and Eddington-scaled accretion rate due to the relative contributions of the hot X-ray corona, the warm Compton region of the disc and the standard disc in their truncated disc model.Additionally, Kubota & Done (2018) observe that an increase in the Eddington-scaled accretion rate of 1 dex should produce a decrease in log X of 0.5 dex (for constant  2500 ).The large uncertainties on our BHM and  Edd measurements preclude a more detailed discussion on the effect of the BH properties and it is unclear if BH mass and/or Eddington rate are responsible for the observed redshift dependence of the  2500 - X relation in our sample.In the future, it would be valuable to extend the Bayesian hierarchical modelling to consider the underlying Eddington-scaled accretion rate and black hole masses, enabling a more direct comparison with Kubota & Done (2018).
Quasar SEDs likely evolve with BHM and  Edd .This true evolution coupled with the flux-limited nature of SDSS leads to quasar samples that have masses and  Edd distributions that appear to evolve with redshift.This is implied by the shift to larger BH masses and Eddington ratios as  2500 increases from left to right in the panels of Fig. 12.However, these measurements help provide insight into the physical origins of the observed trends.Modelling of the optical selection effects will be the focus of a future work.Nevertheless, the lack of significant correlations between the BH properties -masses and  Edd -and redshift does not appear to be consistent with BHM and/or  Edd being responsible for our observed evolution; however, we caveat this with the fact that our measurement uncertainties are large and increase with redshift.Another selection effect in the optical quasar sample or another physical parameter could be responsible for the observed evolution of the relations with redshift.
One property of the BHs that we have not considered is their spin.Higher spins have been found to lead to greater X-ray emission relative to the UV emission (Temple et al. 2023).For spin to produce our observed increase in  X as  increases would require the BH spins of our optically-selected sample to increase with redshift.Cosmological simulations suggest that BH spins can evolve through mergers and/or accretion resulting in higher spins at high redshifts (e.g., Dubois et al. 2014); however, Temple et al. (2023) suggest that the spins of SDSS quasars at  ∼ 2 are generally low.Regardless of the cause of the redshift-dependence, a single  2500 - ox relation with a constant slope across all redshifts is almost certainly not correct and so physical models should not be trying to reproduce a non-evolving  ox relation.

CONCLUSIONS
We have carefully inferred the intrinsic X-ray luminosity distribution as a function of UV luminosity and redshift of the optically-selected SDSS quasars in the Stripe 82 and XXL fields using a sophisticated Bayesian hierarchical modelling approach.We have crossmatched the optical SDSS sample to the XMM point sources with Nway (Salvato et al. 2018).We have combined XMM-detected quasars with Bayesian sensitivity curves calculated with the custom xmmpype pipeline (Georgakakis & Nandra 2011) in order to extract information from the X-ray undetected quasars.Our main findings are: (i) The xmmpype reductions produce log -log  curves that are consistent with previous works (Fig. 4 and Section 2.2) (ii) The  2500 - X relation can be modelled as a Gaussian function with mean  which depends on the log 2500 and width  (Section 3.2).
(iii) There is some redshift dependence of the  2500 - X relation with  increasing with redshift., on the other hand, decreases as  increases (Section 3.5 and Fig. 8).
(iv) For a constant , our fitted  2500 - X relation has a slope of  ∼ 0.3.The slope in the observed relation of  ∼ 0.6 found by previous works is reproduced when considering the joint redshift and  2500 distribution of the optically-selected SDSS quasar sample (Section 4 and Fig. 9).
(v) Measurements from stacked X-ray data should be considered with caution when deriving quantities from a log-normal distribution (Sections 4 and 5).
(vi) Attempting to correct the X-ray luminosity distribution in  X bins to account for the undetected quasars can lead to underestimated source counts and is limited to only X-ray luminosity ranges that have been detected (Fig. 6 and Section 3.1).A more sophisticated estimation of the  X distribution is implemented via the Bayesian hierarchical modelling approach used throughout the rest of the paper.
(vii) We produce a relation to describe  ox that is now a function of  2500 and redshift.When marginalising over redshift in our SDSS sample, the  ox relation we recover has a slope consistent with the literature but with a lower normalisation (Section 5 and Fig. 11).
We have made the first steps to understand the intrinsic relationship between the X-ray and UV luminosity by considering the opticallyselected SDSS quasar sample.The next step is to approach the problem from an X-ray selected sample in order to parametrize the optical selection.The X-ray selected sample from eROSITA (Merloni et al. 2012;Predehl et al. 2021) with follow-up spectroscopy from SDSS-V (Kollmeier et al. 2017) will be beneficial for this work and support a broader goal of obtaining a full characterisation of the UV and X-ray emission properties of the AGN population and the underlying physical structure of the accreting system that produce them.(B7) The first term of equation B7 is the probability of observing   photons given that for source  we have measured data   and propose that it has an X-ray luminosity  X , which does not depend on  2500  or  and thus reduces to (  | X ,   ,   ).This term captures the uncertainty in the value of  X of the source based on the fact that an integer number of counts,   , were detected; we marginalise over the range of possible  X .With a change of notation, (  | X ,   ) is equation 6 with  exp as the expected number of photons from a source with  X and is calculated via equation 7. The second term of equation B7 is the prior expectation for  X given the observed  2500  and   , and model parameters describing the distribution of  X , independent of the X-ray data (and thus we can drop   ).The resulting term is the model we aim to fit (equation 4).Thus, we find Equations B8 and B14 have the same second term (other than different  and  subscripts).The total likelihood for all of our data is given by the product of L  of all detected objects, and L  of all undetected object: (B15) The log-likelihood can therefore be written as, which is equation 5.

APPENDIX C: HARD BAND ANALYSIS
We run our analysis using the same optically-selected sample but with the hard-band X-ray detections and adopting the corresponding sensitivity curve for non-detections.Of the 2292 parent sample, 255 (11 %) are detected in the hard band.Due to the much smaller X-ray detected sample, we only run the maximum likelihood estimation for the unbinned models.Table C1 contains the AICs and Multinest evidence for this run.Note that the AIC and relative log 10  values cannot be compared across different bands as they use different data and have been normalised using their respective best-fit models.The most successful model is model (vii) which is the same as for the full band.Models (vi) and (viii) also prove to fit the data well with the former being preferred by the Multinest Bayesian evidence and the latter by the AIC.Importantly, all three of these models allow redshift evolution.In Fig. C1 the model parameters for model (vii) using the hard band data are consistent with those for the full band data in the main paper.The consistency between our hard-band results, presented here, and the full-band analysis used in the main paper demonstrates both the robustness of our Bayesian analysis when applied to smaller, shallower samples (as is the case for the hardband sample) and that X-ray absorption effects are not driving these results and our observed redshift evolution (as any impact would be severely reduced when using a harder band).

APPENDIX D: DERIVATION OF 𝜶 OX RELATION
By combining model (vii) and equation 13 we derive the relationship between the peak of the distribution of  ox and the  2500 and redshift.The peak of the  X distribution, in units of erg s −1 , is given by  =   log 10  2500 − 30 +    +   .The monochromatic 2 keV

Figure 1 .Figure 2 .
Figure 1.Example reconstructions of quasar spectra using the ICA technique for restframe 1260-3000 Å (top) and the residuals (observed spectrum flux − reconstruction) normalised by the noise (bottom).The left panels contain a  = 2.08 quasar with full coverage of the 2500 Å region.The right panels demonstrate the extrapolation of the reconstruction for a  = 3.47 quasar without coverage of the 2500 Å region.The red shaded area is the 1- uncertainties on the reconstruction.

Figure 3 .
Figure 3. Sensitivity curves for the full band (0.5-10 keV) across the three regions in our sample.

Figure 4 .
Figure 4. Top: cumulative number counts as a function of full band XMM flux and comparison to ExSeSS (Delaney et al. 2023), CDWFS (Masini et al. 2020), the reduction of Stripe 82 by LaMassa et al. (2016), and the model ofGeorgakakis et al. (2008).Bottom: differential number counts with the Euclidean slope removed.Errors are Poisson errors based on the number of sources and scaled accordingly.All measurements have been converted to the 0.5-10 keV band assuming a consistent Γ = 1.4 spectrum.

Figure 5 .
Figure5.Hardness ratios versus redshift for our X-ray detected quasar sample.The median errors are plotted in the bottom right.The horizontal lines mark the average hardness ratios for Γ = 1.4,1.6, and 1.9 assuming appropriate H i column densities of 2 × 10 20 and 3 × 10 20 cm −2 and PN and MOS detectors.

Figure 6 .
Figure 6.Distribution of  X binned by  2500 , .Each panel is a different  2500 (columns, increasing to right) and  (rows, increasing towards bottom) bin.The X-ray detected quasars are presented as the filled histograms.The binned corrected counts and associated Poisson errors are represented by the open histograms and error bars (see Section 3.1).The black histograms are a random sample drawn from the assumed Gaussian distribution of  X with parameters determined by the maximum likelihood estimation with the best-fitting model (vii) (see Section 3.2).From both the binned corrected counts and the MLE results, it is clear that the X-ray detected sample is skewed towards the high  X sources.In the majority of the  2500 and  bins the stacked  X from the MLE results agrees with the stacked data (black and coloured vertical arrows with 1- error bars).Only the ( 2500 , ,  X ) bins populated with X-ray detected quasars can be corrected via the binning method outlined in Section 3.1, providing motivation for the MLE detailed in Section 3.2).

Figure 7 .
Figure 7. Model parameters  (top) and  (bottom) from equation 4 and estimated via MLE as a function of  2500 .The parameters for model (i), which requires running the MLE on data binned by  2500 and , are shown by the circles and 1- error bars.Colours represent  bins.Not all (,  2500 )bins contain data and so some bins are missing from the analysis.The lines and shaded regions correspond to the  and  parameters obtained when modelling a linear dependence of  2500 on  and no dependence for  (model (ii)).Significant trends with  2500 exist for  in all redshift bins and the results from models (i) and (ii) match well.Within the errors, models (i) and (ii) agree with a  2500 -independent  but a clear decreasing  with increasing redshift.

Figure 8 .
Figure8.Breakdown of the parameters  (top and middle) and  (bottom) as a function of redshift for the -dependent models.The top two panels are the gradient   and the value of  at log 10 ( 2500 ) = 30.The squares are the parameter values used to produce the straight-lines in Fig.7from model (ii).The grey lines and shaded regions are the parameter values obtained with model (vii) where redshift evolution is modelled as a linear dependence of  on   ,   , and .The grey squares represent the higher resolution redshift binning used to check the cause of the systematically higher binned points in   .

Figure 9 .Figure 10 .
Figure9.The 2500 - X plane with the X-ray detected sample plotted as black circles.Left: The solid blue and red lines are the peaks of the  X distribution function with constant  = 1, 2, 3 and the shaded areas are the 0.5- width of the distribution.The mock sample is represented by the red contours.The black square is the result from our XMM stacking analysis and red cross is the equivalent stacking of the mock sample.Middle: The same relations and data as in the left panel, now overplotted with the stacked data (squares) and stacked mock sample (crosses) in redshift and  2500 bins.Right: The orange line is produced by fitting a straight line to the mock sample, and the blue points and error bars are the median  X in  2500 bins and the 1- of the  X distributions, respectively.

Figure 12 .
Figure12.Black-hole mass (top) and Eddington ratio (bottom) versus redshift for our optically-selected quasar sample in bins of increasing  2500 (panels left to right).Black hole masses are derived from either the C iv1550 or Mg ii2800 emission lines (blue and orange points, respectively) depending on whether the line falls within the spectral coverage at a given redshift.The numbers in the legend refer to the total number of quasars in each luminosity bin with coverage of the C iv and Mg ii lines.The error bars mark the standard error on the median (solid) and the standard deviation of the sample in each bin (dotted).We observe no obvious trend in either BHM or  Edd with redshift for a given  2500 .

µ
Figure D1.Posterior distributions of the model parameters.Some parameters show correlations, namely   with   , and   with   .

Table 1 .
Number counts for final samples of X-ray detected and undetected sources for the S82, S82X and XXL fields.The first column contains the total number of point sources with detections in the full band extracted with xmmpype.The second and third columns contains the number of opticallyselected quasars that have X-ray counterparts (are X-ray detected) and those that do not (undetected).

Table 2 .
The models fitted with MLE.Columns are, in order, model number for referring to in text; the binning required of the model and whether or not there is redshift evolution for the completely unbinned models; the parameters of the model; the number of dimensions of the model which takes into account the number of redshift and  2500 bins; the AIC values; ΔAIC is the difference between the AIC for that model and the lowest AIC value.log 10  is the PyMultiNest Bayesian evidence for each of the unbinned models normalised by the most likely model (highest evidence).Model (vii) with the lowest AIC is presented as bold.
identified toward higher redshifts and thus in our binned results a steeper relation between  X and  2500 (i.e. a steeper   ) is recovered to account for this redshift evolution.The intercept,   , does not have such a strong dependence on the width of the redshift bin.
), Max-Planck-Institut für Extraterrestrische Physik (MPE), National Astronomical Observatories of China, New Mexico State University, New York University, University of Notre Dame, Observatário Nacional / MCTI, The Ohio State University, Pennsylvania State University, Shanghai Astronomical Observatory, United Kingdom Participation Group, Universidad Nacional Autónoma de México, University of Arizona, University of Colorado Boulder, University of Oxford, University of Portsmouth, University of Utah, University of Virginia, University of Washington, University of Wisconsin, Vanderbilt University, and Yale University.Intrinsic X-ray luminosity distribution 17 the detection was first made (usually with different extraction radii), then this extra step necessitates that the detection probability should remain.In our case, the xmmpype reduction uses the same data to determine if a detection meets the detection criteria and then calculates luminosities.The second term can be expanded as follows: (  ,   ,  2500  ,   |) = (  |  ,  2500  ,   , ) × (  ,  2500  ,   |), (B3) but as   ,  2500  , and   do not depend on  we can drop the second term of equation B3.Next, we introduce the marginalisation over  X via,  X .We have, L  = (  |  ,  2500  ,   , ), ∫ (  ,  X |  ,  2500  ,   , ) d X (B6) = ∫ (  | X ,   ,  2500  ,   , )( X |  ,  2500  ,   , ) d X .
∫(  | exp )( X | 2500  ,   , ) d X .(B8)Movingonto the likelihood for an undetected source , L  is the probability that object  is undetected with data D  for the model parameters :L  = (det, D  |) (B9) = (det,   ,  2500  ,   |),  X ,   ,  2500  ,   |) d X (B11) = ∫ (det,  X |  , 2500  ,   , )(  ,  2500  ,   |) d X   ,  2500  ,   , )( X |  ,  2500  ,   , ) d X .The second term in equation B12 can be dropped for the same reasons as the second term in equation B3.Equation B13 results from expanding the first term of equation B12.Equation B13 is introduced in order to consider the probability of object  remaining undetected for a proposed  X .The first term of equation B13 can be simplified as (det| X ,   ,   ) and is calculated from the sensitivity curves (see equations 8 and 9) for simplicity we remove   since   is absorbed in the sensitivity curves.The second term reduces to ( X | 2500  ,   , ) (dropping   ) and it is again the prior expectation of  X for a given  2500  and   .Explicitly, L  ∝ ∫ (det| X ,   )( X | 2500  ,   , ) d X .