{'fraction_non_alphanumeric': 0.03149114364351688, 'fraction_numerical': 0.019772723435151882, 'mean_word_length': 4.582956259426847, 'pattern_counts': {'":': 1, '<': 20, '<?xml version=': 0, '>': 7, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Disclosing the Radio Loudness Distribution Dichotomy in Quasars: An Unbiased Monte Carlo Approach Applied to the SDSS-FIRST Quasar Sample

We investigate the dichotomy in the radio loudness distribution of quasars by modelling their radio emission and various selection effects using a Monte Carlo approach. The existence of two physically distinct quasar populations, the radio-loud and radio-quiet quasars, is controversial and over the last decade a bimodal distribution of radio loudness of quasars has been both affirmed and disputed. We model the quasar radio luminosity distribution with simple unimodal and bimodal distribution functions. The resulting simulated samples are compared to a fiducial sample of 8,300 quasars drawn from the SDSS DR7 Quasar Catalog and combined with radio observations from the FIRST survey. Our results indicate that the SDSS-FIRST sample is best described by a radio loudness distribution which consists of two components, with 12+/-1 % of sources in the radio-loud component. On the other hand, the evidence for a local minimum in the loudness distribution (bimodality) is not strong and we find that previous claims for its existence were probably affected by the incompleteness of the FIRST survey close to its faint limit. We also investigate the redshift and luminosity dependence of the radio loudness distribution and find tentative evidence that at high redshift radio-loud quasars were rarer, on average"louder", and exhibited a smaller range in radio loudness. In agreement with other recent work, we conclude that the SDSS-FIRST sample strongly suggests that the radio loudness distribution of quasars is not a universal function, and that more complex models than presented here are needed to fully explain available observations.


INTRODUCTION
An important property of Type 1 AGNs (i.e. broad line quasi stellar objects; QSOs) is the existence of radio-loud (RL) and radio-quiet (RQ) populations. One of the more controversial topics in QSO studies is whether these RL and RQ quasars 9 form two physically distinct populations of objects. On the one hand, both types of quasars are likely powered by similar physical mechanisms (e.g. Barthel 1989;Urry & Padovani 1995;Shankar et al. 2010), and their radio loudness has been shown to anticorrelate with accretion rates (in Eddington luminosity units) onto their central supermassive black holes (e.g. Sikora et al. 2007). On the other hand, it has been demonstrated that, relative to RQ quasars, RL quasars are likely to reside in more massive host galaxies (e.g. Peacock, Miller, & Longair 1986;Sikora et al. 2007), and harbor more massive central black holes (e.g. Laor 2000;Lacy et al. 2001;McLure & Jarvis 2004). However, mainly due to two problems, i) severe selection biases, such as incompleteness in observed QSOs samples, and ii) the overwhelmingly high fraction of RQ quasars, it is still unclear whether RL and RQ quasars form two distinct populations of objects, or a continuous sequence (as suggested by e.g. Lacy et al. 2001). Here we focus on this problem.
Whether RL and RQ quasars form two distinct populations can be studied by investigating the relation between their radio and optical emissions. This illuminates the relative importance of the likely dominant sources of electromagnetic radiation in these two wavelength windows, namely the relativistic jet and the accretion disk. Strittmatter et al. (1980) first pointed out that the radio-to-optical flux density ratio for optically selected QSOs appears bimodal (the so-called 'quasar radio dichotomy'), which suggests that QSOs can be divided into two distinct, RL and RQ, populations of objects. Other authors found additional evidence for a dichotomy (e.g. Kellermann et al. 1989;Miller et al. 1990;Ivezić et al. 2002;White et al. 2007;Zamfir et al. 2008). However, several studies (e.g. White et al. 2000;Cirasuolo et al. 2003b;Lacy et al. 2001) disputed its existence. A bimodal distribution in quasar radio loudness would point to distinct physical properties of radio-loud and -quiet QSOs, such as a different physical origin of radio emission (jet vs. corona; Laor & Behar 2008), different black hole masses (e.g. McLure & Jarvis 2004), accretion rates (e.g. Sikora et al. 2007;Hamilton 2010) and/or spins (e.g. Blandford & Znajek 1977;Blandford 1990;Garofalo et al. 2010), as well as host galaxy properties (e.g. Sikora et al. 2007;Lagos et al. 2009;Kimball et al. 2011). The existence of two distinct quasar populations may then be linked to hierarchical structure growth in a ΛCDM universe in which dark matter halos and galaxy mergers play an important role (e.g. Hughes & Blandford 2003;Lagos et al. 2009;Garofalo et al. 2010;Hamilton 2010;Fanidakis et al. 2010).
The radio-to-optical flux density ratio as a measure of radio loudness was initially proposed by Schmidt (1970). A division between RL and RQ quasars at a ratio of ∼10 was found by Kellermann et al. (1986) using restframe 5 GHz (6 cm) radio and B-band optical flux densities. An alternative radio loudness definition 10 is purely based on radio luminosity, e.g. L 1.4 GHz 10 25 W Hz −1 (Miller et al. 1990). Studies of large QSO samples have yielded that the most RL QSOs are ∼ 10 3 times more luminous in radio than in optical (e.g. Ivezić et al. 2002). On the other hand, much deeper radio studies of smaller samples showed that RQ QSOs typically have a factor of ∼ 10 weaker radio than optical emission (e.g. Kellermann et al. 1989). Using SDSS (York et al. 2000) and FIRST (Becker, White & Helfand 1995) data, Ivezić et al. (2004b) have shown the existence of a peak at the high end of the radio-to-optical flux ratio distribution. However, because of the detection limits of both surveys, they could not detect a significant number of RQ objects at radio wavelengths, and hence argued the existence of a secondary peak at the low end (and thus bimodality) based on two arguments: i) the majority of SDSS quasars (∼ 90%) were not detected in the FIRST survey, thus they should lie at the low end of the radio loudness distribution, and ii) deep radio studies showed that typical RQ QSOs have a factor of ∼ 10 4 weaker radio emission than typical RL QSOs, thus a secondary peak must exist.
One of the most recent results on the quasar radio dichotomy comes from a stacking analysis of FIRST 20 cm snapshot images at the positions of ∼ 40, 000 quasars from the SDSS DR3 catalog (White et al. 2007). They showed that a shallow minimum exists between the RL and RQ parts of the radio loudness distribution; however, it was stressed that optical selection effects probably dominate the observed distribution of radio loudness. Specifically, the authors discuss the difficulty of identification of SDSS quasars at 2.4 < z < 3, where their colors are very similar to stellar colors (see e.g., Richards et al. 2002), and support their claim by noting the strengthening of the bimodality for a subsample of quasars in that redshift range. Cirasuolo et al. (2003b) have used Monte Carlo simulations of the quasar population in order to model the intrinsic radio loudness distribution. They have compared their simulated samples to three quasar samples (2dF QRS, LBQS and PBQS; see Cirasuolo et al. 2003b for details) with FIRST survey data and radio observations by Kellermann et al. (1989). The three samples probe a wide range of radio loudness, but contain fewer than 200 radio-detected quasars in total. The authors reported that the best-fit radio-to-optical flux ratio distribution is a double-Gaussian function with ∼97% of the quasars in the RQ, and ∼3% in the RL component. However, they conclude that there is no minimum between the two Gaussians, thus no bimodality exists (but see Ivezić et al. 2004b for a different interpretation).
The main problems causing the discrepancies in literature regarding the existence or absence of a quasar radio loudness dichotomy lie either in the ambiguity of quasar selection, low number statistics, or severe selection biases linked to flux-limited samples (because radio loudness is defined as a radio-to-optical flux ratio, i.e. a ratio of two quantities drawn from flux limited samples). To properly address all of these biases, in the work presented here we utilize a Monte Carlo based approach similar to that in Cirasuolo et al. 2003b, but more robust and using a much larger sample of observed quasars (the largest QSO database currently available is the SDSS quasar catalog).
The outline of the paper is as follows. In Section 2 we present the data and our sample. In Section 3 we outline our methodology and the simulation algorithms. We present our results in Section 4, and discuss them in Section 5. We summarize our conclusions in Section 6. Throughout the paper we report magnitudes in the AB system, for the optical as well as for the radio. We use standard cosmology (H 0 = 70 km s −1 Mpc −1 , Ω M = 0.3, Ω Λ = 0.7) and quasar continuum spectrum defined as f ν ∝ ν α , where ν denotes frequency, f ν is flux density and α is the spectral index.

