© ESO 2020

1. Introduction

At the very low and high end of the electromagnetic spectrum, in the radio and gamma-ray regime, the extragalactic sky is dominated by highly beamed jets from active galactic nuclei (AGN, see, e.g., Becker et al. 1995; Ackermann et al. 2015). Blazars are a particularly variable and luminous AGN group and the subject of multiple monitoring campaigns at all accessible wavelengths (e.g., Villata et al. 2008; Ojha et al. 2010; Jorstad & Marscher 2016). Aside from their prominent occurrence at radio and gamma-ray energies, blazars, and especially low-peaked blazars, are notoriously hard to detect in the hard X-rays, just above 10 keV–20 keV. This is largely due to the characteristic spectral minimum between the low-energy synchrotron bump and the high-energy (HE) emission bump at MeV to TeV energies. Also, currently operating hard X-ray instruments such as INTEGRAL and Swift/BAT are usually background-dominated and do not reach the flux sensitivity of soft X-ray observatories by several orders of magnitude (see, e.g., Bottacini et al. 2012; Wang et al. 2016).

Because of the proximity of the keV to the gamma-ray range, multiwavelength studies of blazars rely on this sampling point in order to model the start of the HE bump of a blazar’s spectral energy distribution (SED). Particularly, the modeling of SEDs of the enigmatic group of gamma-ray-faint blazars requires good coverage of the neighboring spectral range. Past studies (e.g., Giommi et al. 2006; Ajello et al. 2009; Draper & Ballantyne 2009) suggest that between soft X-ray and MeV energies, blazars become the dominating contributor to the diffuse cosmic X-ray background (CXB).

Previous all-sky surveys and studies at hard X-rays with the Swift/BAT instrument featured a moderately high significance cutoff at approximately 5σ, effectively hiding a large number of blazars, and consequently leading to small sample sizes (see, e.g., Baumgartner et al. 2013; Krivonos et al. 2015). The 105-month BAT survey catalog (Oh et al. 2018) counts 158 “beamed AGN”, including blazars and flat-spectrum radio quasars (FSRQ), out of 1632 sources, including 114 unclassified AGN and an additional 129 unidentified sources. The recent 4FGL catalog of the gamma-ray all-sky survey from Fermi/LAT (The Fermi-LAT Collaboration 2020) alone contains 1102 BL Lac type and 681 FSRQ type blazars out of a total of 2940 associated blazars and blazar candidates.

In this work, we perform a detailed analysis of X-ray data from the most recent Swift/BAT survey maps for the radio-selected MOJAVE-1 (Monitoring Of Jets in Active galactic nuclei with VLBA Experiments) beamed AGN sample. We report the hard X-ray characteristics in the 20 keV–100 keV band of this statistically complete extragalactic jet sample, which is mostly composed of blazars.

This paper is structured as follows: Sect. 2 gives a short introduction to the MOJAVE program as well as the Swift mission and describes the object samples and data basis we used. The analysis, including the applied methods and assumptions, is presented in Sect. 3. The results, including the statistical properties of the entire sample and distributions for flux and luminosity, are presented in Sect. 4. Several key aspects of this study are discussed in Sect. 5, such as the peculiar hard X-ray properties of individual sources, the potential Fermi/LAT detection of gamma-faint blazars, and the log N–log S distribution of the MOJAVE-1 blazars. A short summary of the main results and conclusions can be found in Sect. 6. The derived data set for all spectral fits for all sources is presented in Appendix A.1.

2. Object samples

2.1. The MOJAVE-1 sample and Swift

The MOJAVE program¹ provides continuous interferometric measurements of the radio-brightest AGN in the northern hemisphere with VLBI (Lister et al. 2009a). We concentrate on the flux-limited MOJAVE-1 sample, counting 135 sources (Lister & Homan 2005). The selection criteria for this sample were: δ ≥ −20°, Galactic latitude $|b|\ge 2\stackrel{°}{.}5$, a total 2 cm (15 GHz) flux density greater than 1.5 Jy at any epoch between 1994.0 and 2004.0, and more than 2 Jy for sources below the celestial equator. This northern AGN sample can be considered statistically complete in terms of 15 GHz radio flux density due to strict selection criteria applied over a large amount of time, as well as the extensive flux density database that was used (Lister et al. 2009a, and references therein). Following the optical classification scheme of Véron-Cetty & Véron (2003), the sources of this sample can be divided into 104 FSRQs, eight radio galaxies, 21 BL Lacs, and two unidentified objects with no known optical counterpart. The source 1219+044 has been re-classified as a NLSY1 galaxy by Yao et al. (2015). We keep the original FSRQ classification (Hovatta et al. 2014) since the intrinsic X-ray luminosity and photon index are notably atypical for the MOJAVE-1 sources classified radio galaxy, and very similar to the FSRQ classification. In general, for the naming of all sources the IAU B1950 coordinate format is applied.

In November 2004, the Neil Gehrels Swift Observatory (Gehrels et al. 2004) was launched by NASA, bringing a multiwavelength observatory in low-Earth orbit that is dedicated to the detection and study of gamma-ray bursts (GRBs). The Swift satellite is equipped with the optical/UV telescope UVOT (Roming et al. 2005, 170 nm–600 nm), the narrow-field X-ray telescope XRT (Burrows et al. 2005, 0.2 keV–10 keV), and the wide-field coded-mask system, the burst alert telescope (BAT, Barthelmy et al. 2005, 14 keV–195 keV). The purpose of BAT, aside from GRB detection, is the monitoring of hard X-ray light curves for a list of transient and variable sources, as well as a continuous blind all-sky survey. The most recent catalog of galactic and extragalactic hard X-ray sources from Swift/BAT for the first 105 months of operation reaches a flux level of 7.24 × 10⁻¹² erg s⁻¹ cm⁻² over 50% of the sky and 8.40 × 10⁻¹² erg s⁻¹ cm⁻² over 90% of the sky at very uniformly distributed exposure times (Oh et al. 2018).

In previous versions of the Swift/BAT all-sky survey catalog, only sources above a Crab-weighted signal-to-noise ratio (S/N) of 4.8σ found in a blind search were included due to the dominance of local noise below this threshold (Tueller et al. 2008; Baumgartner et al. 2013; Oh et al. 2018). Consequently, sources that are expected to have low emission in the hard X-ray range of BAT are likely to be under-represented in the catalogs.

We performed the extraction of spectra and Crab-weighted S/Ns using the Swift/BAT 8-band all-sky mosaic images of 105 months of survey data. Following the procedure described by Baumgartner et al. (2013), the extraction of fluxes of each sample source is performed with an exclusion zone of $40\stackrel{\prime }{.}5$ (15 pixels) around every target using the ftools software batcelldetect. The associated errors are determined by the rms value in the surrounding map area with the radius of $4\stackrel{°}{.}5$ (100 pixels). Crab-weighted S/N values are computed for all source coordinates using Crab-weighted versions of the all-sky images.

2.2. Comparison with other catalogs

We compare the MOJAVE-1 sample with other catalogs in the hard X-ray and gamma-ray regime in terms of detection statistics in order to assert the nature of spectral properties in the hard X-rays, as well as selection biases of other surveys. For this purpose we use the Swift/BAT 105-month source catalog (Oh et al. 2018), the INTEGRAL/IBIS AGN survey catalog (Malizia et al. 2012), the INTEGRAL/IBIS 11-year survey catalog (Krivonos et al. 2015), the Fermi/LAT 4-year AGN catalog (3LAC, Ackermann et al. 2015), and the Fermi/LAT 8-year source catalog (4FGL, The Fermi-LAT Collaboration 2020).

In order to determine the number of common sources in the BAT 105-month source catalog and the MOJAVE-1 sample, the empirical relation of 90% error radius and signal strength is applied as a guideline (Oh et al. 2018),

$$\begin{array}{c}\hfill r_{\mathrm{error}}(\prime )={({\left(\frac{30.5}{S/N}\right)}^2+0.1^2)}^{\frac{1}{2}}.\end{array}$$(1)

Within the source-specific X-ray error radius we find 36 sources that are common between the two samples. We estimate false and missing associations as very unlikely, since there is a significant difference (factor of about 100) in angular separation between the coordinates of BAT counterpart coordinates and MOJAVE source positions, for sources marked associated and unassociated, respectively. The largest separation of an associated source is $0\stackrel{\prime }{.}03$ (2201+315), while the mean separation of unassociated MOJAVE-1 sources to the nearest BAT catalog source is $3\stackrel{°}{.}11$, with a standard deviation of $1\stackrel{°}{.}65$. Based on the angular separation of 2201+315 we calculate the probability of 4 × 10⁻⁸ of a pure chance association between both sub-samples. Objects with moderately high BAT S/N values, larger than approximately 4σ, are registered in both catalogs².

The 36 MOJAVE-1 sources in the 105-month BAT catalog consist of 28 FSRQs, five radio galaxies, and three BL Lacs. The vast majority of the 36 common sources between both catalogs is classified as “beamed AGN” in the BAT 105-month catalog. Only two sources, both radio galaxies, are listed as Seyfert galaxies: 0415+379 (3C 111) and 1957+405 (Cygnus A) according to the classification in the catalog.

We compare the positional information of the MOJAVE-1 sources with another deep hard X-ray survey, the INTEGRAL/IBIS 11-year survey (Krivonos et al. 2015), consisting of all-sky measurements above 100 keV. From 35 AGN that have been detected by INTEGRAL above 4σ in the 100 keV–150 keV band, six sources are shared with the MOJAVE-1 sample: three FSRQs and two radio galaxies with BAT S/N values all larger than 37σ, and the FSRQ 1219+044 with a BAT S/N of 13σ. In general, common sources in the MOJAVE-1 sample and the INTEGRAL survey have the highest BAT S/N values in the MOJAVE-1 source list. The distribution of BAT S/N values of the MOJAVE-1 sample is described in Sect. 4.1 in more detail.

The INTEGRAL/IBIS AGN catalog (Malizia et al. 2012) lists 272 AGN that were observed in the X-ray bands of 2 keV–10 keV and 20 keV–100 keV, of which 57% are classified as Type 1 AGN, including blazars, Seyferts 1–1.5 as well as several sources of mixed classification such as 3C 273 (classified as Sy1/QSO). Fourteen sources are present in both the MOJAVE-1 and the INTEGRAL/IBIS AGN catalog, five of them categorized as radio galaxies in the MOJAVE-1 sample and Seyfert or QSO Type 2 in the INTEGRAL sample. The nine remaining sources are classified as blazars.

The most comprehensive all-sky survey at gamma-ray energies has been performed by the Fermi/LAT instrument, resulting in the 4FGL catalog (50 MeV–1 TeV), which includes 5098 sources above a significance level of 4σ, 2940 of them blazars or blazar candidates. The 4FGL catalog shares 112 sources with the MOJAVE-1 sample, including 86 FSRQs, 21 BL Lacs, three radio galaxies, and two unidentified types at a broad range of BAT S/N values.

