Bimodal radio variability in OVRO-40m-monitored blazars

TL;DR

This paper introduces a new maximum likelihood method to analyze bimodal radio variability in blazars, distinguishing quiescent and flaring states to better understand their behavior and link to gamma-ray emission.

Contribution

The novel approach models blazar flux-density as 'off' and 'on' states, enabling detailed analysis of variability characteristics and their relation to gamma-ray detection.

Findings

BL Lacs are more variable with larger outbursts than FSRQs.

Gamma-ray detected blazars show different variability properties from non-detected ones.

Flaring ratios, duty cycles, and modulation indices distributions are fitted for population studies.

Abstract

Blazars are known to show periods of quiescence followed by outbursts visible throughout the electromagnetic spectrum. We present a novel maximum likelihood approach to capture this bimodal behavior by examining blazar radio variability in the flux-density domain. We separate quiescent and flaring components of a source's light curve by modeling its flux-density distribution as a series of "off" and "on" states. Our modeling allows us to extract information regarding the flaring ratio, duty cycle, and the modulation index in the "off"-state, in the "on"-state, as well as throughout the monitoring period of each blazar. We apply our method to a flux-density-limited subsample from the Owens Valley Radio observatory's 15 GHz blazar monitoring program, and explore differences in the variability characteristics between BL Lacs and FSRQs as well as between $\gamma$-ray detected and non-detected sources. We find that: (1) BL Lacs are more variable and have relatively larger outbursts than the FSRQs, (2) unclassified blazar candidates in our sample show similar variability characteristics as the FSRQs, and (3) $\gamma$-ray detected differ from the $\gamma$-ray non-detected sources in all their variability properties, suggesting a link between the production of $\gamma$-rays and the mechanism responsible for the radio variability. Finally, we fit distributions for blazar flaring ratios, duty cycles, and on- and off- modulation indices that can be used in population studies of variability-dependent blazar properties.

Tables

Table 1. Table 1: Mean parameter values and probabilities for the Wilcoxon and K-S test that the S 𝑆 S and m 𝑚 m parameters estimated by the two models (single-Gaussian and bimodal) are consistent (drawn from the same distribution).

Parameter	Single	Bimodal	Wilcoxon	K-S
	(mean)	(mean)	(%)	(%)
$S$	1.11	1.11	90.0	99.5
$m$	0.18	0.20	0.28	1.26

Table 2. Table 2: Mean values for the duty cucle, flaring ratio and “off”- and “on”-state modulation indices for the entire sample

Parameter	$f_t$	$R$	$m_{\mathrm{off}}$	$m_{\mathrm{on}}$
Mean	0.47	1.56	0.10	0.10
Standard deviation	0.26	1.55	0.08	0.08

Table 3. Table 3: MLE exponential distributions for R 𝑅 R , m on subscript 𝑚 on m_{\rm on} , m off subscript 𝑚 off m_{\rm off} , for the entire flux-density limited sample.

Quantity	Mean	error on mean
$R$	1.485	$\pm$ 0.006
$m_{\mathrm{off}}$	0.080	$\pm$ 0.004
$m_{\mathrm{on}}$	0.089	$\pm$ 0.004

Table 4. Table 4: Mean parameters and Wilcoxon and K-S test probability values that the two population are consistent.

Parameter	BL Lacs	FSRQs	Wilcoxon	K-S
	(mean)	(mean)	(%)	(%)
$f_t$	0.46	0.47	92.3	86.1
$R$	1.71	1.55	0.3	0.8
$m_{\mathrm{off}}$	0.11	0.10	32.8	64.6
$m_{\mathrm{on}}$	0.14	0.10	${10}^{-3}$	0.2

Table 5. Table 5: Mean parameters and Wilcoxon and K-S test probability values that the two populations are consistent.

Parameter	Blazars	U-Rs	Wilcoxon	K-S
	(mean)	(mean)	(%)	(%)
$f_t$	0.47	0.50	32.9	13.3
$R$	1.57	1.44	0.08	0.01
$m_{\mathrm{off}}$	0.10	0.08	9.0	30.8
$m_{\mathrm{on}}$	0.11	0.10	46.7	33.1

Table 6. Table 6: Mean parameters and Wilcoxon and K-S test probability values that the two populations are consistent.

Parameter	BL Lacs	U-Rs	Wilcoxon	K-S
	(mean)	(mean)	(%)	(%)
$f_t$	0.46	0.50	34.4	33.9
$R$	1.71	1.44	0.007	0.009
$m_{\mathrm{off}}$	0.11	0.08	0.3	0.4
$m_{\mathrm{on}}$	0.14	0.10	5.9	27.0

Table 7. Table 7: Mean parameters and Wilcoxon and K-S test probability values that the two populations are consistent.

Parameter	Fermi	non-Fermi	Wilcoxon	K-S
	(mean)	(mean)	(%)	(%)
$f_t$	0.47	0.48	65.4	52.8
$R$	1.60	1.51	${10}^{-10}$	${10}^{-8}$
$m_{\mathrm{off}}$	0.11	0.09	0.009	0.26
$m_{\mathrm{on}}$	0.13	0.09	${10}^{-7}$	0.001