Choice of the Source Catalog
For a robust study of the radio loudness of quasars it is necessary to use a large, statistically significant radiooptical sample of quasars that covers the broadest possible range in radio loudness. As quasars are rare objects (e.g. the space density of SDSS quasars with i < 19 is ∼11 deg. −2 ) a statistical sample can be assembled via large-area observations. On the other hand, achieving radio sensitivities over these fields deep enough to probe far into the radio-quiet regime with the current generation of radio interferometers is challenging, and unfeasible over fields larger than a few degrees.
Because of the scaling of both radio and optical depth limits with surveyed areas, state-of-the-art surveys, such as e.g. SDSS-FIRST (∼9380 deg. 2 , i < 22.5, S 1.4GHz > 1 mJy ;Schneider et al. 2010;Becker, White & Helfand 1995) Bondi et al. 2003) cover a comparable radio loudness range, while smaller area surveys suffer from small number statistics in their optical quasar samples (see e.g. Smolčić et al. 2008). Hence, for the analysis presented here we have utilized the largest available sample of quasars with radio coverage, drawn from the SDSS and FIRST sky surveys.
For the reasons outlined above, we use the SDSS Seventh Data Release Quasar Catalog (Schneider et al. 2010). It contains 105,783 spectroscopically confirmed quasars taken from ∼9,380 deg. 2 of the sky. The catalog consists of objects with spectroscopy that 1) have reliable redshifts, 2) have at least one emission line with FWHM greater than 1000 km s −1 or interesting/complex absorption features, and 3) are more luminous than restframe M i = −22.0 (in our cosmological model, assuming a power-law continuum with α = −0.5). The selection is based on sources' i magnitude and position in the multi-dimensional color-space based on five SDSS photometric bands (see Richards et al. 2002 for details). The quasars in the sample have 15.0 < i < 21.2, but the majority of quasars are brighter than i ≈ 19 because the flux limit for the main spectroscopic sample is i < 19.1. In addition to the main spectroscopic sample, the catalog contains serendipitously identified quasars and sources selected for their proximity (within 2") to a 1.4 GHz radio source drawn from the FIRST survey. The rest of the quasars in the sample are matched to the FIRST catalog to within a 2" radius. The radio flux densities at 20 cm (1.4 GHz) are given in "AB radio magnitudes": where f rad is the radio flux density at 20 cm in Jy.
FIRST is a radio survey at 1.4 GHz/20 cm, conducted with the Karl G. Jansky Very Large Array' (VLA) in B-configuration (Becker, White & Helfand 1995). It has mapped approximately 10,000 deg. 2 of the North Galactic Cap with a beam size of 5.4" and a typical RMS sensitivity of 0.15 mJy/beam. Flux density values used here are integrated over the two-dimensional Gaussians fitted to each source. Due to the lack of very short spacings in the VLA B-configuration and the nature of the Gaussianfitting source detection algorithm, fluxes for extended objects larger than about 10" are likely underestimated and split into multiple components in the FIRST source catalog. We estimate that this does not significantly affect our sample as most of the radio-detected quasars are expected to be unresolved at the angular resolution of the FIRST survey. Multiple component sources are rare and in most cases radio-loud, which puts them high above the interesting transition region between RQ and RL regimes (e.g. Jiang et al. 2007). For a discussion of the distribution of integrated and peak flux densities and the extended source bias of FIRST-detected SDSS quasars we refer the reader to Kimball & Ivezić (2008). Nonetheless, we do take this effect into account statistically, using the FIRST survey completeness correction derived specifically for SDSS quasars (see Figure 1 in Jiang et al. 2007, and references therein), taking into account both the source size and flux distribution.