Whereas the fraction of sources in the MOJAVE-1 sample with a counterpart in the 4FGL catalog is 83%, the fraction in the beamed AGN sub-sample of the BAT 105-month catalog is only 65%. A significant correlation between radio flux densities and gamma-ray flux measured by Fermi/LAT has been found since its operation (e.g., Ackermann et al. 2011; Mufakharov et al. 2015). The relation of the detection statistic in the LAT survey of X-ray-bright blazars is discussed further in Sect. 4.2. For the sake of comparing gamma-ray detections, we also include data from the previous 3LAC catalog, which shares 100 sources with the MOJAVE-1 sample but also features an integration time frame that is much more compatible with the BAT data (2008–2012 and 2004–2013, respectively). All common sources between the MOJAVE-1 sample and the aforementioned catalogs are marked in Table A.1.

3. Analysis

3.1. Spectral fitting

The full spectral range of the BAT survey comprises 14 keV–195 keV. Due to the background-dominated character of the BAT instrument, especially for sources with low S/N values around 2σ − 3σ, the first spectral bin (14 keV–20 keV) is characterized by very low or even negative count rates. Likewise, the last two bins (100 keV–150 keV, 150 keV–195 keV) feature the very low count rates at high energies. Consequently, we only fit the reduced energy range of 20 keV–100 keV for all sample sources. We model the Swift/BAT spectra, extracted from the 105-month survey maps, with the XSPEC model pegpwlw that is a simple power law with norm K in the aforementioned fixed energy band, expressed as

$$\begin{array}{c}\hfill f(E)=K{E}^{-\mathrm{\Gamma }},\end{array}$$(2)

with the photon index Γ. Spectral fitting of the BAT spectra in the previous object catalogs (Tueller et al. 2008; Baumgartner et al. 2013; Oh et al. 2018) made use of the χ² fitting statistic. Since this fitting statistic assumes Gaussian-distributed data with at least about 30 counts per energy bin, fitting of sources with lower numbers of counts is not feasible. Many of the BAT spectra that are analyzed in this work are low count spectra and, therefore, need to be treated in a different way. Fitting spectra with very low count rates is typically done using the maximum likelihood-based statistic for Poisson-distributed data according to Cash (1979). However, this only accounts for counts that can be described by Poisson statistics with no background, which can usually be modeled separately. BAT spectra are count rate spectra already subtracted by the dominant instrument background. In order to model Poisson-distributed counts of the source itself and a Gaussian-distributed background we use the recommended profile likelihood statistic PGSTAT from the XSPEC statistics appendix³ with zero background counts,

$$\begin{array}{c}\hfill PG=2{\displaystyle \sum _{i=0}^N{t}_{\mathrm{s}}(m_i+f_i)-S_{i}ln(t_{\mathrm{s}}{m}_i+t_{\mathrm{s}}{f}_i)+\frac{1}{2{\sigma }_i^2}{t}_{\mathrm{b}}^2{f}_i^2-S_i(1-ln{S}_i),}\end{array}$$(3)

where

$$\begin{array}{c}\hfill f_i=\frac{-t_{\mathrm{s}}{\sigma }_i^2-t_{\mathrm{b}}^2{m}_i+d_i}{2t_{\mathrm{b}}^2},\end{array}$$(4)

and

$$\begin{array}{c}\hfill d_i=\sqrt{{(t_{\mathrm{s}}{\sigma }_i^2+t_{\mathrm{b}}^2{m}_i)}^2-4t_{\mathrm{b}}^2(t_{\mathrm{s}}{\sigma }_i^2{m}_i-S_i{\sigma }_i^2)},\end{array}$$(5)

with the observed rates S_i, the corresponding uncertainties σ_i, and the predicted rates m_i per channel i. The exposure times for source and background are t_s and t_b, respectively.

A spectral fit is performed for 77 of all 135 sources in the energy range of 20 keV–100 keV (five energy channels). The spectra of all remaining sources possess at least one energy channel with a negative number of counts, which is the result of the background subtraction of these already faint sources. The second logarithmic term in Eq. (3) does not permit a spectral fit in these cases. Instead, the flux is determined using a simulated spectrum with the equal number of counts of the real spectrum and a photon index that has been frozen. The applied photon index is derived as a weighted mean of all photon indices of the fitted spectra in the sample of the same object class (FSRQ, BL Lac, or radio galaxy). For the group of unidentified sources the index of FSRQs is applied, which is the largest sub-group in the sample. Since sources of very high significance and low uncertainty heavily dominate the weighted mean only sources below a significance of 20σ were used in calculating the mean. The template indices are listed in Table 1.

Error bars corresponding to 90% uncertainty ranges for flux and photon indices are calculated for each source via a Monte Carlo approach: the derived values for X-ray flux, photon index, and the Swift/BAT survey response file⁴ are used to simulate a new ideal BAT spectrum. A new and randomized spectrum is then created using the ftools program batphasimerr. Both the flux and photon index of a new power-law fit are then determined. This procedure is repeated 2000 times for every source. The resulting distributions of flux and photon index are Gaussian-shaped. The 90% error ranges are then determined using a Gaussian fit function.

The uncertainties of the spectral parameters for sources with negative spectral counts are derived from the count rate error bars of every energy channel as well as the uncertainty of each template photon index. If the flux of a source in any case is less than the error bar of the flux corresponding to 3σ, it is considered an upper limit.

Spectral contamination of three nearby X-ray sources also necessitates the calculation of an upper limit flux value. The ratio of contamination rate (estimated count rate of all nearby sources from batcelldetect) and count rate is defined as the contamination ratio. We adopt the convention by Baumgartner et al. (2013), where a source is considered confused or contaminated if the contamination fraction surpasses a value of 0.02.

3.2. Source counts and CXB contribution

The number-count diagram, or log N–log S, conveys a number of statistical properties of the flux emission of any given sample as well as possible selection effects. The cumulative log N–log S distribution with the number of sources per square degree with fluxes greater than F_j of source j is defined as:

$$\begin{array}{c}\hfill N(>F_j)=\frac{{\mathrm{\Omega }}_{\mathrm{survey}}}{{\mathrm{\Omega }}_{\mathrm{sky}}}{\displaystyle \sum _{i=0}^N\frac{1}{{\mathrm{\Omega }}_j},}\end{array}$$(6)

with the sky area Ω_j corresponding to the flux F_j, the total sky area Ω_sky, and the total sky area covered by the MOJAVE-1 survey Ω_survey⁵. Based on Poisson statistics the error bars are calculated via:

$$\begin{array}{c}\hfill N_{\mathrm{err}}(>F_j)=\frac{N(>F_j)}{\sqrt{N}}·\end{array}$$(7)

We use the log N–log S statistic to compare the sample source density distribution in space to the Euclidean distribution of flat space, that is described by N = AF^−α with α = 1.5. Any divergence from this shape indicates an evolution of the source flux in this sample and/or a source detection bias (see, e.g, Longair 1966; Mateos et al. 2008).

From the fitted log N–log S diagram we also derive the contribution of blazars to the CXB that has been studied in a number of cases in the past (e.g. Gilli et al. 2007; Beckmann et al. 2006; Ajello et al. 2009). Depending on the specific energy range, the CXB can be resolved into different AGN sub-classes. Studies of the blazar contribution above 10 keV suffer from low sample sizes and are often limited by the low flux of these sources (see, e.g., Ajello et al. 2008a, 2009). Other studies rely on the extrapolation of the contribution to the CXB from lower energies (see, e.g., Giommi et al. 2006). Results therefore vary considerably with the used blazar sample and analytic approach. A contribution of approximately 1.5% in the 14 keV–170 keV range and 10%–20% for the 15 keV–55 keV range has been reported by Ajello et al. (2008a) and Ajello et al. (2009), respectively. The role of blazars in the creation of the CXB above its maximum flux at approximately 30 keV (Gruber et al. 1999) remains subject of discussion. Previous studies have addressed samples of X-ray-bright blazars usually detected with S/N values larger than 5σ (e.g., Sazonov et al. 2007; Ajello et al. 2009; Bottacini et al. 2012). Here, we incorporate a well studied blazar sample, defined by long-term radio surveillance and a broad X-ray S/N distribution. Based on the known blazar positions of the MOJAVE-1 sample we measure their flux distribution to the CXB in the energy range of 20 keV–100 keV.

In order to determine the flux $F_{\text{CXB}}^{{\text{deg}}^{\text{2}}}$ per sky area of the X-ray background itself in this range we integrate the CXB spectrum given by Eq. (5) of Ajello et al. (2008b) and obtain:

$$\begin{array}{c}\hfill F_{\mathrm{CXB},20-100\mathrm{keV}}^{{\mathrm{deg}}^2}=3.14\times {10}^{-11}\mathrm{erg}{\mathrm{s}}^{-1}{\mathrm{cm}}^{-2}{\mathrm{deg}}^{-2}.\end{array}$$(8)

To compare the results of this analysis with the flux data of the beamed AGN class in the most recent 105-month BAT source catalog we also calculate the X-ray background flux in the full BAT range of 14 keV–195 keV, and obtain:

$$\begin{array}{c}\hfill F_{\mathrm{CXB},14-195\mathrm{keV}}^{{\mathrm{deg}}^2}=4.61\times {10}^{-11}\mathrm{erg}{\mathrm{s}}^{-1}{\mathrm{cm}}^{-2}{\mathrm{deg}}^{-2}.\end{array}$$(9)

The contributing flux of a given source sample can be estimated by integrating the differential log N–log S distribution, that is, the derivative of the function N(>F), multiplied with the source flux (see, e.g., Giommi et al. 2006),

$$\begin{array}{c}\hfill F_{\mathrm{contrib}}^{{\mathrm{deg}}^2}={\displaystyle {\int }_{F_{\min }}^{F_{\max }}\frac{\mathrm{d}N}{\mathrm{d}F}F\mathrm{d}F}\end{array}$$(10)

with F_min and F_max usually describing the lower limit and upper limits of source fluxes in the sample. The integral represents the contribution of flux of all sources (of this population in the sky) brighter than the minimum flux F_min. This integration limit is usually equal to the smallest source flux in the sample, or can be extrapolated to lower fluxes, albeit under the assumption of a constant slope of the log N–log S distribution.

In order to determine the log N–log S distribution, we compute the empirical sky survey area as a function of the minimal flux that corresponds to a detection. We extract the function from the BAT 105-month survey maps for a minimal flux corresponding to a detection at 1σ, 3σ, and 5σ in the 20 keV–100 keV band (graphed in Fig. 1). For the calculation of the log N–log S distribution of the MOJAVE-1 sample in the BAT band we choose the 1σ curve that has been extracted from the survey maps. Although a measured signal equivalent to 1σ above the background can hardly be called a detection this criterion only applies to blind surveys, unlike the procedure described here: the position of every source in the sample is already known and, therefore, the effective survey sky area at a given flux must be higher in any case.