Table 8. Table 8: Variability parameters for the OVRO monitored sources. Columns: (1) OVRO name; (2) duty cycle ( f t subscript 𝑓 t f_{\rm t} ); (3) 1 σ 𝜎 \sigma error on f t subscript 𝑓 t f_{\rm t} from the information matrix; (4) upper 1 σ 𝜎 \sigma error on f t subscript 𝑓 t f_{\rm t} from the likelihood slice; (5) lower 1 σ 𝜎 \sigma error on f t subscript 𝑓 t f_{\rm t} from the likelihood slice; (6) off-state flux-density ( S off subscript 𝑆 off S_{\rm off} ); (7) 1 σ 𝜎 \sigma error on S off subscript 𝑆 off S_{\rm off} from the information matrix; (8) upper 1 σ 𝜎 \sigma error on S off subscript 𝑆 off S_{\rm off} from the likelihood slice; (9) lower 1 σ 𝜎 \sigma error on S off subscript 𝑆 off S_{\rm off} from the likelihood slice; (10) standard deviation of the off-state flux-density ( σ off subscript 𝜎 off \sigma_{\rm off} ); (11) 1 σ 𝜎 \sigma error on σ off subscript 𝜎 off \sigma_{\rm off} from the information matrix; (12) upper 1 σ 𝜎 \sigma error on σ off subscript 𝜎 off \sigma_{\rm off} from the likelihood slice; (13) lower 1 σ 𝜎 \sigma error on σ off subscript 𝜎 off \sigma_{\rm off} from the likelihood slice; (14) on-state flux-density ( S on subscript 𝑆 on S_{\rm on} ); (15) 1 σ 𝜎 \sigma error on S on subscript 𝑆 on S_{\rm on} from the information matrix; (16) upper 1 σ 𝜎 \sigma error on S on subscript 𝑆 on S_{\rm on} from the likelihood slice; (17) lower 1 σ 𝜎 \sigma error on S on subscript 𝑆 on S_{\rm on} from the likelihood slice; (18) standard deviation of the on-state flux-density ( σ on subscript 𝜎 on \sigma_{\rm on} ); (19) 1 σ 𝜎 \sigma error on σ on subscript 𝜎 on \sigma_{\rm on} from the information matrix; (20) upper 1 σ 𝜎 \sigma error on σ on subscript 𝜎 on \sigma_{\rm on} from the likelihood slice; (21) lower 1 σ 𝜎 \sigma error on σ on subscript 𝜎 on \sigma_{\rm on} from the likelihood slice. The parameter values for S off subscript 𝑆 off S_{\rm off} , σ off subscript 𝜎 off \sigma_{\rm off} , S on subscript 𝑆 on S_{\rm on} and σ on subscript 𝜎 on \sigma_{\rm on} are in Jansky. For uncertainty values in columns (3), (7), (11), (15), and (19) equal to zero, the information matrix failed to provide an error estimate. In this case, use upper and lower 1-sigma uncertainties from the likelihood slice. The full-version of the table is available online.

Name	$f_t$	${\sigma }_{f_t}$	upper-${\sigma }_{f_t}$	lower-${\sigma }_{f_t}$	$S_{\mathrm{off}}$	${\sigma }_{S_{\mathrm{off}}}$	upper-${\sigma }_{S_{\mathrm{off}}}$	lower-${\sigma }_{S_{\mathrm{off}}}$	${\sigma }_{\mathrm{off}}$	${\sigma }_{{\sigma }_{\mathrm{off}}}$	upper-${\sigma }_{{\sigma }_{\mathrm{off}}}$	lower-${\sigma }_{{\sigma }_{\mathrm{off}}}$	$S_{\mathrm{on}}$	${\sigma }_{S_{\mathrm{on}}}$	upper-${\sigma }_{S_{\mathrm{on}}}$	lower-${\sigma }_{S_{\mathrm{on}}}$	${\sigma }_{\mathrm{on}}$	${\sigma }_{{\sigma }_{\mathrm{on}}}$	upper-${\sigma }_{{\sigma }_{\mathrm{on}}}$	lower-${\sigma }_{{\sigma }_{\mathrm{on}}}$
(1)	(2)	(3)	(4)	(5)	(6)	(7)	(8)	(9)	(10)	(11)	(12)	(13)	(14)	(15)	(16)	(17)	(18)	(19)	(20)	(21)
J0001-1551	0.168	0.0373	+ 0.0369	- 0.0369	0.201	0.0018	+ 0.002	- 0.002	0.0	0.0048	+ 0.0037	- 0.0	0.222	0.0	+ 0.0044	- 0.0044	0.044	0.0	+ 0.0022	- 0.0022
J0001+1914	0.839	0.0208	+ 0.0252	- 0.0252	0.229	0.0017	+ 0.0023	- 0.0023	0.027	0.0014	+ 0.0016	- 0.0013	0.311	0.0	+ 0.0031	- 0.0031	0.009	0.0	+ 0.0023	- 0.0021
J0003+2129	0.978	0.012	+ 0.0196	- 0.0196	0.077	0.0006	+ 0.0008	- 0.0008	0.008	0.0005	+ 0.0006	- 0.0005	0.083	0.0	+ 0.0132	- 0.0124	0.033	0.0	+ 0.0115	- 0.0072
J0004-1148	0.205	0.0237	+ 0.0224	- 0.0216	0.366	0.0061	+ 0.0059	- 0.0059	0.04	0.0046	+ 0.0045	- 0.004	0.631	0.0	+ 0.0069	- 0.0069	0.113	0.0	+ 0.0055	- 0.0052
J0004+2019	0.59	0.0303	+ 0.0277	- 0.0283	0.289	0.0016	+ 0.0014	- 0.0014	0.015	0.0015	+ 0.0014	- 0.0013	0.374	0.0	+ 0.0049	- 0.0049	0.049	0.0	+ 0.0034	- 0.0032
J0004+4615	0.068	0.0121	+ 0.0129	- 0.0116	0.049	0.0021	+ 0.0025	- 0.0025	0.0	0.0044	+ 0.0038	- 0.0	0.207	0.018	+ 0.0041	- 0.0041	0.041	0.0134	+ 0.0017	- 0.0017
J0005+0524	0.997	0.0026	+ 0.01	- 0.01	0.107	0.0006	+ 0.0011	- 0.0011	0.007	0.0005	+ 0.0006	- 0.0006	0.172	0.0	+ 0.0103	- 0.0103	0.0	0.0	+ 0.0124	- 0.0
J0005-1648	0.526	0.0801	+ 0.0421	- 0.0421	0.131	0.0037	+ 0.0026	- 0.0026	0.014	0.0024	+ 0.0014	- 0.0013	0.159	0.0	+ 0.0016	- 0.0016	0.008	0.0	+ 0.0013	- 0.0012
J0005+3820	0.839	0.0181	+ 0.0176	- 0.0185	0.509	0.003	+ 0.0031	- 0.0031	0.053	0.0024	+ 0.0024	- 0.0023	0.698	0.0	+ 0.0049	- 0.0049	0.029	0.0	+ 0.0049	- 0.0042
J0006-0623	0.795	0.023	+ 0.0223	- 0.0239	2.243	0.0089	+ 0.0112	- 0.0112	0.141	0.0068	+ 0.0075	- 0.0071	3.413	0.0	+ 0.0819	- 0.0785	0.572	0.0	+ 0.0595	- 0.0526

Equations

$$P(S|S_{\mathrm{obs}},{\sigma }_{\mathrm{obs}})=\frac{1}{{\sigma }_{\mathrm{obs}}\sqrt{2\pi }}\exp \left[-\frac{(S-S_{\mathrm{obs}}{)}^2}{2{\sigma }_{\mathrm{obs}}^2}\right],$$