Our Main Sample
Following the SDSS DR7 Quasar Sample notation, hereafter we adopt the apparent radio magnitude (t) for radio flux densities. The catalog lists 8,630 quasars with a radio detection within 2" from the optical source position. For 'radio luminosity' here we use the 1.4 GHz/20 cm absolute radio magnitude (M t ), derived from the apparent radio magnitude using a K-correction of the form −2.5(1 + α) log(1 + z) (e.g. Richards et al. 2006), assuming α = −0.5 for each quasar (e.g. Kimball & Ivezić 2008). For the 'optical luminosity' we use the emission-line-corrected absolute magnitude in the restframe SDSS i band (λ ef f = 7471Å). The reason for this choice is that the SDSS i-band was used to construct the flux-limited DR7 quasar sample and that i magnitudes are available even for quasars at the highest redshift. Prior to computing the absolute i magnitudes for the quasar continua, we have corrected the cataloged apparent magnitudes for galactic extinction using the extinction maps of Schlegel, Finkbeiner & Davis (1998). To get continuum magnitudes, we subtract the contribution of emission lines to the i band using a function derived by Richards et al. (2006). Using the mean SDSS quasar spectrum, they have calculated the contributions of major quasar emission lines (above the power-law continuum) as a function of redshift, up to z ≈ 5. After the subtraction of this contribution to the i-band apparent magnitudes, we use the canonical K-correction with α = −0.5. Given the high uncertainty of this correction for z ≥ 5, hereafter we exclude 56 quasars with z ≥ 5 from our sample. For details of the calculation see Section 5 and Figure 15 in Richards et al. (2006). In order to get a 'cleaner' optical selection, we exclude the quasars selected only on the basis of their proximity to a FIRST source. After the exclusion of those quasars, the z > 5 quasars and those with uncertain photometry, we are left with 8,307 quasars in our main radio-optical subsample. Our main optical sample consists of 98,544 quasars.

SIMULATIONS
If one had a complete volume-limited quasar sample with both rest-frame optical and radio luminosities, studies of the radio loudness distribution would be simple. For example, if radio and optical luminosities were linearly correlated (regardless of the physics of such a setup), the radio-to-optical luminosity ratio distribution would be very narrow and the radio versus optical luminosity diagram would show a straight narrow band. Another possible situation could be a broad distribution of the radio-to-optical luminosity ratio, with perhaps evidence for its two-component nature (dichotomy), with or without evidence for a local minimum (bimodality). In this case, the radio versus optical luminosity diagram would show a broadly dispersed line or two possibly overlapping lines, one for the RQ and one for the RL quasars (see the bottom two panels of Figure 1 for an illustration). If the sampled redshift range was sufficiently broad, one could even perform such studies for subsamples selected from narrow redshift slices, and search for evidence of evolution in the inferred properties with cosmic time.
Unfortunately, such a complete sample does not exist. The main difficulty when working with real samples is that both optical and radio data are truncated at finite flux levels, and this strong selection effect must be properly accounted for. There are additional effects, such as K-corrections and optical variability, which can be handled to some extent, as discussed further below. In order to utilize existing samples, various statistical methods have been used in published work to account for the observational effects, as briefly reviewed in Section 1.
Here we use Monte Carlo simulations, with the following main logical steps: 1. The main sample is defined at optical wavelengths, with various selection effects, and the main observable is the optical apparent magnitude. The sample also contains redshifts for all quasars and radio magnitudes for a subset of radio-detected quasars. The selection function for the sample is well known.
2. Only a subset of real optical sources from the main sample are detected at radio wavelengths, compris-ing the main radio-optical subsample. They have a certain distribution of radio magnitudes and, at least in principle, a different distribution of optical magnitudes than the main optical sample. The distributions of radio and optical magnitudes for the radio-optical subsample are the main constraints on models for the relationship between the radio and optical luminosities.
3. Given a quasar from the main optical sample with certain optical luminosity, a parametrized model generates its radio luminosity (as detailed below). Using its implied apparent radio magnitude and the radio selection function, this object is either retained or rejected from the simulated radio-optical sample.
4. Starting with the main optical sample, the simulation produces a subset of quasars with simulated radio magnitudes -a simulated radio-optical subsample. The distribution of simulated radio magnitudes for this subsample, as well as its corresponding distribution of optical magnitudes, is compared to observed radio and optical magnitude distributions of the main radio-optical subsample, and utilized in a χ 2 minimisation procedure to get the best-fit model parameters.
Although the Monte Carlo sample utilizes observed optical magnitudes, its radio magnitudes are generated stochastically. Therefore, there is no object-to-object correspondence between the real radio-optical subsample and the simulated radio-optical subsample. For a good model, we expect a statistical agreement for two main quantities: i) optical magnitude distribution, and ii) radio magnitude distribution between the observed and simulated radio-optical subsamples. For simulating the radio luminosities we initially use four models independent of redshift and optical luminosity and the complete main sample presented in Section 2.2. However, once we select the best model, we investigate the redshift evolution and dependence on optical luminosity for its parameters. The following subsection defines radio loudness and describes the models considered in this work. The rest of this section describes the simulation algorithm (Section 3.2), the optimization strategy (Section 3.3) and the evaluation of our method on purely artificial data (Section 3.4).

The Radio Luminosity Models
We use four different models for the relationship between the radio continuum luminosity at 1.4 GHz (L rad ) and optical continuum luminosity in the SDSS i band (L opt ). Their relationship is parametrized using radio loudness (R K ) defined as: (2) Here K is either 0 or 1 and it selects one of the two different families of models that correspond to two common definitions of radio loudness found in the literature. For K = 0 the distribution of radio loudness simply equals the distribution of radio luminosities (thus RL quasars are defined as being more luminous than some threshold radio luminosity; e.g. Peacock et al. 1986;Miller et al. 1990, also examined in Ivezić et al. 2002 andJiang et al. 2007), while for K = 1 radio loudness is defined as the logarithm of the radio-to-optical luminosity ratio (with RL quasars having this ratio greater than some threshold; e.g. Kellermann et al. 1989;Ivezić et al. 2002Ivezić et al. , 2004bCirasuolo et al. 2003a,b;Jiang et al. 2007). Each definition is suitable for a specific model of the relationship between radio and optical luminosities. We have considered four simple models for the radio luminosity: two in which it is independent of the optical luminosity and two in which it is directly proportional to the optical luminosity. In the former case radio loudness is properly defined with K = 0, R 0 = log (L rad ) and radio luminosity is modeled as L rad = 10 R0 . In the latter case, the proper definition of radio loudness is K = 1, R 1 = log (L rad /L opt ), hence the radio luminosity model is L rad = L opt × 10 R1 . In both cases, we examine single-and double-Gaussian distributions of radio loudness. The designations of the models in this paper and their basic descriptions (where PDF stands for probability density function) are: Examples of all four models are plotted in Figure 1. Models M1 and M3 are two-parameter models, with the free parameters being x (the mean of the Gaussian) and σ (the width of the Gaussian). Models M2 and M4 have five free parameters: two for each Gaussian (x 1 , σ 1 , x 2 and σ 2 ) and an additional parameter f , which determines the ratio between the integrals of the two Gaussians. In all cases the overall normalization is determined automatically by the requirement that the number of radiodetected quasars in the simulated radio-optical samples matches the number in the observed SDSS-FIRST sample.
Initially, we use these four models to simulate radiooptical samples based on the complete main optical sample of 98,544 SDSS quasars and match them to the radiooptical subsample. Later we search for redshift evolution and dependence of the radio loudness distribution shape on optical luminosity not by adjusting the models, but rather by separating the main sample into smaller samples constrained in redshift and optical luminosity. Possible trends can then be inferred by observing how best-fit model parameters change between different subsamples.