Fig. 1. Swift/BAT 105 month survey sky coverage as the function of minimal flux corresponding to a detection at a level of 1, 3, and 5σ in the band 20 keV–100 keV, from the left to the right. The dashed line indicates the full sky.

3.3. X-ray and radio luminosity function

The shape of the log N–log S diagram of any given source sample is related to the distribution of the sources in space. However, the data are possibly influenced by a number of factors such as selection effects or an intrinsic evolution of the emission in the observed band. In order to understand the shape of the log N–log S distribution of the MOJAVE-1 blazar sample and possibly responsible biases and intrinsic flux evolution we calculate the hard X-ray and radio luminosity functions (LF). The differential luminosity function reveals the density of the sample sources per unit redshift and luminosity. Our approach closely follows the methods implemented by Ajello et al. (2009), Ebrero et al. (2009), and Miyaji et al. (2015). The goal of this part of the analysis is to compare models of positive evolution of luminosity with a no-evolution scenario for the X-ray and radio bands.

The binned LF is determined by the number of sources N in a logarithmic luminosity bin L_bin, min − L_bin, max divided by the co-moving volume dV/dz (see, e.g., Hogg 1999), integrated over the total luminosity range of the sample L_min − L_max and the redshift range z_min − z_max:

$$\begin{array}{c}\hfill \frac{\mathrm{d}\mathrm{\Phi }}{\mathrm{d}logL}=\frac{N}{{\int }_{L_{min}}^{L_{max}}{\int }_{z_{min}}^{z_{max}}\frac{\mathrm{d}V}{\mathrm{d}z}\mathrm{d}z\mathrm{d}logL}·\end{array}$$(11)

The double integral in Eq. (11) is solved using a numerical approach: the co-moving volume is calculated on a discrete log L–z grid with the size of 50 by 50 steps. Each volume element is then multiplied by the step size in both log L and z and summed up.

We also fit an analytic model of the LF directly to the data using a maximum likelihood method. A single and a double power-law model of the differential LF are fitted to the data. The present day LFs (z = 0) are expressed by:

$$\begin{array}{c}\hfill \frac{\mathrm{d}\mathrm{\Phi }}{\mathrm{d}logL}=A{[{\left(\frac{L}{L_{\ast }}\right)}^{{\gamma }_1}+{\left(\frac{L}{L_{\ast }}\right)}^{{\gamma }_2}]}^{-1},\end{array}$$(12)

with the normalization A, along with power-law indices γ₁ and γ₂ for a double power law. The parameter L_* functions as a break luminosity for the double power-law model, that is also fitted. For single power-law models, L_* is fixed and the second term in Eq. (12) is 0. The evolution of the LF with z is given by the term e with the two parameters k and g:

$$\begin{array}{c}\hfill e={(1+z)}^{(k+gz)}.\end{array}$$(13)

It is applied for two different scenarios: pure luminosity evolution (PLE) and pure density evolution (PDE):

$$\frac{\text{d}\Phi }{\text{d}\log L}(L,z){|}_{\text{PLE}}=\frac{\text{d}\Phi }{\text{d}\log L}(L/e,z=0),$$(14)

$$\frac{\text{d}\Phi }{\text{d}\log L}(L,z){|}_{\text{PDE}}=\frac{\text{d}\Phi }{\text{d}\log L}(L,z=0)\cdot e.$$(15)

We can further divide the models by setting the evolutionary parameter g = 0 for one set of PLE and PDE models and allow the parameter to vary for others, named PLEg and PDEg. The maximum likelihood algorithm that is utilized, following, for instance, Miyaji et al. (2015), minimizes the expression of:

$$\begin{array}{c}\hfill S=-2{\displaystyle \sum _{i=0}^{N}log\frac{\mathrm{d}\mathrm{\Phi }}{\mathrm{d}logL}{|}_i-log{\int }_L{\int }_z\frac{\mathrm{d}\mathrm{\Phi }}{\mathrm{d}logL}\frac{\mathrm{d}V}{\mathrm{d}z}\mathrm{d}z\mathrm{d}logL,}\end{array}$$(16)

for all sample sources i. For the minimization of the function and the calculation of the best fit parameters and errors we use the MINUIT software package (James 1994). In order to determine a best-fit model, we apply the Akaike information criterion (AIC, Akaike 1973):

$$\begin{array}{c}\hfill \mathrm{AIC}=2k-2log\mathcal{L},\end{array}$$(17)

with the number of free parameters k and the MINUIT return value ℒ (maximum likelihood value). A minimal AIC indicates the best fit LF model. As outlined in, for example, Burnham & Anderson (2004) the difference of AIC values Δ_j = AIC_j − AIC_min of model j to the model with the lowest AIC value determines if the model is equally or less probable. For practical purposes, the likelihood, p_j, of model, j, can be calculated. It expresses the relative probability, compared to AIC_min/the best fit, that model j minimizes the AIC:

$$\begin{array}{c}\hfill p_j=e^{-{\mathrm{\Delta }}_j/2}.\end{array}$$(18)

4. Results

4.1. Signal-to-noise ratio and BAT source detection

Other than the past BAT blind survey catalogs with a significance threshold of 4.8σ, we derive the significance values and spectra at known source positions, independent of signal strength. For each source in the MOJAVE-1 sample, the S/N value has been extracted from the Crab-weighted 105-month survey maps in the full energy range of 14 keV–195 keV. In Fig. 2 we present the distribution of the BAT S/N values for all source types in the sample. Hatched areas represent sources for which only upper limits for the X-ray fluxes could be determined. The highest value at about 192σ belongs to the FSRQ 1226+023 (3C 273) and has been omitted in the histogram. 27 sources have S/N values between 4.8σ and 20σ, and can be considered clear detections, according to the blind survey conventions of Oh et al. (2018), Baumgartner et al. (2013), and Tueller et al. (2008). A tail of very X-ray bright and highly significant sources (three FSRQs and four radio galaxies) shows higher values still. The brightest source in the BL Lac category is BL Lac itself (2200+420) with a S/N of 16.8σ. All remaining BL Lac objects in the sample have S/N values under 5.4σ. The sub-set of radio galaxies appears relatively bright compared to the other source types, which is due to the close proximity of these sources. The average redshift of the radio galaxies in the sample is 0.06, whereas BL Lacs and FSRQs are characterized by noticeably higher redshifts of 0.39 and 1.16, respectively.

Fig. 2. BAT S/N distribution of MOJAVE-1 sample for different source classes. The brightest source 3C 273 at 192σ has been omitted for better readability.

The majority of 101 sources of the sample lies below the 105-month catalog threshold of 4.8σ, although with a clear offset from 0σ towards positive S/N values. Negative S/N values are viable since the survey maps that we used are already subtracted with a dominant instrument background count rate value. Hence, areas with no bright source emission can be over-subtracted, leading to negative S/N values. Since BAT is a background-dominated coded mask instrument, low signal measurements from source positions might also be of random nature. In order to estimate the number of signals from true sources we compare a pure noise distribution with the S/N distribution of our sample. Figure 3, top panel, shows the S/N characteristic of 1000 pointings over the entire sky which have at least a 100 pixel ($4\stackrel{°}{.}5$) distance from any known source brighter than 4.8σ in the 105-month survey maps. The resulting distribution is centered around 0σ and ranges from about −3.3σ–3.5σ. It can be described very well with a Gaussian function that is centered at 0σ with a width of 1σ, as indicated by the red line (reduced Cash statistic of 0.86 for 19 degrees of freedom). The 99.73% confidence intervals (3σ) of the fit function are shown by the red shaded area.

Fig. 3. Top panel: distribution of Swift/BAT S/N values of 1000 random blank sky positions. The red line indicates a Gaussian fit with center of 0σ and width 1σ. The red shaded region indicates 3σ uncertainty ranges for the norm of the fitted function. Middle panel: distribution of Swift/BAT S/N values of the MOJAVE-1 sample. The red line indicates a Gaussian fit with center of 0σ and width 1σ to the distribution with values smaller than 0σ. The red shaded region indicates 3σ error ranges. The 4.8σ threshold of the 105-month BAT survey is shown by the dashed line. Bottom panel: resulting histogram after subtracting the MOJAVE-1 S/N distribution by the distribution described by the Gaussian fit. Arrows indicate the cutoff at 15σ (for better readability).

We want to stress the difference in the idea of “noise” in terms of signal strength and noise contribution to the signal distribution of a source sample. The former relates to the S/N of the count rate of a central pixel of a source’s coordinates in relation to the mean instrumental background. The latter describes the contribution of low-significance sources (low or even negative S/N values) to the overall S/N distribution of a sample.

A direct comparison of the top and middle panel of Fig. 3 clearly demonstrates that the MOJAVE-1 S/N characteristic is significantly different from random noise⁶. In order to determine the fraction of a source sample whose hard X-ray emission is not compatible with random map noise we subtract a Gaussian model fit from the S/N distribution of the sample. The Gaussian, with a center at 0σ and a variance of 1σ, is modeled using a fit to all negative S/N values in the histogram. The model function and its set parameters for center and width follow the measured pixel significance distribution by Tueller et al. (2008). We apply this method to the measured MOJAVE-1 sample S/N distribution.

The S/N distribution of the MOJAVE-1 sample is shown in the middle panel of Fig. 3. The errors are purely statistical. We fit a Gaussian function with a center of 0 and width of 1σ to the negative part of the distribution, which describes the map noise contribution. The source distribution is then subtracted by the values indicated by the Gaussian fit function. In this way we obtain the number of sources not compatible with random map noise. It remains a number of $121.{94}_{-11.16}^{+4.39}$ sources, representing 82.06%–93.58% of the MOJAVE-1 sources, which is the percentage of the sample that is emitting hard X-rays with respect to the achievable flux limit of the BAT instrument after 105 months of observations.

We note that this global consideration does not yield an immediate statement regarding which individual sources of the low-S/N sources are real detections (i.e., belonging to the subtracted sample within the error margins). However, within a given significance bin, we can determine the probability that a given source is a real detection by comparing the number of expected signals through background fluctuations to the actual number of counts in this bin. For example, in the (0.5σ–1σ) bin, we expect up to approximately five counts from the fitted background distribution (within the 3σ uncertainty range of the fitted Gaussian; see Fig. 3, middle panel). We find thirteen sources in this bin. Each individual of these thirteen sources thus has a probability of less than 5/13 of being associated with a background fluctuation. Within the higher significance bins, this chance-fluctuation probability drops quickly. For example, in the (2σ–2.5σ) bin, we expect only up to 0.5 counts but detect fourteen sources. Each of those thus has a probability of larger than 96% of not being associated with a background fluctuation.

The most recent Swift/BAT catalog (Oh et al. 2018) was only able to identify 36 out of 135 radio-selected MOJAVE-1 AGN, which are for the vast majority blazars. We show that a lower S/N threshold for BAT detections leads to a large increase in number of detected hard X-ray emitting blazars.

4.2. Spectral shape – photon index Γ

4.2.1. General description