$$P(S|f_{\mathrm{t}},S_{\mathrm{off}},S_{\mathrm{on}},{\sigma }_{\mathrm{off}},{\sigma }_{\mathrm{off}})=\frac{1-f_{\mathrm{t}}}{{\sigma }_{\mathrm{off}}\sqrt{2\pi }}\exp \left[-\frac{(S-S_{\mathrm{off}}{)}^2}{2{\sigma }_{\mathrm{off}}^2}\right]$$

$$+\frac{f_{\mathrm{t}}}{{\sigma }_{\mathrm{on}}\sqrt{2\pi }}\exp \left[-\frac{(S-S_{\mathrm{on}}{)}^2}{2{\sigma }_{\mathrm{on}}^2}\right],$$

$$l_{\mathrm{obs}}$$

$$\mathcal{L}=\prod _{j=1}^N{l}_{\mathrm{obs},\mathrm{j}}\Rightarrow \log \mathcal{L}=\sum _{j=1}^N\log l_{\mathrm{obs},\mathrm{j}}.$$

$$⟨S⟩=(1-f_{\mathrm{t}})S_{\mathrm{off}}+f_t{S}_{\mathrm{on}},$$

$$\mathrm{Var}=[(1-f_{\mathrm{t}})({\sigma }_{\mathrm{off}}^2+S_{\mathrm{off}}^2)+f_t({\sigma }_{\mathrm{on}}^2+S_{\mathrm{on}}^2)]-⟨S{⟩}^2.$$

$$l_i$$

Full text

1 Introduction

BL Lac objects (BL Lacs) and Flat Spectrum Radio Quasars (FSRQs) constitute a sub-class of active galactic nuclei (AGN) called blazars. Blazars are known for their powerful and highly relativistic jets, which are pointed close to our line of sight (Readhead et al., 1978; Blandford & Königl, 1979; Scheuer & Readhead, 1979; Readhead, 1980). Due to the alignment of the jet, their emission, from radio to the highest energy $\gamma$-rays, is dominated by relativistic effects such as boosting of the observed flux and compression of timescales. They show a compact one-sided core-dominated jet morphology, as well as apparent superluminal motion of the radio jet components propagating downstream of the core as seen through VLBI observations (Readhead et al., 1978; Scheuer & Readhead, 1979; Readhead, 1980). Such radio components have shown superluminal motions of up to about $50c$ (Lister et al., 2009, 2013), while they have also been associated with numerous other phenomena such as the production of $\gamma$-rays and rotations of the optical polarization plane (Marscher et al., 2008, 2010).

Variability in the radio regime was one of the first identified characteristics of blazars (Dent, 1965). Although there are cases of quasi-periodicity (Carrasco et al., 1985; Valtaoja et al., 1985; King et al., 2013), generally the variability is erratic. However, it cannot be fully explained by a stochastic process, due to the appearance of outbursts, often in several frequency bands simultaneously (Aller et al., 1999; Hovatta et al., 2008; Max-Moerbeck et al., 2014; Fuhrmann et al., 2014; Blinov et al., 2015; Angelakis et al., 2016; Hovatta et al., 2016), followed by periods of relatively low activity. It has been suggested that such outbursts can be characterized by an exponential rise and decay with the ratio of the respective timescales to be approximately 1.3 (Valtaoja et al., 1999; Lähteenmäki et al., 1999; Lähteenmäki & Valtaoja, 1999; Hovatta et al., 2009). The timescale and flux-density amplitude of an outburst can vary from days to months and from comparable, to orders of magnitude higher than the quiescent one (Aller et al., 1999). Thus the radio flux-density curve of a blazar can generally be described by a quiescent level (the minimum radio output of the source), and a series of consecutive aperiodic outbursts on top of that minimum output.

Modeling the variability of blazars has been the subject of several studies. Valtaoja et al. (1988) attempted to separate the quiescent from the flaring flux through multi-wavelength flux-density curves (several frequencies from 4.8-90 GHz). The authors examined the spectrum of each source at periods of minimum flux between outbursts (what they considered as the “constant” flux of the jet), and subtracted it from the source spectrum during an outburst in order to obtain the “variable” flux. Lister (2001) chose to model the blazar flux-density curves (4.8 and 5 GHz) with a shot-noise process (Lister, 2001, and references therein). Using Poisson statistics [and the exponential profile for the flares described in Valtaoja et al. (1999); Lähteenmäki et al. (1999); Lähteenmäki & Valtaoja (1999)] the authors investigated the variability duty cycles and sample selection biases due to variability in blazars.

The interest in the subject is well-motivated: understanding the general variability properties of blazars in different frequencies can provide important information on their emission mechanisms as well as the location of the emission region. Although rest-frame time delays between frequencies are necessary to constrain the location of emission regions on a source-by-source basis (Fuhrmann et al., 2014; Max-Moerbeck et al., 2014), differences in the variability properties at different frequency bands could be indicative of the spatial connection of their respective emission regions on a population level.

The OVRO 40m monitoring program (Richards et al., 2011; Richards et al., 2014) provides a unique opportunity for studies of blazar radio variability, thanks to its unprecedented sample size and cadence (about 1800 sources observed twice weekly on average over 8 years). OVRO data from 2 years of monitoring were used in Richards et al. (2011) to examine blazar variability and amplitudes using a likelihood approach. Modeling each source’s flux-density distribution by a single Gaussian distribution, the “intrinsic” mean flux-density and modulation indices (what one would have observed in the limit of infinite accuracy and sampling) were estimated. The authors used these results to uncover a statistically significant discrepancy between the radio variability (as quantified by the intrinsic modulation index) of $\gamma$-ray loud and $\gamma$-ray quiet sources that were otherwise similar. For this reason, the radio intrinsic modulation index has since been used to select $\gamma$-ray quiet sources that are as similar in their radio properties as possible to $\gamma$-ray loud blazars, for monitoring in other frequencies (e.g. Pavlidou et al., 2014; Angelakis et al., 2016).