The Main Optical Sample
For all radio loudness models that we have considered, the simulated radio magnitude (t) was calculated for each optically detected quasar in the main optical sample presented in Section 2.2. This sample contains 98,544 optically detected quasars with all of the selection biases introduced by the DR7 Quasar Catalog, except for the radio-only selection which was eliminated from the catalog for this work. Since we compare our simulated radiooptical subsamples to the radio-optical subsample of that -Examples of the radio models (M1-4) considered in this paper plotted in the M i -Mt plane. The radio loudness distribution is a single Gaussian for models M1 and M3 and a double Gaussian for models M2 and M4. Note that, in general, the Gaussians forming the two-component radio loudness distribution may exhibit complete, partial or no overlap, regardless of the model. In models M1 and M2 the radio absolute magnitude (Mt) is independent of the optical absolute magnitude (M i ), while in models M3 and M4 Mt is proportional to M i for every quasar. For clarity axis labels are shown only for the bottom-left panel, but they are the same for all panels. The dotted lines show where Mt = M i . catalog (i.e. to the sources also detected in the radio), this biased sample is indeed the most valid optical sample for our simulations.
Another option for the optical sample would be a simulated sample produced from an empirical luminosity function. However, such an approach suffers from two main problems: i) a simulated optical sample would be free of selection biases (assuming selection bias was properly corrected in the construction of the luminosity function) and therefore we would additionally need to simulate the SDSS optical selection; and ii) observed optical counts are difficult to simulate with proper uncertainties included (e.g. from photometry uncertainty, conversion between photometric bands, K-corrections, optical variability etc.) unless the same uncertainties have been considered during the generation of the luminosity function (see e.g. La Franca & Cristiani 1997 for an approach that does involve substantial consideration of uncertainties).

Assignment of Apparent Radio Magnitudes and Application of the Radio Selection Function
In each simulation, we assign an apparent radio magnitude (t) to each quasar from the main optical sample. The apparent radio magnitude t is calculated from the absolute radio magnitude M t , which is modelled according to Equation (2) as where K takes on values of 0 or 1 and R K is a random variable drawn from either a Gaussian or a double-Gaussian probability distribution, depending on the model (see Section 3.1).
To convert from M t to t we use a K-correction of the form −2.5(1 + α) log(1 + z) (e.g. Richards et al. 2006), assuming α = −0.5 for each quasar (e.g. Kimball & Ivezić 2008), and including corresponding uncertainties. Thus, where DM (z) is the distance modulus and K c,rad (z) is the K-correction for α = −0.5, both dependent on quasar's redshift. The last term is a normally distributed random variable accounting for uncertainties in K-corrections and photometry. Its mean is zero and the standard deviation (σ t ) can be estimated by considering the scatter in the measured radio spectral indices (e.g. Kimball & Ivezić 2008). As the standard deviation in spectral indices is approximately 0.3, it follows from the definition of the K-correction that σ t ≈ 0.35 for z = 2, where the number density of quasars is maximal. This uncertainty completely dominates the small photometric uncertainty (typically 0.03-0.04 mag). For our K = 0 models (M1 and M2) Equation (3)  simplifies The absolute optical magnitude is where DM (z) is the distance modulus, K c,opt (z) is the K-correction and the remaining two terms account for the scatter in K-correction, photometry (normal distribution, n(σ i )), and variability (exponential distribution, e(σ v )). The optical spectral index distribution with an average of −0.5 and a standard deviation of 0.3 (Richards et al. 2006) yields σ i ≈ 0.35, taking into account that the optical photometric errors (typically 0.03 magnitude) are negligible. Optical variability scatter has an empirically determined exponential distribution with zero mean and σ v ≈ 0.2 (Ivezić et al. 2004a). Finally, inserting the expression for M i into Equation (6), the distance moduli cancel out, as well as the mean Kcorrections (because of our assumption of the same mean spectral index in the optical and radio bands), and the normally distributed uncertainties σ t and σ i add up in quadrature, yielding a convenient expression for t in our K = 1 models: Our simulation algorithm assigns radio magnitudes to optically detected quasars by drawing random numbers from appropriate distributions specified by Equations (5), for M1 and M2, and (8), for M3 and M4. Each simulated radio-optical sample (consisting of 98,544 quasars with apparent optical and radio magnitudes) is then subjected to the radio selection function in order to produce a smaller subsample of 'radio-detected' quasars corresponding to the observed radio-optical sample from SDSS and FIRST. All quasars fainter than the FIRST flux limit at 1 mJy are rejected. Some of the quasars above the FIRST flux limit are randomly rejected in order to mimic survey incompleteness. This is performed by randomly choosing and excluding a number of quasars in narrow t-magnitude bins, with the fraction of excluded objects being determined by the completeness function of the survey up to its flux limit (see Figure 1 in Jiang et al. 2007). At the end of this procedure, each simulated radio-optical subsample is one particular realization of the real radio-optical subsample drawn from SDSS and FIRST. The simulated radio-optical subsamples after radio selection have to consist of approximately 8,300 quasars to at least roughly match the observations, which sets the overall normalization for our models.