The BAT spectra of 77 out of all 135 sources were modeled with power law continua. This includes 46 sources below the 105-month catalog 4.8σ threshold. Spectra with at least one negative count rate bin (58 sources) due to low signal strength have been processed using template spectra (see Sect. 3.1). No photon index has been derived for the three sources 0917+624, 1502+106, and 1928–179, which have sufficient count rates for fitting, but are spectrally contaminated by other nearby sources.

In a canonical spectral energy distribution (SED) of a blazar in a νF_ν plot the photon index of Γ = 2 indicates a flat spectrum. Smaller values correspond to a rising spectrum, typically the left flank of the HE bump at hard X-rays. The majority of all measured blazar photon indices in the sample are smaller than 2, which locates the BAT band at the rising part of the HE bump. Radio galaxies show a different behavior and concentrate around Γ = 2 with the largest value at 3.2 (0316+414/3C 84).

The distribution of derived photon indices Γ is shown in Fig. 4. FSRQs form a quasi-Gaussian distribution around 1.6, with values ranging from about 0.5 to 2.7. The distribution for the low number of BL Lacs shows a broader and less peaked shape, ranging from 0.8 to 3. In the plot the fainter half of the sources (S/N <  4σ) is also marked. The brighter half of the shown source sample (omitting 3C 84) is characterized by a smaller spread compared to the fainter sources, approximately 0.8–2.1 vs. 0.5–3.0. Sources that are less significant show large uncertainties of the photon index, approximately 0.4 for brighter sources and up to 1.7 for the faintest of the fitted spectra. Consequently, the distribution of the photon index becomes noticeably broader.

Fig. 4. Distribution of photon index Γ for all fitted MOJAVE-1 sources. Hatched boxes indicate sources with S/N values below 4σ (fainter 50% of the sample).

4.2.2. Fermi/LAT detections

The slope of the hard X-ray spectrum pinpoints in most cases the rising part of the HE emission bump and therefore holds information about the behavior of the SED at higher energies. Figure 5 displays the distribution of the photon index for all fitted MOJAVE-1 blazars, and compares gamma-bright sources which are listed in the 3LAC catalog (top panel) and gamma-faint sources (middle), respectively. For this part we choose to incorporate the Fermi/LAT detection statistics from the 4 year catalog (3FGL/3LAC), since the observation time slot of 2008–2012 is close to the BAT 105-month survey time range. We compare this with the later 4FGL Fermi/LAT source catalog at the end of this section. The MOJAVE-1 blazars in the 3LAC catalog (sub-sample A) have noticeably harder photon indices than gamma-faint sources (sub-sample B). We perform a two-sample KS-test for the two distributions. The null hypothesis of both data sets drawn from a common distribution can be rejected at a level α = 0.008 with a test statistic of D_sample = 0.497. Additionally, we apply a one-sided Z-test, comparing the difference of means of A and B, depending on the sample size. The resulting p-value of 0.0028 also strongly implies a distinct difference of distribution of values.

Fig. 5. Distribution of BAT photon index Γ for all fitted MOJAVE-1 blazars, detected and not detected by Fermi/LAT, respectively. Bottom panel: added Gaussian distributions for each photon index. The sharp peak at 1.75 stems from the brightest blazar in the sample, 3C 273.

Even so, the tests ignore the partially substantial uncertainties of the photon indices. In order to take into account the (Gaussian-distributed) uncertainties for every source we model the probability distribution of each index as a Gaussian curve with a width corresponding to its uncertainty range. Then, the distributions of probability for every source in each data set A and B are added up. The resulting distributions are shown in the bottom panel of Fig. 5. Since the area of every individual Gaussian curve is equal, bright sources with small uncertainties translate into sharp peaks, whereas sources with large uncertainties equal low and broadened curves. The distributions are plotted as fine histograms on a grid of 1000 bins, normalized to their respective sample size. From the summed distribution of all sources (black line) 10⁵ times two test samples A_simul and B_simul with the respective sample size of A and B are randomly drawn and each time a KS-test is performed. The probability of a certain photon index value to be drawn is determined by the height of the probability function/summed distribution at that value. The goal of the analysis is to find the number of randomly chosen sample pairs, which give a KS-test statistic D that is equal or higher than the test statistic D_sample of the original samples A and B. This percentage indicates how likely it is that the photon index distributions of A and B are distinct form each other.

The fraction of A_simul and B_simul with a KS-test statistic D equal or greater than the test statistic D_sample of A and B is 0.9%. This is equivalent to a significance of 2.61σ for a Gaussian distribution of test statistic values. Thus, sub-samples of Fermi/LAT-detected and non-detected blazars in the MOJAVE-1 sample depict two significantly different groups of sources in terms of spectral shape in the BAT energy range: Fermi/LAT-detected sources are harder, while non-detected sources tend to be softer.

Including the remaining two Fermi/LAT-detected and four non-detected radio galaxies in the MOJAVE-1 sample gives an even more distinctive result. We obtain a sample test statistic of D_sample = 0.55, equal to a rejection level of the null hypothesis of α = 0.002. The fraction of KS-tests with D >  D_sample is 0.05%. This effect however is largely due to the photon index of the radio galaxies near Γ = 2 in sub-sample B, emphasizing the difference of both distributions.

A difference in the photon index or spectral slope in an energy band at the rising part of the SED’s HE bump would correlate with the position of the SED along the frequency axis. In Fig. 6 we show the photon indices of both Fermi/LAT-detected and non-detected sub-samples against the HE emission peak frequencies, taken from Chang (2010) and transformed to the rest frame. Whereas Fermi/LAT-detected MOJAVE-1 blazars tend to concentrate around HE bump peak frequencies $log{\nu }_{\mathrm{HE}}^{\mathrm{peak}}$ of 22.5–23.5 non-detected blazars span a range of approximately 19–24. Taking into account the previous result of 3LAC non-detected sources having noticeably harder BAT photon indices than non-detected sources (see testing procedure above), this indicates a possible correlation of spectral shape in the BAT 20 keV–100 keV band and the probability for a detection at GeV energies. Additionally, the sources which have a 4FGL catalog entry are plotted using filled circles in Fig. 6. The two sources 1458+718 (3C 309.1) and 2145+067, previously undetected by Fermi/LAT, are now sharing the same space in the diagram with the majority of detected sources at high HE peak frequencies.

Fig. 6. SED HE peak frequency of MOJAVE-1 blazars (rest frame) against BAT photon index Γ.

4.3. Hard X-ray flux

The histogram of hard X-ray flux measurements is presented in Fig. 7. We show both the 77 fitted sources and the 58 remaining sources. The vast majority of 127 of 135 sources lies below 30 × 10⁻¹² erg s⁻¹ cm⁻² with four FSRQs and four radio galaxies having higher significance values. Most of the sample, especially BL Lacs and unidentified sources, is characterized by low hard X-ray flux, as shown by the median of 4 × 10⁻¹² erg s⁻¹ cm⁻², also taking into account upper limits. Following the flux calculation described in Sect. 3, we obtain a number of 59 upper limits, including three sources that have been fitted, but are also designated upper limits because of spectral contamination of nearby X-ray sources.

Fig. 7. Distribution of hard X-ray flux of the MOJAVE-1 sample.

Figure 8 shows the relation between X-ray flux and S/N. Seven data points at negative S/N values as well as source 1504–166 (PKS 1504–167) at 0.01σ have been omitted. The 105-month catalog threshold value of 4.8σ is indicated by the dashed line. Beginning at about 2σ, the slope of the narrow distribution in Fig. 8 is about unity, making the relation of flux and S/N linear, as expected. Below about 2σ the flux value saturates around (3−4)×10⁻¹²erg s⁻¹ cm⁻² at a point where almost all flux values are described as upper limits.

Fig. 8. Relation of hard X-ray flux to S/N. The dashed line represents the S/N value of 4.8σ. The source 1504–166 at 0.01σ has been omitted for better readability.

4.4. Hard X-ray luminosity

We calculate the hard X-ray luminosities including the K-correction,

$$\begin{array}{c}\hfill L=\frac{1}{{(1+z)}^{2-\mathrm{\Gamma }}}4\pi d_{\mathrm{L}}^{2}F,\end{array}$$(19)

where d_L is the luminosity distance, and assuming the cosmological parameters: H₀ = 70.0 km s⁻¹ Mpc⁻¹, Ω_M = 0.30 and Ω_λ = 0.70. The resulting distribution of the hard X-ray luminosity is presented in Fig. 9. Four sources with missing redshift information (BL Lacs and unidentified objects) were not included. This quasi-Gaussian distribution centers around (10⁴⁶ − 10⁴⁷) erg s⁻¹, with all BL Lacs and radio galaxies below 10⁴⁶ erg s⁻¹. BL Lacs have minimum values of about 10⁴³ erg s⁻¹. The brightest sources are all FSRQs with luminosities up to 10⁴⁸ erg s⁻¹ (0836+710) and upper limits up to 5 × 10⁴⁷ erg s⁻¹. Radio galaxies possess lower luminosities around 10⁴¹ erg s⁻¹–10⁴⁵ erg s⁻¹.

Fig. 9. K-corrected X-ray luminosity distribution for the energy range of (20–100) keV of the MOJAVE-1 sample.

The two different blazar classifications, FSRQ and BL Lac, clearly do not follow a common distribution. FSRQs exhibit a pronounced maximum, which is more than an order of a magnitude higher compared to BL Lacs. A KS-test for both data sets reveals a test statistic of D = 0.88. We can reject the null hypothesis that the distributions are equal at a level of α = 3.2 × 10⁻⁶.

In Sect. 4.2 we have presented data that suggests a correlation of spectral shape in the hard X-ray regime and the probability of a detection at gamma-ray energies. In order to test for a difference of the intrinsic luminosity of both Fermi/LAT-detected and non-detected sub-samples, we apply a KS-test to both data sets. We obtain a KS test statistic of D = 0.22. The null hypothesis that both luminosity distributions of Fermi/LAT-detected and non-detected MOJAVE-1 blazars are the same can be rejected at a level of α = 0.35. Both distributions are not significantly different from each other in the BAT band.

4.5. Source count statistics

4.5.1. log N–log S distribution

Figure 10 (left panel) shows the log N–log S distribution of the BAT flux values of the MOJAVE-1 blazars (open squares), excluding all upper limit sources, leaving 70 sources. The distribution follows a power law, which is significantly flatter than −1.5 for all fluxes F ≳ 2.2 × 10⁻¹² erg s⁻¹ cm⁻². We fit the data using a least squares approach and obtain the parameters for the best fit with A_Moj, BAT = (8.25 ± 0.57)×10⁻³ deg⁻² and α_Moj, BAT = 1.13 ± 0.04 for the fitting function AF^−α, where F is the flux in 10⁻¹² erg s⁻¹ cm⁻². Below a flux of approximately 3 × 10⁻¹² erg s⁻¹ cm⁻², the distribution of the data points is saturated.