The likelihood formalism used by Richards et al. (2011) provides a robust way to account for observational uncertainties and finite cadence as well as calculate uncertainties for the estimated quantities, but at the expense of model-independence. In this case, model-dependence enters through the assumption that a single-Gaussian is a good description of the distribution of flux-densities. If that assumption is not valid for some sources, this can affect the results in two ways. First, any estimated quantities that are sensitive to the assumed underlying family of flux-density distributions may exhibit systematic offsets. Most importantly, these offsets could not be accounted for by the calculated uncertainties, which are statistical and based on the assumption that the underlying distribution model holds (see e.g., Mouschovias & Tassis, 2010). Second, important information regarding the detailed behavior of the source (such as bimodality in the flux-density distribution) is lost through the simplified treatment.

As discussed in Richards et al. (2011), some of the OVRO 40m-monitored blazars show a flux-density distribution that can be well described by a single Gaussian, while others show a bimodal distribution. To remedy the first of the two problems described above (i.e. cases where a single Gaussian model is not a good description of the flux-density distribution) Richards et al. (2011) used, and reported, only quantities that their method estimates robustly. Such quantities are the mean flux and the modulation index, in contrast, for example with the most likely flux of a source. A typical example of such a source is BL Lac, the flux distribution of which (using the full 8 years OVRO dataset) is shown in Fig. 1. The maximum-likelihood single Gaussian for this data set is overplotted with the dashed line. As BL Lac is very bright, and the observational uncertainties are very small compared to the typical flux ($<2\%$), we do not expect the uncertainties to widen the flux distribution significantly, or, conversely, the maximum-likelihood distribution to be appreciably narrower than the data. Indeed, this is the case in Fig. 1. The source exhibits bimodality. For this reason, the maximum-likelihood single Gaussian is not a good description of the underlying flux-density distribution. As a result, while the single-Gaussian mean and spread are clearly reasonable representations of the corresponding parameters of the data distribution, the most-likely flux is missed by the single Gaussian, and information such as the flaring ratio (ratio of the mean “on”-state flux-density to “off”-state flux density $R=S_{\rm on}/S_{\rm off}$) is lost. On the other hand, a model that can account for the bimodality of the sources can recover this important information. Such a model is particularly important given the large size, good sampling, low noise and statistical completeness of the dataset at hand.

In this work, we extend the formalism of Richards et al. (2011) in order to model the observed blazar variability as an alternation between a low-flux state and a high-flux state with the possibility of different degrees of variations about the mean flux of each state. The physical interpretation of such a model is that blazars spend a fraction of their time in a state characterized by relatively low or no activity, and the remaining time in a state of increased variability and outbursts. This way we attempt to separate the quiescent and flaring components of the flux-density distribution of each source, and examine the general variability characteristics, as well as differences between different blazar subclasses. The paper is organized as follows. In section 2 we present the maximum likelihood formalism we developed in order to model the blazar flux-density distributions. In section 3 we present the sample used in this work. In section 4 we compare the results from our model with that of a single Gaussian distribution and test the bimodality of the flux-density distributions of the sources in our sample. In section 5 we derive best-fit distributions of all quantities characterizing bimodal variability for the blazar population and we compare the variability characteristics of different blazar subsamples. Finally in section 6 we summarize and discuss our findings and conclusions.

The cosmological parameters we adopt in this work are $H_{0}=71$ ${\rm km\,s^{-1}\,Mpc^{-1}}$, $\Omega_{m}=0.27$ and $\Omega_{\Lambda}=1-\Omega_{m}$ (Komatsu et al., 2009).

2 Maximum likelihood analysis

At any time of observation the flux-density ($S$) emitted by a source, given the observational uncertainty (which we assume to be Gaussian), can be described as

[TABLE]

where $S_{\rm obs}$ is the observed flux-density, and $\sigma_{\rm obs}$ is the error on the measurement. We model the flux-density distribution of blazars as a sequence of two states dubbed “off” and “on”. Each state can be assumed to follow a Gaussian distribution similarly to Richards et al. (2011). The “off”-state describes the period of time a blazar spends in quiescence, while the “on”-state the period of time a blazar spends flaring. Then the probability density function of the blazar emitting a certain flux-density will be given by,

[TABLE]

where $f_{\rm t}$ is the duty cycle, $S_{\rm off}$ is the mean flux-density of the “off”-state, $\sigma_{\rm off}$ is the standard deviation of the mean flux-density of the “off”-state, $S_{\rm on}$ is the mean flux-density of the “on”-state, and $\sigma_{\rm on}$ is the standard deviation of the mean flux-density of the “on”-state.

Marginalizing over all possible flux-densities S (limits of the integral from $-\infty$ to $\infty$, see Venters & Pavlidou (2007); Richards et al. (2011) for a detailed derivation of the integral), the likelihood of observing $f_{\rm t}$, $S_{\rm off}$, $S_{\rm on}$, $\sigma_{\rm off}$, and $\sigma_{\rm on}$ given a flux-density $S_{\rm obs}$ and error $\sigma_{\rm obs}$ is,

[TABLE]

The first term of Eq. 3 describes the “off”-state, while the second the “on”-state. For a number of observations $N$ ($j=1,2,3...N$), the likelihood function is,

[TABLE]

We use the Simplex algorithm (also known as Nelder-Mead, Nelder & Mead, 1965) implemented in the scipy.optimize.minimize python package to minimize the negative log-likelihood, and thus to maximize the likelihood. After several trials with different methods and minimizing routines implemented in various programming environments, the Simplex algorithm was proven to be the most consistent and stable, given the task at hand.

Given the complexity of the likelihood function, in order to avoid numerical instabilities and possible sensitivity to initial conditions, we perform the minimization $10^{3}$ times, each time with random initial conditions, and choose as the best-fit parameter estimates those that gave the minimum function value. Figure 1 shows the maximum likelihood fit of our method to the 15 GHz flux-density distribution of BL Lacertae while Fig. 2 shows the OVRO light curve. The horizontal lines mark the mean flux-densities of the “off”- and “on”-states and the grey shaded areas their 1-standard deviation.

To derive the error on our estimated parameter values, we use the Fisher information matrix. The information matrix gives a measure of the amount of information every parameter carries on the curvature of the likelihood at the best-fit values. From the information matrix, we calculate the variance-covariance matrix. The error on each parameter estimate is the square root of the corresponding diagonal element of the variance-covariance matrix. It is possible that the matrix is not positive definite. Since the likelihood is multi-parametric, this would suggest that two or more parameters are anti-correlated. In such a case, we estimate the error of that parameter using a slice of the likelihood surface parallel to the axis of the parameter of interest, with the values of the other parameters set at their best-fit values (i.e., a slice passing through the maximum likelihood point). Using that slice, we determine the values of the parameters that reduce the likelihood by factor $e^{-1/2}$ (in a Gaussian slice, this would be the $\pm 1\sigma$ points). These two values set the 1$\sigma$ uncertainty on that parameter.