Evaluation of Goodness of Fit for the Simulated Samples
Each of the simulated radio-selected subsamples ("realizations" hereafter) is binned into 10 bins in t, i, z and R ′ = 0.4(i − t) distributions. We chose to examine these four distributions to provide better constraints for our models; in principle, weaker constraints, always consistent with those derived by using all four distributions, can be can be obtained with any subset of these distributions that includes either t or R ′ . Note that the R ′ distribution is a proxi to the radio loudness distribution in the case of K = 1; otherwise it is just an additional constraint on the relation between t and i magnitudes 11 . The bin values are added to a pool of realizations simulated with identical model parameters. The radio magnitude assignment, the selection procedure and the binning need to be repeated a large number of times for a given model and a set of its parameters in order to properly account for the stochasticity of particular realizations. The aim of this Monte Carlo procedure is to calculate the mean distributions of t, i, z and R ′ for a given set of model parameters and derive their expected variance. We calculate the mean and the standard deviation of the number of simulated objects in each bin (SN j and σ j in the equation below) from a large number of realizations and compute the total χ 2 with respect to the binned t, i, z and R ′ distributions of the real SDSS-FIRST data. The χ 2 is defined as where the sum runs over all 40 bins, while RN j is the number of real quasars in a particular bin and SN j and σ j are the mean and the standard deviation for that bin, determined from the simulations as described above.
Although the four binned distributions are not entirely independent of each other (in particular, R ′ is a linear combination of i and t), the tests on artificial data (see Section 3.4) have shown that the number of degrees of freedom may be assumed as if the bins were independent. Therefore, the number of degrees of freedom equals N − p − 1, where N is the total number of bins used for evaluation and p is the number of free parameters of a particular model. Hereafter we refer to one evaluation 11 Also note that this distribution is severely biased and cannot be directly used to infer the distribution of radio loudness; see Ivezić et al. (2002 for an explanation. of χ 2 from its distribution given a large number of realizations, as a "simulation" and abbreviate the number of realizations per simulation with N r . It is a general property of Monte Carlo simulations to converge to an average result only after the experiment is repeated a large number of times. We test the convergence by examining how the means χ 2 , standard deviations and minimum/maximum values for different simulations with same N r . This is illustrated in Figure 2 for model M4 (analogous results follow from any of the models considered here) with a fixed arbitrarily chosen set of its parameters and a range of different N r . Clearly, a larger number of realizations per simulation reduces the probability of getting outlying χ 2 values that significantly deviate from the mean χ 2 for the simulated model, but at the expense of computing time. We have found it more time-effective to work with a relatively small number of realizations in our optimization procedure and correct for the scatter in the χ 2 values just before deriving the marginal probability density functions for model parameters (see Section 3.3). Due to computing time constraints, we chose to calculate 1000 realizations per simulation as a reasonable compromise between accuracy and computing time. The calibration shown in Figure  2 shows that for N r = 1000 one should expect variations with standard deviation of ∼2%. For N r = 3000 the variation drops to 1.2%, which does not represent a significant improvement compared to N r = 1000. The computing time scales linearly with the number of realizations per simulation, so N r = 3000 would require three times more computing time than N r = 1000.

The Optimization Strategy
Since there is no a priori indication that the parameter space for our models is simple, e.g. that there is a single minimum of the total χ 2 , we began to search the parameter space with a random walk across a wide range of parameters. The algorithm we used to initially sample the parameter spaces with 2 to 5 dimensions, is a Metropolis algorithm (Metropolis et al. 1953). Several runs of the algorithm were used to identify the area around the global minimum of χ 2 for each of the models. The program would then proceed to evaluate the region on a regular grid, so that marginal probability distributions could be estimated. We have used three levels of refinement of the grid in order to be able to sample well the narrow region around the global minimum of χ 2 .
Since variations displayed in Figure 2 are expected to occur in the evaluation process, we needed to take into account the fact that some simulations resulted in χ 2 values as low as the minimum χ 2 value, although their mean χ 2 would be much larger. For all parameter space points for which χ 2 was found to be less than 5 σ χ 2 (recall that this is ∼2% for N r = 1000) above the lowest value found, the optimization would proceed to map the surrounding parameter space in more detail (typically, with a factor of 3 more resolution in all dimensions).
For computation of the marginal probability density functions for model parameters, we assigned them the probability normally associated with χ 2 values: where x is a vector (a set of coordinates) in the parameter space of some model. The marginal probability of a certain parameter value is then the sum of the probabilities over all other parameter space dimensions. The probability density functions were normalized a posteriori so that the sum of probabilities of each bin in the parameter values equals unity. To cope with the problem of the internal scatter in χ 2 values (see Section 3.2.3), we used a special procedure to compute realistic probability density functions which take into account that the computed χ 2 may vary according to Figure 2. First, we estimate the standard deviation of χ 2 from the calculations illustrated in Figure 2 and described in the caption; e.g. the fractional standard deviation for N r = 1000 is ∼2%. Then, we calculate the corrected marginal probability density functions by adding a normal random variate with standard deviation equal to 2% of the lowest χ 2 value to the previously computed χ 2 values, calculating the probability distributions in each case and repeating this procedure many times. This Monte Carlo method requires ∼1000 repetitions for the resulting mean marginal probability density function to converge. We derive the asymmetric error bars on parameter values as intervals containing ±34% of the total probability around the median of the marginal probability density function for each parameter.

Verification of the Approach on Artificial Datasets
In order to evaluate our method and better understand the possible uncertainties and biases it involves, we first simulated artificial datasets resembling real data in all aspects, differing only in the fact that their radio loudness distributions are exactly known. Treating them the exact same way as the observational data, we have re- Fig. 3.-The distribution of parameters derived from our parameter optimization procedure (X derived ) relative to the input parameters of our artificial datasets (X real ), normalized to the derived standard deviation of the reconstructed parameters (σ derived ). The data is displayed for 20 optimizations on artificial samples constructed with model M4, which has 5 free parameters. The vertical error bars represent Poissonian noise. The solid black line marks the Gaussian fitted to the data, which is a close match to the ideal Gaussian with σ = 1 (dashed line), verifying that with our method it is possible to reconstruct the model parameters as expected from statistics.

constructed their model parameter values and compared them to the input ones.
For a set of artificial datasets constructed with the radio models presented in Section 3.1, the parameters were reconstructed to within 3σ uncertainties in all trials. The distribution of the reconstructed parameters, normalized to their respective standard deviations, is displayed in Figure 3. The lowest χ 2 values reached in optimizations were χ 2 ≈ 35, yielding a reduced χ 2 of ∼1 if the problem has the number of degrees of freedom as if the fitted distributions were independent (see Section 3.2.3).
We have also tried fitting an M4 (K = 1) model to artificial datasets constructed with model M2 (K = 0). A good fit should not be possible in that case. The lowest χ 2 values in most of these cases stayed above ∼1000. The fitted parameters present in both models (e.g. x 1 , x 2 ) would often be several standard deviations from their correct values. We consider this result to be an indication of what may be expected if none of the models considered in this work yield an adequate description of the observed data. Tests performed with other models are consistent with the above results.

Models Fitted to the Complete Main Sample
We attempt to reproduce the observed SDSS-FIRST sample by performing simulations to the whole main sample described in Section 2.2. The results of these optimizations are given in Table 1 for all four models. Based on our analysis, the least likely models to fit the observations are M1 and M2, in which the radio luminosity of quasars is independent of their optical luminosity (K = 0). The minimum total χ 2 values reached in our optimization procedure were ∼4000 and ∼2500 for models M1 and M2, respectively. Each of the models has 37  TABLE 1 Best-fit median values and 1σ uncertainties of model parameters for our four models, obtained from fits to the complete main sample of quasars. Value of K is a propery of the models and parameters describe the Gaussian peaks (x 1 and x 2 ), widths (σ 1 and σ 2 ), and their relative normalization (f ). See Section 3.1 for details. The parameters x 1 and σ 1 set the radio-quiet portion of the radio loudness distribution, f is the fraction of radio-loud quasars and x 2 and σ 2 set the radio-loud part. The distributions were computed through Monte Carlo simulations that compensate for the instability of the χ 2 computed using a low number of realizations per simulation (see Section 3.3 for details).
degrees of freedom. For both models, the largest contribution to the total χ 2 comes from the disagreement with the observed redshift distribution. They also fail to correctly reproduce the most populated bins in the radio and optical magnitude distributions. Model M2, however, fits the bright end of the radio magnitude distribution considerably better and this is clearly reflected in the lower χ 2 value. Values of χ 2 generally fluctuate with σ χ 2 < 100 (see Figure 2 and Section 3.2.3), making the difference of 1500 in total χ 2 of the two models highly significant.
Models M3 and M4, both involving a proportionality between the optical and the radio luminosity of quasars (K = 1), were found to provide better fits to the main sample compared to models M1 and M2. Their lowest χ 2 values of ∼2000 and ∼1000, respectively, are significantly lower than for models M1 and M2. Both models can reproduce the observed redshift distribution much more correctly (especially the lowest redshift bin), but the overall match is not a statistically good one since both have only 34 degrees of freedom. The most significant contributors to the total χ 2 are the optical magnitude bins with the highest numbers of quasars, with the bright end of the radio magnitude distribution contributing considerably to the higher χ 2 value for model M3. The non-optimal matching in the radio and optical magnitude distributions results in discrepancies in the R ′ = 0.4(i − t) distribution as well.

The Best-fit Model for the Complete Main Sample
Among the models considered in this paper, the lowest χ 2 values for the complete main sample were achieved for model M4 (K = 1 and a double Gaussian radio loudness distribution). We plot the marginal probability distributions for the parameters of the best-fit model M4 in Figure 4. Note that the expected scatter in χ 2 values between simulations (due to a limited number of realizations) was taken into account prior to constructing the probability density distributions (see Figure 2 and Section 3.2.3 for details). The median values and 1σ uncertainties derived from the marginal probability distributions for all five parameters of the model are given in Table 1. In Figure 5 we plot the t, i, z and R ′ distributions for the observed SDSS-FIRST data (the main sample) and the median values obtained from a simulation (1000 individual realizations) performed with model M4 and the best-fit set of parameters. The error bars on each bin are the minimum and the maximum value occurring in that bin when parameters are shifted randomly within ±1σ from their respective best-fit medians.
In Figure 6 we show the radio loudness distribution for our best-fit model M4, as well as the associated 1σ uncertainties. We compute the probability of bimodality of the best-fit radio loudness distribution by deriving the fraction of the M4 model parameters which result in a bi-  modal distribution. We consider a distribution bimodal if there is at least a slight minimum between the RL and RQ regimes. We infer a probability of ∼20% when we take the 1σ confidence intervals of each parameter into account. Thus, we conclude that the likelihood of bimodality in the radio loudness distribution of our best-fit model is relatively small. Note, however, that based on our simulation tests (Section 3.4) the total χ 2 value of our best-fit model (1000 for 34 degrees of freedom) indicates that it is not a statistically acceptable representation of the data. In order to find a better-fitting model, in the next section we test for possible dependence of the M4 model parameters on redshift and optical luminosity.

Dependence of the Best-fit Model Parameters on
Redshift and Optical Luminosity In order to investigate whether a certain parametrization of dependence on redshift or optical luminosity could be added to our M4 model, we have divided our initial main sample into 4 bins in redshift and 4 bins in apparent magnitude and repeated the optimization procedure for each of them. Since most of these 16 subsamples span fairly narrow ranges in redshift and optical magnitude, the fits were done using only radio magnitude (t) and R ′ = 0.4(i − t) distributions and hence the number of degrees of freedom for this case is 14 per subsample. The choice for the binning was such that every subsample contains ∼6000 optical and ∼500 radio quasars, which is still large enough to perform fits to a distribution divided into 10 bins. The binning is given in Table 2, along with the results for each subsample.
Each of the 16 subsamples was best fitted with model M4, although for the optically faintest subsamples model M3 was almost equally well fit. χ 2 values for individual subsamples range from 10 to 32 for 14 degrees of freedom and the lowest total χ 2 value summed over all 16 subsamples is 291 for 224 degrees of freedom (reduced χ 2 is 1.3). This is much lower than the lowest χ 2 value for the M4 model fitted to the complete main sample, which is 1053 for 34 degrees of freedom. This result further confirms that a simple M4 model, which represents the distribution of the radio-to-optical ratio as a universal function independent of redshift and/or optical luminosity, is not a satisfactory representation of the real quasar sample from SDSS and FIRST.
The results are plotted in Figure 7 using median redshift and absolute optical magnitudes of the subsamples, showing that there are no clear trends in parameters x 1 and σ 1 and that possible trends exist in the remaining three parameters as a combination of dependence on redshift and/or optical luminosity. Results for the optically brightest subsample of each redshift bin were used to plot the continuous change of the radio loudness distribution in Figure 8. With the current results it is not possible to disentangle the two dependencies or to quantify them, but we do discuss the tentative trends in the following section.

Implications from the Parameters of the Best-fit
Model Our best-fit model hints at the possibility that there might exist two distinct populations of quasars, assuming that the double-Gaussian parametrization is an appropriate one. Considering the fraction of radio loudness distributions that show bimodality within 1σ confidence intervals of our best-fit model parameters, the probability that the global radio loudness distribution is bimodal is ∼20%. Therefore, we can neither confirm nor firmly exclude that the radio loudness distribution is bimodal. A more robust result, however, is that the radio loudness distribution of SDSS-FIRST quasars can be described much better with two Gaussians than with a single one. This would imply that RL quasars obey a different relationship between radio and optical luminosity compared to RQ quasars. In this paper we parametrize our models so that there is either no relationship between radio  2 The bin limits (i min , i max , z min , z max ), median optical magnitudes (i med and M i,med ) and median redshifts (z med ) for the 16 subsamples used to examine changes in best-fit parameters of the M4 model with redshift and optical luminosity. Best-fit parameters of the M4 model and their associated 1σ uncertainties are given for each of our 16 subsamples.   Table 2). Colors and symbols mark different redshift bins. The dashed horizontal lines mark best-fit values for the complete main SDSS-FIRST sample and the dotted lines indicate 1σ uncertainties. See Section 5.1 for a discussion of possible trends. and optical luminosities (K = 0 models, M1 and M2), or L 1.4GHz = L i-band × 10 R1 (K = 1 models, M3 and M4) for both types of quasars. The latter parametrization, for which we find a significantly better fit than for the former one, implies that the constant terms of the relationship (parameters x 1 and x 2 , locations of the two Guassian peaks in R 1 ) and its scatter (σ 1 and σ 2 , widths of the Gaussians in R 1 ) are different for the two classes of quasars. In general, the proportionality between the radio and optical luminosity might be different for each of the two quasar classes, but this type of a model was not considered in the work presented here.
Investigating a possible dependence on redshift or optical luminosity in our SDSS-FIRST sample, we have found a statistically good fit when the full sample is divided into 16 subsamples, each fitted independently (reduced χ 2 =291/224=1.3, summed over the 16 subsamples). The parameters describing the RQ Gaussian (x 1 and σ 1 ) were found not to vary significantly, which is not unexpected since this Gaussian is constrained only by its fall-off towards the RL regime. A possible trend with the absolute optical magnitude is suggested for the parameters describing the RL Gaussian (x 2 and σ 2 ) and the radio-loud fraction (f ). The optically bright quasars appear to be better described by an RL Gaussian which is 'louder' (larger x 2 ), a smaller dispersion (σ 2 ) and a lower radio-loud fraction (f ) than the optically faint quasars. There is a slight possibility that the trends are biased by the redshift-luminosity correlation inherent in all fluxlimited surveys, but our method was designed specifically to avoid that kind of bias. Our tests performed on artificial datasets (see Section 3.4) lend confidence that model parameters can be recovered correctly to within statistical uncertainties under a variety of different conditions. Figure 8 shows a possibly significant change of the radio loudness distribution shape implied from our results for quasars with i < 18.5 (the brightest subsample in each redshift bin). This apparent shift of the RL peak and the change in its width could be indicative that the Fig. 8.-Radio loudness distribution as a function of redshift for flux-limited subsambles of quasars with i < 18.5. Note that bimodality seems to have been more prominent at high redshift and that it flattened out at low redshift. The same general trend is detectable in fainter subsamples, but the shape of the distribution does not change in a monotonous manner. radio loudness distribution becomes increasingly bimodal at high redshift. A monotonous change in the radio loudness distribution shape is clearly visible for our subsamples with i < 18.5, indicating that radio-loud quasars are rarer, 'louder' and less spread out in their loudness at high redshift. The same general trend is roughly detectable for subsamples with i > 18.5, but the change of the radio loudness distribution shape is not nearly as clear and continuous as for the i < 18.5 subsamples. We wish to emphasize that we understand selection effects for the i < 18.5 subsample much better than for the i > 18.5 subsample. First, the quasar targeting is complete only for i < 19 candidates. Second, the quasar variability will introduce an RMS scatter of 0.2-0.3 magnitudes and thus blur this supposedly sharp flux limit. Hence, our choice of i < 18.5 binning defines the faintest sample for which the simple SDSS optical selection function is highly reliable. The lack of clear and continuous change in fainter samples may be reflecting the lower reliability of the i > 18.5 objects in the SDSS quasar sample, or indicate that the dependence of the radio loudness distribution on redshift or optical luminosity is not monotonous. We plan to investigate this in the future with more complex model and different statistics (e.g. maximum likelihood), so that more information may be extracted from the existing data.