Fig. 10. Left: cumulative log N–log S distribution of the MOJAVE-1 sample excluding radio galaxies, omitting all upper limits (flux for 20–100 keV, blue squares) and log N–log S for the beamed AGN sub-sample of the BAT 105-month source catalog (flux for 14 keV–195 keV, red circles) and the same sub-sample excluding all MOJAVE-1 sources (14 keV–195 keV, purple x). Right: log N–log S for 15 GHz flux density VLBI measurements from Lister et al. (2015) of all blazars in the MOJAVE-1 sample (black circles) and MOJAVE-1 blazars from left graph (blue squares, norm with a factor of 0.1 for better readability). The black lines represent power-law fits, the gray area shows the corresponding error range. The radio data have been fit for all fluxes higher than 2.5 Jy. The dashed lines indicate a slope of −1.5.

The faintest source in the plot (1417+385) is not compatible with the power-law fit within the error margins. The fitted power law at the lowest flux in the sample amounts to a cumulative source number of (3.33 ± 0.22)×10⁻³ deg⁻², or 88 ± 6 sources in the survey part of the sky. The difference to the 70 sources of the fitted sample likely stems from the missing contribution of at least some of the upper limits which have not been included in the plot.

Within the uncertainty range of α_Moj, BAT the log N–log S distribution for the hard X-ray data is not compatible with a regular Euclidean distribution. For comparison we also show the log N–log S distribution of the beamed AGN sub-sample of the 105-month Swift/BAT source catalog (Oh et al. 2018), indicated by open circles, and for flux values in the complete BAT range of 14 keV–195 keV. The BAT catalog is compiled purely on the basis of BAT significance and does not suffer from any selection effects that might be introduced by a catalog at a different wavelength. The original sub-sample of objects designated beamed AGN includes 158 sources. We exclude the brightest source, Centaurus A, which is more commonly classified as a radio galaxy (see, e.g., Steinle 2006, for a review). Also, all sources with a significance smaller than 5σ are excluded, leaving 143 beamed AGN. A power-law fit of the distribution yields a slope of α_cat, BAT = 1.24 ± 0.02. A fit to the log N–log S distribution of the entire BAT 105-month AGN sample (955 sources fitted, not shown) gives a very similar power-law index of α_BATAGN = 1.27 ± 0.01. The slope of the log N–log S distribution of the MOJAVE-1 hard X-ray emission is thus not compatible with the most recent BAT source catalog and the beamed AGN sub-sample (within a 2σ uncertainty range of both results). However, the close proximity of the determined power-law indices (1σ errors) and the significantly greater distance of all indices to the expected Euclidean slope of α = 1.5 suggests a similar underlying reason for this behavior of blazar samples in the hard X-ray domain. In terms of sky area density the BAT catalog sample shows approximately 1.6 times the density of beamed AGN compared to the MOJAVE-1 blazars (non upper limits) regarding the X-ray-faintest source in both distributions. Fitting the log N–log S distribution of the BAT catalog beamed AGN sample excluding all MOJAVE-1 sources results in α_{cat, BAT, noMoj} = 1.38 ± 0.03, bringing the slope closer to the Euclidean case.

Additionally, we fit the MOJAVE-1 log N–log S distribution of available VLBI radio flux density data from Lister et al. (2015), using median flux density values. The observation time of the measurements spans 2008–2012, reasonably close to the BAT survey time. The right panel of Fig. 10 displays the log N–log S distribution of the MOJAVE-1 blazars in the 15 GHz band from Lister et al. (2015), where the black circles denote the full sample of MOJAVE-1 blazars (122 sources, all with known redshift). The blue squares indicate the same sub-set of X-ray bright blazars as in the left graph. Both distributions are plotted and fitted in order to determine possible selection effects in the log N–log S distribution of the BAT data. Fitting the complete flux range can introduce a new data bias, however. Because the MOJAVE-1 sample is just a sub-sample of the compiled flux-limited compact object radio samples of Lister et al. (2015) we choose a cutoff at 2.5 Jy, below which all data are ignored. The vast majority of additional sources to the MOJAVE-1 sample in both larger samples are located below this cutoff. By ignoring all data below the cutoff we avoid using a sample with a flux distribution that is not representative of a flux-limited survey. The normalizations and power-law slopes of all samples that we used are listed in Table 2.

The power-law fits of the radio flux densities larger than 2.5 Jy (122 source and 70 source sub-sample) are described by $a_{\mathrm{Moj},15\mathrm{GHz}}^{122}=1.{74}_{-0.12}^{+0.13}$, and $a_{\mathrm{Moj},15\mathrm{GHz}}^{70}=1.{61}_{-0.14}^{+0.15}$, respectively. The smaller sample of 70 sources lies well within the Euclidean distribution, and is not compatible with the slope of the BAT data within their 1σ uncertainty range. The fit of the larger sample of 122 sources in the radio band produces a slope that is even larger than for a Euclidean distribution. Since the difference in slopes from radio to X-ray band is likely not directly due to a selection effect, we investigate the influence of variability and intrinsic evolution of the emitted X-ray flux in the following sections. The influence of different contributions to the slopes from FRSQs and BL Lacs individually cannot be examined in detail because of the low number of only four BL Lacs in the reduced sample above 2.5 Jy.

To test for possible biases from selection effects of the MOJAVE-1 sample, we analyze the extended MOJAVE-1.5 sample (Lister et al. 2013) and the corresponding log N–log S distribution for radio and BAT fluxes. Contrary to the main sample of this study all sources with a minimum 15 GHz flux density of 1.5 Jy and with δ ≥ −30° and the Galactic plane are included. However, the source confusion problem is not trivial in the hard X-ray band because of the large instrument PSF and high source density, especially within the Galactic plane.

The larger fitted sample of blazars in the radio band reveals a steep slope of $a_{\mathrm{Moj}1.5,15\mathrm{GHz}}^{170}=1.{72}_{-0.10}^{+0.11}$, which is compatible with the results of the MOJAVE-1/Lister data set. As expected, the normalization of the fitted power law for the BAT fluxes is somewhat higher because of the small number of additional sources. Only 22 newly added blazars from the MOJAVE-1.5 sample were fitted because of the high number of upper limits around F = 3 × 10⁻¹²erg s⁻¹ cm⁻².

The derived slope in the hard X-ray band is practically identical to the MOJAVE-1 data set with α_{Moj1.5, BAT} = 1.13 ± 0.03. It can thus be concluded that the slightly different selection criteria of the earlier MOJAVE-1 sample are not a significant influence on the relatively flat slope of the X-ray log N–log S distribution. Both samples share the behavior of flux distribution in space. It is therefore highly suggested that the samples depict the same blazar population.

4.5.2. Influence of flaring sources

Since blazars are known to exhibit strong flaring behavior on a variety of time scales and because the BAT fluxes are derived from an integration over 105 months, it is reasonable to assume that at least a part of the blazar sample is characterized by fluxes that are distinctively higher compared to their quiet state. These increased fluxes may distort the log N–log S distribution, moving data points toward higher fluxes, and, consequently, change the distribution’s slope. A long-term variability study of AGN at hard X-rays, conducted by Soldi et al. (2014), revealed significant variation of source flux in the 14 keV–195 keV band of the first 66 months of light curve data collected by Swift/BAT. The authors calculated the variability estimator S_v, which expresses the intrinsic variability of a light curve of a source. The estimator is based on the maximum-likelihood estimate of the variability parameter σ_Q, which depends on the individual light curve flux values and measurement errors (see Soldi et al. 2014, Eq. (1)). This parameter is then normalized by the average flux of the specific source. An average value ⟨S_v⟩ in per cent for a sample of sources can then be calculated. The blazar sub-sample in their study showed the highest variability (⟨S_v⟩ = 33%±2%) as well as an exceptionally high variability of two gamma-ray bright blazars (⟨S_v⟩ = 90%). About 80% of the surveyed AGN are found to have variability on the time scale of months to years.

In order to estimate the influence of randomly increased flux within our sample we assume that 80% of a sample are variable with the maximum flux increase factor equal to the variability estimator S_v, for instance, a factor of 1.33 for S_v = 33%. From our sample of 70 MOJAVE-1 blazars in the log N–log S distribution we calculate S_v for the available monthly binned Swift/BAT light curves in the 105-month time span, and additionally the 70-month time span⁷. The results are presented in Table 3.

Based on the number of available light curves, it can be stated that the variability amplitude estimator S_v for Fermi/LAT-detected sources (on the basis of the 3FGL catalog) is approximately 33%–55%. The results for Fermi/LAT-non-detected sources differ significantly from approximately 30% in the 70-month data set to almost 80% in the 105-month data set. This, however, is likely due to the very low number of only three light curves per data set. Furthermore, the longer light curves also tend to have larger variability estimators. As a conservative estimate we assume that the variability amplitude does not exceed 100% for all light curves, equal to a maximum flux increase of a factor of two.

We start from a perfectly Euclidean flux distribution of 500 sources. The brightest source is set to 10 arbitrary flux units, the next source at a flux that is determined by following a slope in the log N–log S diagram of −1.5. In this simulation we produce two separate effects. First, the flattening of the log N–log S distribution at the lower end, which is caused by graphing the median radio flux whereas the selection of the sample was purely based on a possibly one time high state of a source, pushing the flux reading over the registration threshold. Second, a random flaring of 80% of all sources in the simulated sample, as outlined for the X-ray band.

A simulated flaring is induced in 400 randomly picked sources with a maximum random factor for the flux increase in the range of [1–2] and a uniform probability distribution. We simulate a large number of instances of the first and the second flaring behaviors and fit a power law to each of the resulting distributions. All data points with fluxes larger than 0.5 of the arbitrary flux units are included in the fit (see Fig. 11). A highly significant majority of 99.99% of all fitted instances have power-law slopes larger than 1.28, which is thus not compatible with the fit of the BAT data (1.13 ± 0.04), within 4σ of the slope’s uncertainty range. This result makes a simple flaring scenario of the X-ray emission unlikely to be solely responsible for the small slope of the X-ray log N–log S distribution.

Fig. 11. Simulated sample of sources with uniform luminosity and ideal Euclidean log N–log S flux distribution (red) and one instance of the same sample with random flaring (black, see text). The shaded areas indicate the error ranges derived by the Monte-Carlo approach, with 99.99% of all fitted power laws lying inside the outer shaded area.

4.5.3. High X-ray and low radio flux sources

For the log N–log S distribution of the BAT data to be steeper (more Euclidean), either faint sources are missing or bright sources are too numerous. Flux values that are only upper limits, and are consequently not in the log N–log S diagram, cannot be responsible for the missing faint sources since the vast majority of all upper limits are located below 4 × 10⁻¹² erg s⁻¹ cm⁻², a region that is still very well described by the power-law fit.

The potential over-abundance of X-ray-bright sources compared to the Euclidean radio measurements is examined in the following. As indicated by Fig. 12, no noticeable correlation of radio against X-ray flux is present in the whole MOJAVE-1 source sample. Also, a Spearman’s rank correlation coefficient of ρ = 0.31 does not indicate any significant link between the flux values in both bands. The figure shows a scatter plot with the flux data of both bands covering little more than two decades. At the same time, the distribution along the flux axes in both bands is dissimilar. Compared to the radio data set, the BAT fluxes (excluding upper limits) are concentrated at low values, below 10⁻¹¹ erg s⁻¹ cm⁻².