The method described above can provide robust results as long as the number of observations is sufficiently large. In the application of our method (see section 3) the number of observations per source is 421 on average with a standard deviation of 89 (minimum 217 and maximum 1108).

3 Sample

We apply our method to data from the Owens Valley Radio Observatory’s (OVRO) 15 GHz blazar monitoring program111http://www.astro.caltech.edu/ovroblazars/ (Richards et al., 2011). OVRO has been monitoring a sample of $>$1800 blazars in support of the Fermi Gamma-ray Space Telescope (Fermi, Acero et al., 2015) with an approximate cadence of twice per week since 2007. The monitoring program began using the 1158 sources from the Candidate Gamma-ray Blazar Survey (CGRaBS) complete sample (Healey et al., 2008). Since then, however, to facilitate monitoring of additional sources detected by Fermi, the sample has increased by approximately 1/2 of its original size. In order to maintain the statistical integrity on the sample under investigation while taking advantage of the additional $\gamma$-ray sources that were added later, we construct a new statistically complete sample (a flux-density–limited one), as follows: we use the maximum likelihood mean flux-density from Richards et al. (2014) to plot the source count distribution for all OVRO monitored sources (Fig. 3). The flux limit is set at 0.354 Jy which is where the distribution shape starts to deviate from a power law. The final sample consists of 584 sources; 435 FSRQs, 81 BL Lacs, and 68 other sources (which we dubbed U-R), 17 of which are classified as radio galaxies and 51 as possible blazar or unidentified. We follow the classification of Richards et al. (2011); Richards et al. (2014) (and references therein) i.e. objects with broad emission lines in optical are classified as FSRQs whereas objects with weak or even no emission lines are classified as BL Lacs. We have included data from 05/01/2008 to 14/02/2016.

4 Single Gaussian versus bimodal

Using the method described in section 2 we fit the flux-density distribution of every source in the flux-limited subsample defined above, in order to extract their variability properties. It is not necessary that the flux-density distribution of a source is described well by a “off”-“on” model. It is very well possible that a source resides in either state for the entire duration of the monitoring period, which is roughly 8 years. In this case, the maximum likelihood best-fit model is degenerate. Models with duty cycle $f_{\rm t}=1$, $f_{\rm t}=0$, or two identical states with $f_{\rm t}=0.5$ are mathematically possible and describe the same single Gaussian distribution. For this reason, in this section we compare the results of our model and the single Gaussian model as in Richards et al. (2011) (fit to the 8-year dataset) in order to understand which model best describes the data, and identify all sources with such degeneracy. For the purposes of our analysis we use the Kolmogorov-Smirnov test (K-S test) and the Wilcoxon rank-sum test (Wilcoxon test). The K-S test gives the probability of two samples being drawn from the same distribution, while the Wilcoxon test gives the same probability but with the alternative hypothesis that the values of one sample are systematically larger than those of the other. For any probability higher than 5% we accept that neither test can reject the null hypothesis that the two samples are drawn from the same distribution. Only 52 sources ($8.9\%$) can be well-described by a single-Gaussian model. For the purposes of the population studies (see section 5) we have excluded these sources from our analysis. The excluded sources also account for cases where the source has no significant variability, since the flux-density distribution of such a source would be consistent with a single narrow Gaussian distribution.

We next proceed to a comparison between the overall mean flux-density and modulation index (in contrast to the mean and modulation index of the individual variability states) of the bimodal model with the corresponding qualities of the single-Gaussian model, in order to evaluate any systematic effects on these quantities due to the single-Gaussian assumption.

From the definition of the mean and variance, the mean flux-density of the bimodal distribution is equal to,

[TABLE]

and the variance is

[TABLE]

Using Eq. 5 and 6 the overall modulation index is equal to $\mathrm{\langle m\rangle=\sqrt{\rm Var}/\langle S\rangle}$. We calculate the overall mean flux-density and intrinsic modulation index for each of our sources, and compare it with the corresponding values from a single Gaussian model from Richards et al. (2011). In order to estimate uncertainties on $\mathrm{\langle S\rangle}$ and $\mathrm{\langle m\rangle}$ we use the uncertainties estimated by the likelihood analysis for ${\rm f_{t}}$, $S_{\rm on}$, $S_{\rm off}$, $\sigma_{\rm on}$, and $\sigma_{\rm off}$, and standard error propagation.

For the mean flux-density the mean of the two samples is the same with both tests, while both K-S and Wilcoxon tests allow for the hypothesis that the samples are drawn from the same distribution (Table 1). Figure 4 shows the distribution of the overall $\mathrm{\langle S\rangle}$ for the single Gaussian (solid black line) and bimodal (dashed green line) models and Fig. 5 shows the comparison of the two models for individual sources (upper panel) and the distribution of the fractional difference $\mathrm{(\langle S_{\rm Gaussian}\rangle-\langle S_{\rm Bimodal}\rangle)/\langle S_{\rm Bimodal}\rangle}$, between models (lower panel). It is clear that the mean is very robust against the single-Gaussian assumption. The distribution of fractional differences is symmetric about zero, indicating that there is no bias in the mean introduced by the single-Gaussian assumption.

For the modulation index, Fig. 6 shows the distribution of the overall $\mathrm{\langle m\rangle}$ for both models and Fig. 7 shows the comparison of the two models for individual sources (upper panel) and the distribution of their fractional difference (lower panel). Both K-S and Wilcoxon test indicate disagreement between the two estimates, although not at extremely high significance (Table 1). The scatter between models is larger than with the $\mathrm{\langle S\rangle}$ (Fig. 7). However, for most sources the values are consistent within uncertainties, and the fractional difference (lower panel of Fig. 7) is generally less than 20% with the exception of very few sources. For this reason, we do not expect the single-Gaussian assumption to have had a strong impact on the statistical comparison between populations in Richards et al. (2011) and Richards et al. (2014). It is nevertheless worth noting that the single-Gaussian assumption introduces a negative bias in the modulation index: the distribution of fractional difference between models is not symmetric about zero (mean$\approx$-0.09, median$\approx$-0.04) with the single-Gaussian model tending to underestimate the modulation index and the overall variability of the sources.