Comparison to Previous Findings
Previous work on the issue of the bimodality in the radio loudness distribution of quasars were either based on small samples of quasars, or large but less reliable ones. For example, White et al. (2000) and Cirasuolo et al. (2003a,b) used spectroscopically confirmed samples of 636, 141 and 195 radio-detected quasars, respectively, while Ivezić et al. (2002Ivezić et al. ( , 2004b used photometric samples with ∼4400 and ∼10,000 radio-detected quasar candidate sources. In comparison to those results, the results presented here have considerably better statistics and reliability. Our result that a two-component model with a direct proportionality between radio and optical  Table 1 for model M4. The thick dashed line shows the result one would obtain with our sample by not accounting for the incompleteness of FIRST survey near its flux limit. luminosities fits the observed data best is consistent with Cirasuolo et al. (2003b), who have used a method similar to ours on a much smaller heterogeneous sample. Using a stacking analysis of the FIRST data to probe faint radio fluxes, White et al. (2007) also found a strong dependence of the median radio luminosity of quasars on the optical luminosity.
Our results can be compared to some of the previous ones as plotted in Figure 9. The prominent minimum separating the RQ and RL Gaussians observed at R 1 ∼1.3 by Cirasuolo et al. (2003b) and Ivezić et al. (2004b) is likely due to the incompleteness of the FIRST survey at its faint end. In the region around R 1 ∼ 1, FIRST is not more than ∼70% complete for the faintest radio sources, but this was not taken into account in the earlier works. The plot also shows how our result would change in the case where we neglect FIRST incompleteness and assume it is 100% complete down to its flux limit -in that case we would have found a weak bimodality with a minimum at R 1 ∼1. Its relative weakness in comparison to previous results can be understood as a difference in the fraction of optically faint quasars in the sample, which is lower in our case (45% here compared to ∼75 in other cases), and which limits the influence of incompleteness to lower values of R 1 . Note that the flux limit in the radio is the same for all three studies compared in Figure 9, so optically fainter samples introduce incompleteness at a higher value of R 1 . Some of the difference with previous results may also be related to the uncertainties arrising from K-corrections and optical variability.
A small degree of discrepancy exists in the fraction of RL quasars between our and previous results; f = (12 ± 1)% (inferred here) compared to (8 ± 1)% from Ivezić et al. (2004b) and (3 ± 2)% from Cirasuolo et al. (2003b). We have defined the radio-loud quasar fraction Fig. 10.-Dependence of the radio-loud fraction, defined with a threshold in radio loudness and calculated from a best-fit model, as a function of apparent i magnitude (grey symbols and lines). For comparison, the same result from Jiang et al. (2007) is plotted in black. M4 model parameters for each subsample were obtained independently and radio-loud fraction was calculated by integrating the analytic form of the best-fit radio loudness distribution function (e.g. it is not the f parameter of model M4).
as the ratio between the areas under the RL and RQ Gaussians, denoted in our models as parameter f . In other words, the fraction of RL quasars is the fraction of quasars whose radio loudness is determined by the RL Gaussian. Note that with this definition the range of radio loudness between the RL and RQ peaks is shared by both types of quasars. Cirasuolo et al. (2003b) have used that definition as well, but Ivezić et al. (2002Ivezić et al. ( , 2004b defined radio-loud quasars as the ones having R 1 > 1. This definition is more practical and a radio loudness may be unambiguously assigned to each individual quasar. With this definition, our RL fraction is (10 ± 1)%, consistent with the results of Ivezić et al. (2004b). Jiang et al. (2007) have found that the fraction of RL quasars depends both on optical luminosity and redshift, practically independent of the exact definition of an RL quasar. They have found that the fraction of RL quasars, when defined with a threshold in radio loudness is higher for optically bright quasars and at low redshift. Our results tentatively confirm that such trends exist, albeit at low significance (see the lowest panel in Figure 7). With our definition of the RL fraction (the fraction of quasars in the radio-loud Gaussian, f ) the trend with respect to the optical luminosity seems to be just the opposite -f tends to be lower for optically brighter quasars. Note, however, that even if the radio loudness distribution shape changes (e.g. as in our i < 18.5 results shown in Figure 8), the fraction of RL quasars defined with a threshold in radio loudness can remain constant or have an opposite trend than parameter f . In order to check for consistency with Jiang et al. (2007), we have performed additional fits to subsamples derived from the original main sample by dividing it into narrow bins in apparent optical magnitude. As displayed in Figure 10, the fraction of RL quasars calculated from the best-fit M4 models for each subsample, using a threshold in radio loudness, matches the data from Jiang et al. (2007) very well.
In broad agreement with previous work by other authors, our current results imply that the radio loudness distribution did not have the same shape at all times in cosmic history. It might be possible to investigate this in the future using deep radio imaging of samples in thin redshift slices. For example, Kimball et al. (2011) have used deep imaging with EVLA to probe a statistically complete sample of quasars in a narrow redshift bin at 0.2 < z < 0.3. They found that the radio loudness distribution can be well explained with a superposition of two QSO populations (radio-loud, which is AGN-dominated and radio-quiet, dominated by star formation in the host galaxy), and their distribution of radio luminosity appears similar to the best-fit two-component model presented here. Since currently available survey data lacks either depth in the radio or significant sky coverage, the data from the newer generation of radio instrumentation (JVLA, SKA) will be essential in order to fully constrain the radio loudness distribution and ultimately, understand its physical origin.

SUMMARY
In this paper we present Monte Carlo simulations of the optically-selected quasar population fine-tuned to study the radio loudness distribution of quasars. We investigate four different models based on single and double Gaussian distributions and different radio-to-optical luminosity relationships. Our aim was to investigate the long-standing ambiguity in the existence of the intrinsic difference between radio-quiet and radio-loud quasars and in particular, the existence of bimodality in their radio loudness distribution. The sample used here, based on SDSS DR7 Quasar Catalog matched to the FIRST survey, is the largest ever analyzed (8307 radio-detected quasars), uniformly selected and reliable (spectroscopically confirmed), and the method properly accounts for uncertainties in the K-corrections, optical variability and survey incompleteness.
We find that the best-fit model for this SDSS-FIRST quasar sample is a two-component model, but the components overlap so that radio loudness bimodality (i.e. a minimum between the radio loud and quiet portions of the distribution) is not apparent. Statistics of our fits indicate that even our best-fit model does not describe the data optimally, although it is significantly better than any other model we used. Our main result, that the radio loudness distribution of quasars likely consists of at least two components, agrees with earlier findings indicating the existence of two distinct populations of quasars. In the framework of the simple two-component models presented here, we conclude that bimodality is not likely in a sample which, like ours, covers a broad redshift and optical luminosity range, even if it is present in narrower redshift and magnitude bins (e.g. at z > 1.5 and i < 18.5). We also conclude that the distribution is unlikely to be universal, i.e. independent of redshift or optical luminosity.
Investigating possible redshift and optical luminosity dependence of the radio loudness distribution, we have found that a monotonous change of its shape is visible for quasars with i < 18.5. It would imply that at high redshift radio-loud quasars were rarer, on average 'louder' and less spread out in radio loudness. The same general trend is marginally detectable for fainter quasars in the