Fig. 12. MOJAVE-1 sample: BAT 20 keV–100 keV flux plotted against 15 GHz flux density. The symbols with arrows to the left represent upper limit values. The shared area indicates a selected number of X-ray-bright while radio-faint sources (see text).

The gray area in Fig. 12 marks X-ray fluxes greater than 10⁻¹¹ erg s⁻¹ cm⁻² and radio flux densities smaller than 5 Jy, comprising 10 sources (9 FSRQs and BL Lac itself). If this set of X-ray bright and radio-faint blazars is removed, the resulting slope of the log N–log S distribution is 1.54 ± 0.06, well compatible with the Euclidean case.

The scenario of a number of (bright) X-ray sources that govern the behavior of the blazar sample in the log N–log S plot is equal to a population that shows intrinsic evolution in luminosity. The following section analyzes the sample in terms of evolution of luminosity per unit co-moving volume.

4.5.4. LF data fit

We can exclude a dominant influence on the difference of the log N–log S shapes through variability. In the following, we present the results of the test for intrinsic emission evolution. In fitting the hard X-ray (XLF) and radio (RLF) luminosity function, we have two goals: first, study the evolution of the respective emission in the MOJAVE-1 blazar sample, and second, answer the question whether the difference in slope in the log N–log S distributions can be attributed to a different evolution in flux emission in both bands.

We analyze the data sets for the hard X-ray and 15 GHz bands for the MOJAVE-1 blazars. AIC values are calculated for all described models as well as the probability p_j that model j also maximizes the AIC value compared to the best fit model. The fitted models are as described in Sect. 3.3 for single power laws: PDE, PLE, PDEg, PLEg, and the PLE case with no evolution (k = g = 0). The double power-law models are generally less likely to produce an optimal fit, but are often not significantly less likely than the best fit (p_i ≈ 0.05−0.70). However, we choose to neglect these models because of the additional number of free parameters which do not introduce significant improvements in the fits. The fitting parameters and all of the following results are also listed in Table 4.

The best fit for the X-ray data set (69 sources⁸) is achieved with the single power-law PDEg and PLEg models. Double power-law models resulted of PLE and PLEg are considerably less probable (p_PLE = 0.247, p_PLEg = 0.320). The model corresponding to no evolution can be excluded (p_{no evol} = 0.0269).

The same sub-sample at radio frequencies gives a somewhat different picture. A simple evolution scenario (PLE, PDE) is preferred, with the second most likely cases of PLEg and PDEg considerably close (p_PLEg = p_PDEg = 0.737). The no evolution scenario can also not be fully excluded (p_no evol = 0.12).

The analysis of the larger sample of 123 MOJAVE-1 blazars⁹ in the radio band shows that the evolutionary models of PLEg and PDEg are the preferred case against no evolution (p_{no evol} = 0.003). Also, the largest available sample from Lister et al. (2015), counting 170 blazars with redshift information, shows that PLEg is the preferred model against the no evolution case p_{no evol} = 0.0049).

Figure 13 shows the binned XLF and best fit model of the XLF for a number of redshifts. Both the binned XLF and the model fit show the general falling trend typical for luminosity functions of this type. Luminosity and density evolution fit the data equally well. In any case, it can be stated that evolution of the emission is favored by the analytic models compared to the static one.

Fig. 13. MOJAVE-1 XLF model PLEg, one of the two best fit models. Binned LF and analytic form are fitted to the 70 source sub-sample for which BAT fluxes are non upper limits.

In Fig. 14, we graph the evolutionary factor e against redshift for all analyzed samples and for most of the luminosity function models in Table 4. For the 69 source X-ray data set both best-fit models show a trend of positive evolution, that is, an increase in X-ray luminosity or density up to z ≈ 1.5 and a negative trend after this point towards higher redshifts where sources tend to be less dense or luminous. Previous studies (e.g., Ueda et al. 2003; Miyaji et al. 2015; Ranalli et al. 2016) regarding the XLF of AGN samples also have applied a model of luminosity-dependent density evolution. However, because of the relatively low samples size, and additional number of parameters in such models, we choose to only use the simplest available models for the XLF.

Fig. 14. Evolution parameter e as a function of redshift for some of the fitted luminosity function models. The identification in the upper left corner corresponds to the samples in Table 4. The models that describe the data best are marked with an * in the legend.

Generally, all analyzed samples prefer an evolutionary scenario against no evolution at all. The preferred RLF models of the data sets of the 123 and 170 source samples show a peak in the at z ≈ 0.6. The likeliest evolution models for both radio samples and the 69 source X-ray data set show that the luminosity output per volume increased up to a point in time for the given population and decreased again. However, the same 69 sample of the hard X-ray bright sources in the radio band shows no such peak. The evolutionary parameter e only decreases with redshift. In any case, it is suggested that in terms of emission, that is, intrinsic luminosity, the MOJAVE-1 blazar sample behaves differently in the radio and hard X-ray regime.

However, with the available data, we cannot fully exclude the contribution of a selection effect that is given by the reduction of the full radio sample to the 69 sources that are analyzed in the X-ray band. In any way, the conclusion of a wavelength-dependent evolution through the analysis of the luminosity functions and the fact that the X-ray sample includes over-proportionally many bright sources, leads us to believe that this explanation is consistent and viable.

Assuming that a truly different evolutionary path for the radio and hard X-ray emission is present, that is, a steady decrease in the radio band with redshift or a small maximum at low redshifts but a large increase in luminosity of the population per volume for X-rays up to z ≈ 1.5, would imply a relatively low proportion of low flux sources in the X-ray data. This relation can thus be seen to be at least partially responsible for the lack of low flux sources and therefore flat slope of the log N–log S distribution.

We also test for evolution of the 69 X-ray bright MOJAVE-1 blazars using a standard V/V_max test (Schmidt 1968). The result of the test value and its mean absolute error ⟨V/V_max⟩ = 0.451 ± 0.184 suggest, although not significantly, a qualitative agreement with our previous results of (negative) evolution of the hard X-ray emission.

4.6. CXB contribution

We calculate the contribution to the hard X-ray background from blazars of the MOJAVE-1 type and the BAT catalog beamed AGN sample for comparison. The MOJAVE-1 blazar sub-sample yields the following results: Using the fit of the hard X-ray log N–log S data set (see Fig. 10), we extrapolate the number of sources brighter than 2.2 × 10⁻¹² erg s⁻¹ cm⁻² (flux of the faintest source) for the entire sky and obtain:

$$\begin{array}{c}\hfill N_{\mathrm{all}\mathrm{sky}}^{\mathrm{MOJAVE}}(>2.2\times {10}^{-12})=N(>2.2\times {10}^{-12})·{\mathrm{\Omega }}_{\mathrm{sky}}\approx 139\pm 9.\end{array}$$(20)

This number exceeds the actual total number of MOJAVE-1 sources for which a hard X-ray flux could be calculated (70 sources; see Sect. 4.5.1) but is consistent with the larger total number of sources contributing to the BAT signal (Sect. 4.1). The integrated flux of radio-selected blazars of the MOJAVE-1 type above 2.2 × 10⁻¹² erg s⁻¹ cm⁻² over the whole sky can thus be calculated from Eq. (10) using the parameters in Table 2:

$$\begin{array}{cc}& F_{\mathrm{contrib}}^{{\mathrm{deg}}^2,\mathrm{MOJAVE}}={\displaystyle {\int }_{2.2\times {10}^{-12}}^{\infty }\alpha ·A·F^{-\alpha }\mathrm{d}F}\hfill \\ \hfill & =(6.{47}_{-1.51}^{+2.83})\times {10}^{-14}\mathrm{e}rg{s}^{-1}c{m}^{-2}de{g}^{-2}\hfill \end{array}$$(21)

which is equivalent to 0.21% of $F_{\text{CXB},\hspace{0.17em}\text{20-100}\hspace{0.17em}\text{keV}}^{{\text{deg}}^2}$, with the uncertainty range of 0.16% to 0.30%. The contribution of MOJAVE-1-type blazars to the CXB increases only slowly with lower integration thresholds.

In order to reach a total blazar contribution to the CXB of 20%, as determined by Ajello et al. (2009) based on X-ray selected blazars in the 15 keV–55 keV band, the lower flux integration limit needs to be F = 1.15 × 10⁻²⁷ erg s⁻¹ cm⁻² for the measured slope of −1.13 of the MOJAVE-1 blazars. A number of 10¹⁴ sources per deg² would be needed to contribute to the signal, which is orders of magnitude larger than the estimated total number of 2 × 10¹² galaxies up to a redshift of z = 8 (Conselice et al. 2016). In contrast, a log N–log S power law, normalized on the brightest source and with a slope of −1.5 would correspond to a more realistic minimum integration flux of F = 5 × 10⁻¹⁵ erg s⁻¹ cm⁻², or about 421 sources per deg².

In the study by Ajello et al. (2009), the blazar log N–log S distribution was likewise characterized by very steep slopes (α= 1.7–2.0) and substantial contribution percentages. In the same study a contribution of Seyfert galaxies in the same energy band (15 keV–55 keV) was determined to be approximately 22% to 55%, depending on the specific evolutionary model, already putting another restriction on the maximum value for blazars. The noticeably different results of blazar CXB contribution of our study and that of Ajello et al. (2009) do not contradict each other because of the selection criteria for each blazar sample, that is, radio- vs. X-ray-bright sources. Finally, we compare the results of the MOJAVE-1 blazar sample to the log N–log S distribution of the BAT 105-month source catalog beamed AGN sub-sample using the catalog flux values in the 14 keV–195 keV band. The analysis shows a contribution to the CXB of 0.31% (0.29%–0.35%) when integrating starting at the faintest source flux in the sample (6.3 × 10⁻¹² erg s⁻¹ cm⁻²). Removing all MOJAVE-1 sources from that sub-sample reveals a contribution of 0.20% (0.18%–0.23%). Both the MOJAVE-1 and BAT beamed AGN samples have a very similar flux contribution to the CXB. For reaching a total blazar contribution of 20% a lower integration limit of F = 2.5 × 10⁻¹⁹ erg s⁻¹ cm⁻² is required, corresponding to a source count of 7.3 × 10⁶ per deg². Both the radio-selected MOJAVE-1 sample and the BAT catalog beamed AGN show two distinct groups of sources that, through their different log N–log S slopes, contribute in a significantly different way to the CXB in terms of necessary numbers.

5. Discussion

5.1. Individual sources with peculiar hard X-ray properties