The discrepancy originates from the fact that, for a source well-described by a bimodal distribution (e.g. Fig. 1), a single Gaussian model would favor the dominant peak (in that case the “off”-state) and attempt to accommodate the alternate state in the tail of the distribution. That is why although the $\mathrm{\langle S\rangle}$ derived by the two models are consistent, the $\mathrm{\langle m\rangle}$ are not.

On the other hand, if we compare the overall $\mathrm{\langle m\rangle}$ from the bimodal model (maximum-likelihood estimate, MLE) with a standard sample $\mathrm{\langle m\rangle}$ (ratio between sample standard deviation and sample mean, Fig. 8) we find that they are in good agreement within uncertainties. The Wilcoxon test yields a $45\%$ and the K-S test a $87\%$ probability of consistency respectively. In order to arrive to such result, two conditions need to be met. First, the model we have assumed (bimodal) has to be a good description of the underlying distribution of the sources, and second the observational errors have to be relatively small compared to the intrinsic variability of the sources. The latter is expected for the sources in our sample, because of the two filters we have used: a relatively high flux-density limit (which ensures that fractional observational uncertainties are low), and the rejection of sources well-described by a single-Gaussian (which ensures that the remaining sources are significantly variable). To test whether this is indeed the case we perform the same analysis, but for the fainter sources observed by OVRO that were not bright enough to be in our flux-density-limited subsample. For the fainter sources ($\mathrm{\langle S\rangle}\leq 0.354~{}Jy$, the flux-density limit) the observational errors are relatively larger and will contribute more to the overall variance. We find that the fractional difference increases (Fig. 9), with both tests (Wilcoxon, K-S) rejecting the null hypothesis that the two estimates of $\mathrm{\langle m\rangle}$ are drawn from the same distribution, as expected. Thus, for the brightest and most variable sources, the ability of the maximum likelihood approach to take errors into account does not provide a significant improvement on the overall estimate of $\mathrm{\langle m\rangle}$ over the simple direct standard deviation-over-mean estimate from the data. However, even for these sources the model introduced in this work has the advantages that: (a) it allows us to separate quiescent and flaring states; and (b) it provides information on the variability of each state, flaring ratio and duty cycle on a source-by-source basis.

5 Variability properties

5.1 Overall variability properties of the sample

The distribution of the flaring duty cycle $f_{\rm t}$, the flaring ratio $R$, the intrinsic modulation index in the “off”-state, $m_{\rm off}$ and the intrinsic modulation index in the “on”-state, $m_{\rm on}$ for the entire sample, are shown in Figs. 10, 11, 12, and 13 respectively with the black solid line. Table 2 shows the sample mean and standard deviation for these quantities.

It has been shown in Richards et al. (2011) that intrinsic modulation index follows a mono-parametric exponential distribution. This appears to be true for the individual “on”- and “off”- state modulation indices as well as the flaring ratio, while the duty cycle appears to be consistent with a uniform distribution for the flux-limited subsample. Parameterizing such quantities could prove a useful tool in modeling blazars at the population level.

For the duty cycle, we test with the use of the K-S and Wilcoxon tests whether it is consistent with a uniform distribution in the [0,1] interval. Both tests cannot reject the null hypothesis that the duty cycle and a uniform distribution in the [0,1] range, are drawn from the same distribution (14% and 19% probability of consistency respectively). For this test, we have not excluded any sources based on their consistency with a single-Gaussian model.

For the flaring ratio and modulation indices, the mean of the exponential distribution ($K_{\rm 0}$, where $K$ is the parameter to be fitted) given a set of observations ($K_{\rm obs,i}$) with Gaussian uncertainty $\sigma_{\rm obs,i}$ , can be calculated using,

[TABLE]

from Richards et al. (2011), where $l_{\rm i}$ is the likelihood and erf is the error function. The error on the quantities to be fitted is calculated through integration of the normalized likelihood. For the modulation indices we exclude all values of $m<0.06$ and re-normalize the likelihood accordingly (see Eq. 30 and the corresponding discussion in Richards et al., 2011).

Table 3 shows parameters of the best-fit distributions for each of the aforementioned quantities which can be used as inputs to population models for blazars that require some treatment of source variability and/or activity state, such as the derivation of luminosity functions from single-epoch surveys. From the fits we have excluded sources J1433$-$1548, J1823+7938, J1808+4542, J1852+4019 for being outliers (either of $R$ or $m$), preventing the maximum likelihood method to achieve a good fit. In subsequent sections we focus on the comparison of variability properties between blazar subclasses.

5.2 FSRQs versus BL Lacs

Figure 10 shows the distribution of the fraction of time spent on the “on” state (the blazar flaring duty cycle) for the whole sample (black solid line), BL Lacs (red dotted line), FSRQs (green dashed line), and U-Rs (blue dashed-dotted line). The distribution is similar to a uniform one with no apparent difference between BL Lacs and FSRQs. This is also confirmed by the K-S test and the Wilcoxon test.

We next explore the flaring ratio of the flux-density of the “off” and“on” states which is unaffected by redshift effects. This way we are not limited by the redshift completeness of our sample.

Figure 11 shows the ratio of the “on”-to-“off”-state ($R$) for the different populations. The BL Lacs appear to have larger $R$ (on average) than the FSRQs (Table 4). Both tests reject the null hypothesis that the samples are drawn from the same distribution with the values of one (BL Lacs) being systematically larger than the other. This would suggest that while flaring the BL Lacs reach, on average, higher flux-densities on the on-state relatively to their off-states than the FSRQs.

Although the ratio $R$ is different for the BL Lacs and FSRQs, the individual flux-densities of the “off” (Fig. 14) and “on” states (Fig. 15) appear to be similar with the K-S ($45\%$ for the “off” and $46\%$ for the “on” state) and Wilcoxon ($83\%$ for the off and $32\%$ for the on-state) tests unable to reject the null hypothesis, suggesting that the BL Lacs exhibit relatively larger outbursts than the FSRQs.

We also explore the characteristic variability of the populations using the intrinsic modulation index defined as $m=\sigma/S$ (Richards et al., 2011). Figures 12 and 13 show the distribution of $m$ for the “off” and “on” states respectively.

For the “off”-state the two tests cannot reject the null hypothesis that the two samples are drawn from the same distribution. However, for the “on”-state the mean for the BL Lacs is significantly larger than the FSRQs both tests rejecting the hypothesis that the two samples are drawn from the same distribution. The fact that the mean “on”-state modulation index value is larger for the BL Lacs suggests that while flaring their flux-density distribution is on average wider. On the other hand, for the FSRQs, there is no difference of the intrinsic modulation index between the two states. We conclude that the BL Lacs are relatively more variable during outbursts than the FSRQs. All the mean parameter and probability values are summarized in Table 4.

5.3 Blazars versus Unidentified

Comparing the blazar sample (FSRQs + BL Lacs) against the U-R sources we find very interesting similarities. The duty cycle is similar, with neither test rejecting the null hypothesis that the samples are drawn from the same distribution (Table 5). The mean flaring ratio $R$ is mildly (but significantly) larger for the blazars than for the U-Rs, in both states the blazars have larger mean modulation index, yet both tests are unable to reject the null hypothesis that the two samples are drawn from the same distribution (Table 5).

We also compare individual classes (i.e., FSRQs , BL Lacs) versus U-Rs. Comparing the FSRQs none of the above results changes in any significant manner, which is to be expected since the FSRQs dominate the blazar sample (75% of the sample). For the BL Lacs the results are somewhat different. The duty cycle as well as the “off”-state modulation index are similar, with both tests not being able to reject the null hypothesis that the two samples are drawn from the same distribution. However, both tests reject the null hypothesis for the flaring ratio $R$ and for the “on”-state modulation index (Table 6).

If the unidentified and blazar candidate sources in the U-R sample were unidentified BL Lacs due to the absence of spectral lines, we would expect the opposite results, i.e., the BL Lac sample to be more consistent with the U-R sample in the variability characteristics. Instead the FSRQs appear to be consistent with the U-Rs suggesting that the majority of the sources in the U-R sample are either unidentified FSRQs or radio galaxies with jets pointed close to our line of sight, yet not close enough to be considered a blazar (often referred to as misaligned blazars). The latter would suggest that the FSRQs and the radio galaxies share similar variability characteristics. All the mean parameter and probability values are summarized in Table 5.

5.4 Fermi-detected versus Fermi non-detected blazars

There are 267 sources in our flux-density-limited sample (50%) that have been associated with Fermi detected sources. Out of the 267 sources, 186 are FSRQs, 62 are BL Lacs, and 19 are U-Rs. Figure 16 shows the distribution of the duty cycle. The solid black line is for the Fermi-detected and the dashed green for the Fermi-non-detected sources. The distributions are very similar in shape with the Wilcoxon and K-S test suggesting the two samples being drawn from the same distribution (Table 7). The high probability of consistency, as well as the almost identical mean values suggest that there is no difference in the time the two samples spent in either state.

However, their flaring ratio $R$ (Fig. 17) is rather different. Both tests rejected the null hypothesis of consistency between populations with high confidence (i.e., small probability value). The mean values combined with the results of the Wilcoxon test suggest that the Fermi-detected sources have systematically higher ratios, which translates to them having larger differences in flux-density between states, i.e., more powerful outbursts.

Examining the intrinsic modulation index in the different states (Fig. 18, 19), we find that in both the “off” and the “on” state, both tests rejected the null hypothesis. In both states the Fermi-detected sources have larger mean values which means that on average they are more variable than the Fermi-non-detected sources. At the same time there is no difference in the variability of the two states in the Fermi-non-detected sources, contrary to the Fermi-detected. This could partially be attributed to the fact that the majority of the BL Lacs (which are more variable in the “on” state, see section 5.2) are in the Fermi-detected subsample, however, the fact that the Fermi-detected are systematically more variable than the Fermi-non-detected in both states, is in support of the Fermi-detected being intrinsically more variable than the Fermi-not-detected sources (see also Richards et al. (2014)). All the mean parameter and probability values are summarized in Table 7.

5.5 Low versus High redshift sources

Richards et al. (2011) and Richards et al. (2014) point out that FSRQs show a negative correlation between radio variability as quantified by their intrinsic modulation index and redshift ($z$). They find that sources with $z<1$ are more variable in radio than sources with $z>1$ suggesting that the radio variability of FSRQs evolves with cosmic time. The trend first appeared in the two year dataset (Richards et al., 2011), and persisted in the four year dataset (Richards et al., 2014). We use the redshift values from Richards et al. (2014) to separate our sources into two subsamples according to their redshift. All the FSRQs in our sample have a known redshift. Separating our samples in low ($z<1$) and high ($z>1$) redshift sources we have 149 and 252 sources respectively.

For the duty cycle and the intrinsic modulation index in the “on”-state, both tests agree that the two subsamples are drawn from the same distributions. For the flaring ratio and the “off”-state modulation index the tests reject the null hypothesis. If instead of $z=1$ we separate our samples according to the mean ($z\approx 1.22$) or median ($z\approx 1.32$) value the results of both test vary for all parameters. Sensitivity to redshift separation makes the individual components of the analysis ($f_{\rm t}$, $R$, $m_{\rm off}$, $m_{\rm on}$) unreliable indicators for cosmic evolution. However, if we compare the overall modulation index ($\mathrm{\langle m\rangle}$) we find that regardless of redshift separation the two samples are inconsistent. Separating the samples at $z=1$ and $z=\langle z\rangle_{\rm median}$ the probability values of both tests are $\leq 10^{-3}$. If we separate the samples at $z=\langle z\rangle_{\rm mean}$ the Wilcoxon test yields a $2\%$, whereas the K-S test yields a $6\%$ probability of consistency which is marginally acceptable. Since the mean $\langle m\rangle$ values for the two samples are different ($\langle m\rangle_{\rm mean}=0.185$ for the low and $\langle m\rangle_{\rm mean}=0.177$ for the high redshift sources) we conclude that the two samples are not drawn from the same distribution. On the other hand, if we correct for the cosmological time dilation ($\Delta t_{\rm obs}=(1+z)\Delta t_{\rm rest-frame}$) and repeat the analysis, but this time using equal redshift-corrected observing lengths for all sources, we find no apparent trend of cosmic evolution. For all the parameters both tests cannot reject the null hypothesis that the two samples are drawn from the same distribution. However, we have not taken into account relativistic compression of variability timescales, which should be larger for higher-z sources (because of flux-density selection effects) and thus affect results in the opposite direction to cosmological time dilation. Given the large span of Doppler factors in blazar jets (Hovatta et al., 2009; Liodakis et al., 2017) even the redshift-corrected observing lengths will be significantly different from the “true” jet rest-frame observing lengths for each source. Until a large enough number of Doppler factor estimates is available that will allow us to confidently correct for the relativistic effects on a source-by-source basis, we are unable to drawn firm conclusions on the cosmic evolution of FSRQ variability properties.