The distribution of photon indices for sources shown in Fig. 4 is very similar to the distribution for the soft X-ray regime (2 keV–10 keV) as obtained by Swift/XRT (Chang 2010). A notable difference between both distributions is a softer spectrum for the radio galaxy sub-set at hard X-rays whereas FSRQs and BL Lacs are in general broader distributed with no significant offset to the soft X-ray data. A peculiar case is set by the X-ray-bright radio galaxy 0316+413 (3C 84) at S/N = 58.2 <  SPSDOUBLEDOLLAR >  σ with Γ = 3.17 ± 0.12. The poor fit statistic of the power-law fit (PG = 14.8), however, indicates that a more complex fit model and measurements of higher spectral resolution is needed. Similar fit results are presented in the 105-month BAT source catalog with Γ = 3.82 ± 0.09 and the reduced fit statistic χ² = 3.8 (Oh et al. 2018). Since the X-ray emission of this source in particular is not dominated by synchrotron or inverse-Compton (IC) processes, but by thermal emission (Churazov et al. 2003) we do not expect a typical blazar spectrum.

Another source with relatively soft photon index is the BL Lac 0754+100 with Γ = 3.04 ± 1.56 at S/N = 3.9 <  SPSDOUBLEDOLLAR >  σ. The fitted SED from Chang (2010) suggests a HE emission bump smaller and flatter than the synchrotron bump in νF_ν scaling with the BAT energy band at the rising part of the HE bump. Considering the large error a BAT photon index Γ <  2 seems reasonable and more realistic.

The BL Lac 0716+714 has the hardest spectrum of all sources in this class with Γ = 0.83 ± 0.46 while being still reasonably bright at 5.4σ. Additionally, the very strong variability of this source (e.g., Wagner et al. 1996) makes meaningful spectral measurements difficult. A number of studies in the hard X-ray band in the past have not been able to derive a spectral index above 10 keV due to variability issues and low significance. The second INTEGRAL AGN catalog (Beckmann et al. 2009) lists 0716+714 as a source of strong variability (>0.5 mag) in the optical photometric V-band. The authors analyzed the data and derived a flux value under the assumption of a photon index of Γ = 2 because of low significance in the observed 18 keV–60 keV band of less than 5σ. In another study Pian et al. (2005) derived the flux of the source in the range of 30 keV–60 keV using high-state-triggered observations. This was also done by freezing the photon index to the known value of the Crab pulsar Γ = 2.1, because of low significance. The very hard photon index of our measurement can be attributed to the position of the BAT band at a very steep rising part of the HE bump. Recent observations of the source with Swift/XRT (0.3 keV–10 keV) and NuSTAR (3 keV–50 keV) in a flaring state (Wierzcholska & Siejkowski 2016) have clearly shown the spectral minimum between both emission bumps in the covered X-ray range at E_br = (8.01 ± 0.56) keV. The spectrum has been modeled with a broken power law with the resulting photon indices of Γ₁ = 2.40 ± 0.01 and Γ₂ = 1.61 ± 0.05. The found break energy E_br has been the highest value ever measured. A previous study by Wierzcholska & Siejkowski (2015) found break energy values between 2 keV and 5 keV during another flare state. Also, a direct correlation of elevated flux levels and break energy was found. The position of the very steep rising part of the HE bump can be placed within the BAT band (20 keV–100 keV) with high confidence. The detection of this source by Fermi/LAT (Lister et al. 2009b) supports this assumption, since the HE bump is comparable to the majority of the sample, shifted to higher energies, allowing Fermi/LAT to detect more high energy photons.

A relatively broad distribution of spectral shapes for low-peaked sources was also reported in a recent study by Marchesini et al. (2019). The authors analyzed a large sample of high- and low-peaked BL Lac type objects in the soft X-ray band and found an especially broad distribution of hardness ratios for the log-peaked sub-set, which has been attributed to the intersecting region of synchrotron and HE bump in the SED.

5.2. Fermi detections of gamma-faint blazars

Despite the substantial integration time of the LAT instrument during the four years (3FGL catalog), some of the most luminous blazars at radio wavelengths are not detected by Fermi/LAT (Sect. 2). Lister et al. (2015), using the extended MOJAVE-1.5 sample, found that the synchrotron bump of the Fermi/LAT-non-detected radio-loud AGN tend to be characterized by a lower peak frequency than the Fermi/LAT-detected part of the sample. This was partially explained by the SED position and consequently lower flux values above 0.1 GeV besides low relativistic boosting of the jet and its velocity. Also, low radio variability levels were associated with a low probability of detection in the gamma-ray range.

In Sect. 4.2.2, we have presented the distribution of the measured BAT photon indices of all fitted MOJAVE-1 blazars, Fermi/LAT-detected and non-detected. The analysis, using the extended 2-sample KS testing method, showed that both Fermi/LAT-detected and non-detected sub-samples exhibit significantly different distributions, suggesting that detected sources (harder) are associated with higher peaked and non-detected sources (softer) with lower peaked SEDs. The shown distribution of HE peak frequencies against BAT photon index in Fig. 6 supports this trend. A very similar relation has been found by Abdo et al. (2010). The authors report a strong correlation of gamma-ray photon index and the synchrotron as well as HE bump peak frequency. The higher (softer) the index, the lower the peak frequency. Furthermore, the gamma-ray photon index was found to be strongly correlated with the X-ray photon index.

Previous to the publication of the forth Fermi/LAT source catalog (The Fermi-LAT Collaboration 2020), Lister et al. (2015) identified three blazars in the extended 1.5 Jy MOJAVE sample which were thought to be strong candidates for a future detection by Fermi/LAT on the basis of their observed properties, that is, high radio modulation index and high apparent jet speed. None of the three mentioned blazars (III Zw 2, PKS 0119+11, and 4C +69.21) are listed in the 4FGL catalog, which includes approximately double the number of blazars compared to the 3FGL/3LAC catalog.

We propose a different approach for predicting which blazars in the MOJAVE-1 sample will most likely be detected in future Fermi catalogs. We compare the HE peak frequencies and BAT photon indices of non-detected with detected sample sources. The two previously non-detected FSRQs (3LAC) 1458+718 and 2145+067 have HE peak frequencies comparable with the majority of detected sources, as seen in Fig. 6. However, when compared with the rest of all MOJAVE-1 blazars, both sources have relatively low apparent maximum jet speeds of 6.6c and 2.8c, respectively. As indicated in the plot both sources have since been detected, reinforcing the hypothesis that HE peak position is likely the dominant factor for a significant gamma-ray signal.

Based on our findings, we can predict likely gamma-ray detections of not yet detected blazars in the MOJAVE-1 sample. The mean of the maximum jet speeds of all Fermi/LAT-detected MOJAVE-1 blazars is 12.8c and a maximum value of 41.8c (Lister et al. 2013). A source with a hard BAT spectrum (Γ = 1.50 ± 0.44) and reasonably high apparent jet speed of 9.9c is the FSRQ 2005+403. Its HE peak frequency is on the order of a magnitude lower than the bulk of all detected MOJAVE-1 blazars, although higher than the majority of all non-detected ones. Also, the source 0742+103 can be considered a detection candidate because of its relatively high peak frequency. This case is of particular interest, since the previously proposed criteria of high jet speed and variability index (v = 2.8c and m = 0.016) are far from ideal considering the mean values of these quantities in the mentioned study (v_mean = 10.1c and m_mean = 0.2). A future gamma-ray detection would certainly emphasize the importance and impact of soft and hard X-ray coverage in order to estimate the HE peak position and, thus, the gamma-ray detection probability for a blazar.

5.3. log N–log S distribution and LF

The log N–log S distribution of the MOJAVE-1 blazars is characterized by a flat and non-Euclidean slope in the 20 keV–100 keV band (α_Moj, BAT = 1.13 ± 0.04). Most X-ray AGN surveys, which are usually dominated by Seyferts (e.g., surveys with INTEGRAL; Beckmann et al. 2006 and Swift/BAT; Ajello et al. 2012), show a slope that is compatible with a Euclidean distribution in space. However, the beamed AGN sub-sample of the hard X-ray BAT source catalog (Oh et al. 2018) is also described by a relatively flat fit function to the log N–log S distribution (α_cat, BAT = 1.24 ± 0.02), reinforcing this peculiar trend of blazars in this energy range. It could be argued that a substantial number of the BAT catalog’s unidentified 243 sources are indeed blazars, which would likely influence the distribution’s slope. The fraction of unidentified sources in the MOJAVE-1 sample, on the other hand, is with the two sources 0648–165 and 1213–172 relatively small.

The radio data set incorporates the full set of all 125 MOJAVE-1 blazars, whereas the BAT data that we used comprise the 70 brightest blazars in the 20 keV–100 keV energy band. One likely effect responsible for the different slopes is the measurement of two different frequency bands. The sample is statistically complete in the radio band by definition. The X-ray fluxes, on the other hand, do not correlate strongly with the radio emission. The plot of BAT flux vs. VLBI MOJAVE flux at 15 GHz (Fig. 12) shows no significant positive correlation of both data sets for any of the sub types. However, a majority of all the sample is concentrated at low X-ray flux (Sect. 4.3) whereas these sources scatter over a wide range of radio fluxes (approximately 0.6 Jy–5 Jy). The position of data points in the 15 GHz log N–log S diagram are consequently shifted towards low X-ray flux, and in this case even below the set flux limit of the X-ray log N–log S plot, lowering the number of sources at the bright end, and thus reducing the distribution’s slope. This on its own could lead to the conclusion that the difference between α_Moj, BAT and α_{Moj, 15 GHz} is due to source selection bias of a sample in a different wavelength regime. However, the very similar slope of the log N–log S distribution for the MOJAVE-1 sample and the beamed AGN sub-sample of BAT 105-month source catalog make a pure selection bias unlikely. The influence of flaring sources can be neglected, as well. A likely conclusion is a different evolutionary path for X-ray and radio emission in the source sample. Whereas the preferred XLF model indicates a maximum of the evolution parameter e at around z = 1.5 the RLFs of the blazars of the full MOJAVE-1 and the extended MOJAVE-1.5 samples have maxima at approximately z <  1.

In a recent study, Ighina et al. (2019) have found a significant difference of the evolutionary trend of soft X-ray to radio luminosity ratio in a high-redshift blazar sample (z >  4) compared to a sample at low redshifts. The ratio was found to be approximately 2–3 times larger for the high-redshift sample. This trend, although of the soft X-ray band (0.2 keV–10 keV), generally supports our results, which indicate an earlier X-ray emission maximum compared to the radio data. The authors tentatively attribute the decreasing fraction of X-ray to radio emission with the interaction of extended jet regions with the CMB by the IC process, which is strongly dependent on the redshift. However, Ighina et al. (2019) also make the argument that, in order to reconcile previous results from radio- and X-ray-selected source samples (Ajello et al. 2009; Caccianiga et al. 2019), the ratio of extended to compact emission has to be non-uniform, likely depending on the specific source and the source selection.

In general, the log N–log S distribution stays consistent to very low fluxes in the BAT band. Thus, the behavior of the fainter sources in the plotted sample is not significantly different to the rest. This consistency reinforces our approach of subtracting the fitted noise component of the S/N distribution from the measured S/N distribution of the blazar sample (Sect. 4.1). Also, the analysis of the extended MOJAVE-1.5 sample shows the same slope of the distribution, excluding any influence that might have been introduced by non-uniform source selection criteria in the original MOJAVE-1 sample.

The contribution of AGN, blazars, and other source types to the CXB has primarily been examined using source samples that are also X-ray or gamma-ray bright (see, e.g., Mateos et al. 2008; Ajello et al. 2009; Bottacini et al. 2012). The results for the blazar contribution to the hard CXB vary greatly in previous works, largely due to different source samples and analytic approaches. In this work we use a sample of radio-bright low-peaked blazars and obtain up to 0.30% of the contribution in the 20 keV–100 keV band above 2.2 × 10⁻¹² erg s⁻¹ cm⁻², which again is compatible with the beamed AGN sub-sample of the 105-month BAT survey catalog.

While the contribution to the approximately 20% of the CXB from blazars (Ajello et al. 2009) can be reached with a low enough X-ray flux integration limit of F = 5 × 10⁻¹⁵ erg s⁻¹ cm⁻² for a Euclidean slope, the measured slope of −1.13 from the MOJAVE-1 sources is too flat for a physically reasonable result (Sect. 4.6). Our results thus suggest that X-ray and radio-selected blazars behave like two two different populations of which the radio-selected population makes only a negligible contribution to the CXB. The parameters of the fits of the luminosity functions to the MOJAVE-1 BAT data show a key difference to comparable blazar studies of in soft and hard X-ray bands. Studies by, for example, Ajello et al. (2009), Sazonov et al. (2007), and Hasinger et al. (2005) revealed very similar parameter values for the evolutionary term of the luminosity functions of larger AGN and blazar samples. However, the slopes of the fitted power law are usually relatively steep, that is, around γ = 2−3. If a double power law was applied the index of the lower luminosity part was flat at around γ = 0.5. In contrast, the most likely models of the fitted MOJAVE-1 BAT data show indices of approximately γ = 1. The analyzed radio data sets reveal even flatter indices of γ = 0.5−0.65. The fitted blazar luminosity functions from Ajello et al. (2009) for the energy range of 15 keV–55 keV show very similar source densities at low luminosities around 10⁴⁴ erg s⁻¹. For higher luminosities of 10⁴⁸ erg s⁻¹ the source density is approximately 10³ times lower compared to our studies. The relative scarcity of low luminosity sources in the MOJAVE-1 sample could also be attributed to the 52 blazars that have been omitted in the analysis because of low signal strength (upper limits). The radio luminosity function of the complete MOJAVE-1 blazar sample of 122 sources shows a similar behavior without the possible influence of missing contributing sources, however. It is therefore highly suggested that the X-ray-selected blazar samples of previous studies are of a different blazar population and evolutionary behavior than the radio-selected MOJAVE-1 sample.

Compared to the MOJAVE-1 blazar sample, the 38 hard X-ray selected sources in the study by Ajello et al. (2009) show a large amount of very low (18 with z <  0.5) and high redshift sources (10 with z >  2). The HE emission bump of high-redshift sources is shifted down to the BAT detection band. Also, the included BL Lac type sources are almost all HBLs. This creates a fundamentally different selection of sources and spectral characteristics in the hard X-ray regime relative to the radio-selected sample with a much more even redshift distribution in this study. Again, this shows that the different results of positive blazar evolution (Ajello et al. 2009) against the negative evolution (MOJAVE-1, hard X-rays) do not directly contradict each other.

6. Summary and conclusions

In the following, we list the major results of our study of the hard X-ray properties of the radio-selected and statistically complete MOJAVE-1 AGN sample, consisting mainly of low-peaked blazars. We present the spectral characteristics of this under-represented and rare AGN type at hard X-rays in the 20 keV–100 keV band, based on the Swift/BAT 105-month survey maps.

• The X-ray flux, luminosities, and photon indices for the MOJAVE-1 AGN sample are now available for the first time in this spectral band. Out of the 135 sources in the sample, 121 are characterized by a hard X-ray signal that is not compatible with background noise, although most sources have low S/N values.

• The derived distributions of hard X-ray flux, luminosity, and photon index clearly distinguish FSRQs, radio galaxies, and BL Lacs. Furthermore, the distribution of photon indices is significantly different for Fermi/LAT-detected and non-detected sources, indicating that the position of the HE peak in the SED is one likely contributor to the question why certain radio-bright AGN are not also gamma-bright.

• The source count distribution/log N–log S for the 70 X-ray-brightest blazars in the MOJAVE-1 sample shows a relatively flat slope of α = 1.13 ± 0.04, which is clearly not Euclidean. One possible main contributor responsible for the difference in slope, and consequently the scarcity of hard X-ray detected blazars, is the different evolution of hard X-ray emission, peaking at around z = 1.5.

• We derive the contribution of the hard X-ray flux for the MOJAVE-1 and beamed AGN BAT catalog samples to the CXB and find that approximately 0.2%–0.3% of the CXB can be resolved into the sample sources for fluxes greater that their respective survey limits. Compared to earlier studies with X-ray-selected blazar samples this contribution is significantly lower, which suggests a different blazar population due to the selection criteria of the MOJAVE-1 survey, that is, low-peaked and radio-loud compact and beamed sources.

The obtained properties for the various AGN types will become important for further research in areas such as, for example, the blazar sequence and broadband SED modeling, especially near the spectral minimum between synchrotron and HE emission peak. Large and deep future X-ray surveys, using, for example, the recently launched eROSITA X-ray telescope, will further complete the understanding the source types that are responsible for the creation of the CXB. Since the soft and hard X-ray regimes are instrumental for the estimation of future gamma-ray detections of beamed AGN in terms of SED position, complementary studies of both bands are essential for finding interesting new targets for gamma-ray studies using GeV and TeV observatories, such as the currently constructed Cherenkov Telescope Array (CTA).

Acknowledgments

We want to thank Moritz Böck and Takamitsu Miyaji for useful discussion of the analysis methods. We also thank the anonymous referee for many useful comments and suggestions that helped improving the paper. This work was partially funded by Deutsche Forschungsgemeinschaft grant 50 OR 1709. This research has made use of data from the MOJAVE database that is maintained by the MOJAVE team (Lister et al. 2018) and the ISIS functions (ISISscripts), provided by the ECAP/Remeis observatory and MIT (http://www.sternwarte.uni-erlangen.de/isis/). The Swift/BAT survey maps have been provided by the Swift team. Furthermore, this research has made use of the NASA/IPAC Extragalactic Database NED (Jet Propulsion Laboratory, California Institute of Technology/NASA) and of the SIMBAD database (CDS, Strasbourg, France). The work of M. Kreter is supported by the South African Research Chairs Initiative (grant no. 64789) of the Department of Science and Innovation and the National Research Foundation of South Africa. Any opinion, finding and conclusion or recommendation expressed in this material is that of the authors and the NRF does not accept any liability in this regard.

References

Appendix A: Hard X-ray data for the MOJAVE-1 sample

All Tables

All Figures

	Fig. 1. Swift/BAT 105 month survey sky coverage as the function of minimal flux corresponding to a detection at a level of 1, 3, and 5σ in the band 20 keV–100 keV, from the left to the right. The dashed line indicates the full sky.
In the text

	Fig. 2. BAT S/N distribution of MOJAVE-1 sample for different source classes. The brightest source 3C 273 at 192σ has been omitted for better readability.
In the text

Fig. 3. Top panel: distribution of Swift/BAT S/N values of 1000 random blank sky positions. The red line indicates a Gaussian fit with center of 0σ and width 1σ. The red shaded region indicates 3σ uncertainty ranges for the norm of the fitted function. Middle panel: distribution of Swift/BAT S/N values of the MOJAVE-1 sample. The red line indicates a Gaussian fit with center of 0σ and width 1σ to the distribution with values smaller than 0σ. The red shaded region indicates 3σ error ranges. The 4.8σ threshold of the 105-month BAT survey is shown by the dashed line. Bottom panel: resulting histogram after subtracting the MOJAVE-1 S/N distribution by the distribution described by the Gaussian fit. Arrows indicate the cutoff at 15σ (for better readability).

In the text

	Fig. 4. Distribution of photon index Γ for all fitted MOJAVE-1 sources. Hatched boxes indicate sources with S/N values below 4σ (fainter 50% of the sample).
In the text

	Fig. 5. Distribution of BAT photon index Γ for all fitted MOJAVE-1 blazars, detected and not detected by Fermi/LAT, respectively. Bottom panel: added Gaussian distributions for each photon index. The sharp peak at 1.75 stems from the brightest blazar in the sample, 3C 273.
In the text

	Fig. 6. SED HE peak frequency of MOJAVE-1 blazars (rest frame) against BAT photon index Γ.
In the text

	Fig. 7. Distribution of hard X-ray flux of the MOJAVE-1 sample.
In the text

	Fig. 8. Relation of hard X-ray flux to S/N. The dashed line represents the S/N value of 4.8σ. The source 1504–166 at 0.01σ has been omitted for better readability.
In the text

	Fig. 9. K-corrected X-ray luminosity distribution for the energy range of (20–100) keV of the MOJAVE-1 sample.
In the text

Fig. 10. Left: cumulative log N–log S distribution of the MOJAVE-1 sample excluding radio galaxies, omitting all upper limits (flux for 20–100 keV, blue squares) and log N–log S for the beamed AGN sub-sample of the BAT 105-month source catalog (flux for 14 keV–195 keV, red circles) and the same sub-sample excluding all MOJAVE-1 sources (14 keV–195 keV, purple x). Right: log N–log S for 15 GHz flux density VLBI measurements from Lister et al. (2015) of all blazars in the MOJAVE-1 sample (black circles) and MOJAVE-1 blazars from left graph (blue squares, norm with a factor of 0.1 for better readability). The black lines represent power-law fits, the gray area shows the corresponding error range. The radio data have been fit for all fluxes higher than 2.5 Jy. The dashed lines indicate a slope of −1.5.

In the text

	Fig. 11. Simulated sample of sources with uniform luminosity and ideal Euclidean log N–log S flux distribution (red) and one instance of the same sample with random flaring (black, see text). The shaded areas indicate the error ranges derived by the Monte-Carlo approach, with 99.99% of all fitted power laws lying inside the outer shaded area.
In the text

	Fig. 12. MOJAVE-1 sample: BAT 20 keV–100 keV flux plotted against 15 GHz flux density. The symbols with arrows to the left represent upper limit values. The shared area indicates a selected number of X-ray-bright while radio-faint sources (see text).
In the text

	Fig. 13. MOJAVE-1 XLF model PLEg, one of the two best fit models. Binned LF and analytic form are fitted to the 70 source sub-sample for which BAT fluxes are non upper limits.
In the text

	Fig. 14. Evolution parameter e as a function of redshift for some of the fitted luminosity function models. The identification in the upper left corner corresponds to the samples in Table 4. The models that describe the data best are marked with an * in the legend.
In the text