A similar effect as the one discussed here and in Richards et al. (2011); Richards et al. (2014) is true for interstellar scintillation at 5 GHz (Lovell et al., 2008; Koay et al., 2012). The authors found that interstellar scintillation is suppressed a higher redshift and that high redshift sources have steeper spectra (in the 5-9 GHz range). That the interstellar scintillation of high redshift sources is suppressed indicates either a larger apparent angular size, beyond the expected cosmological expansion, or a smaller compact fraction in the high redshift blazars. However, since the OVRO dataset is single frequency we cannot test if a similar effect is present in our sample.

6 Summary

We have presented a novel five-dimensional maximum likelihood formalism in order to characterize the variability properties of blazars and blazar-like sources as a series of two states: an “off” state describing the low activity periods of a source, and an “on” state that describes periods of outbursts. We used our method to fit the 15 GHz flux-density distribution of blazars as seen by OVRO (Richards et al., 2011), and extract their variability properties (Table 8).

For our statistical analysis of the population properties of the OVRO blazars, as well as for the comparison of the behaviors of blazar subclasses, we define a statistically complete subset of all OVRO monitored sources: a flux-density limited subsample, based on average fluxes from Richards et al. (2014). We have tested whether the model assumed in this work is a good description of the underlying flux-density distribution of the OVRO-monitored blazars. If this is the case, then for bright sources (where observational uncertainties do not widen appreciably the observed flux density distribution compared to the intrinsic one) we would expect the value of $m$ obtained from the sample mean and sample standard deviation for each light curve to agree well with our likelihood-derived intrinsic $\mathrm{\langle m\rangle}$. We have verified that the two are, in fact, in excellent agreement.

Having established that our model is a good description of the underlying distribution, we compared the variability characteristics of different subsamples. Our results can be summarized as follows:

• •

BL Lacs are more variable than the FSRQs. This is consistent with the finding of Richards et al. (2011); Richards et al. (2014) that BL Lacs have a higher overall $\mathrm{\langle m\rangle}$. However, we have now established that BL Lacs also exhibit stronger outbursts (have a higher flaring ratio), and that their increase in $\mathrm{\langle m\rangle}$ is dominated by the “on”-state (it is $m_{\rm on}$ that is significantly higher, while $m_{\rm off}$ is similar in BL Lacs and FSRQs). Interestingly, Liodakis & Pavlidou (2015); Liodakis et al. (2017) find that the Doppler factors of FSRQs are on average, significantly higher than those of BL Lacs, so this discrepancy must have its origin in rest-frame properties, rather than differences in boosting between the two classes.

• •

Sources classified as blazars (BL Lacs and FSRQs) have systematically larger flaring ratios (i.e., stronger outbursts) than the U-Rs with otherwise similar variability characteristics. The variability characteristics of U-Rs are similar to FSRQs. Since the majority of U-Rs are blazar candidates, it would suggest either that they are unclassified FSRQs most likely due to lack of multi-wavelength observations, or that they are unclassified radio galaxies, which in turn would suggest that FSRQs and radio galaxies share similar rest-frame variability characteristics.

• •

Fermi-detected sources are intrinsically more variable than the Fermi-non-detected sources. This result agrees with the overall findings of Richards et al. (2011); Richards et al. (2014), however the bimodal model offers an opportunity to trace the origin of this result in the details of the behavior of blazars in the flux density domain. Indeed, Fermi-detected sources have higher flaring ratios and higher modulation indices in both states, with the most significant difference being in the flaring ratio. This results indicates that the mechanisms responsible for the amplitude of radio variations and the $\gamma$-ray loudness of a source may share a physical link.

• •

The overall intrinsic modulation index ($\langle m\rangle$) is consistently (regardless of redshift separation of subsamples) supporting the negative correlation between radio variability and redshift in FSRQs reported in Richards et al. (2011); Richards et al. (2014). Once we accounted for the cosmological time dilation we found no evidence for such negative correlation. However, since we are not yet able to take into account the relativistic effects compressing blazar timescales, we are not able to come to firm conclusions regarding a possible cosmic evolution of variability properties of FSRQs.

We caution the reader that for an analysis such as the one presented here, there is a dependence of the derived variability parameters to the length of the monitoring program. Short (in time) monitoring programs may not be able to sample the entire flux-density distribution of a blazar. However, given enough time, estimates will converge to their “true” values. For our well sampled sources, we find that by splitting the monitoring period in half (4 years) the difference in the derived variability estimates from the two periods is $<50\%$ in the majority of cases.

The overall radio modulation index as calculated by Richards et al. (2011) is one of the properties used to select samples for monitoring of other blazar properties (such as their optopolarimetric behavior e.g., King et al., 2014; Pavlidou et al., 2014). Given that we find the bimodal flux distribution to be a much better description than the single Gaussian used by Richards et al. (2011) we would advice using the overall $\mathrm{\langle m\rangle}$ values from this work for sample characterization.

The tools presented in this work were used to explore the variability properties of a 15 GHz selected flux-density-limited sample. They are, however, not restricted to any particular frequency since the formalism is based on statistics alone and is independent of any emission mechanism or other physical arguments. Thus, are suitable for all wavelengths and sources that can be well described by a bimodal Gaussian distribution. However, one should keep in mind that using a large number of observations for the fitting (in our case 421 on average) is critical to ensure robust results for the estimated parameters and their errors.

Acknowledgments

The authors would like to thank the referee Hayley Bignall and Vassilis Karamanavis for constructive comments and suggestions that helped improve this work. I.L. thanks the Caltech Astronomy Department for their hospitality during the completion of this work. This research was supported by the European Commission Seventh Framework Program (FP7) through grants PCIG10-GA-2011-304001 “JetPop” and PIRSES-GA-2012-31578 “EuroCal”. T.H. was supported in part by the Academy of Finland project number 267324. This research has made use of data from the OVRO 40-m monitoring program (Richards et al., 2011) which is supported in part by NASA grants NNX08AW31G, NNX11A043G, and NNX14AQ89G and NSF grants AST-0808050 and AST-1109911.