Constraining the Limiting Brightness Temperature and Doppler Factors for the Largest Sample of Radio-bright Blazars

Blazar jets are known to show extremely fast variability, boosted emission, and apparent superluminal motion of jet components. These, as well as the other unique features seen in blazars, are due to the relativistic effects dominating the emission from the jet. The relativistic effects arise from the preferential orientation of the jet typically within <20° from our line of sight (Readhead et al. 1978; Blandford & Königl 1979; Scheuer & Readhead 1979; Readhead 1980). Radio emission in blazars is produced by relativistic electrons accelerated in the magnetic field of the jet-emitting synchrotron radiation. It is characterized by a flat spectrum from the centimeter and up to, in some cases, millimeter wavelengths thought to be the result of the superposition of multiple synchrotron self-absorbed spectra. Although the radio emission is considered to be less variable than, e.g., optical or γ-rays, it can still be exceptionally variable and bright, with some sources reaching radio core brightness temperatures of >10¹³ K (e.g., Kovalev et al. 2016). Since the intrinsic brightness temperature of a jet is expected to be of the order of ∼5 × 10¹⁰ K (Readhead 1994), this would suggest that the jets continue to be highly relativistic on very large scales far from the supermassive black hole. Quantifying the beaming properties of the jets is then necessary in order to understand their energetics at large scales. These relativistic effects are quantified by the Doppler factor (δ), which is a function of the velocity of the jet and the angle to the line of sight $\delta ={[{\rm{\Gamma }}(1-\beta \cos \theta )]}^{-1}$, where Γ is the Lorentz factor (${\rm{\Gamma }}=1/\sqrt{1-{\beta }^{2}}$), β is the velocity of the jet in units of the speed of light (β = u_j/c), and θ is the viewing angle. The Doppler factor, although a crucial parameter in the blazar paradigm dictating all of the observed properties of blazars, is notoriously difficult to estimate since there is no direct method to measure either β or θ. For this reason, many indirect methods have been proposed in order to estimate δ, which usually involve different energetic (e.g., Ghisellini et al. 1993; Mattox et al. 1993; Fan et al. 2013, 2014) and/or causality arguments (e.g., Lähteenmäki & Valtaoja 1999; Jorstad et al. 2005, 2017; Hovatta et al. 2009) or fitting the spectral energy distribution (SED; e.g., Ghisellini et al. 2014; Chen 2018) of γ-ray emitting blazars.

However, different methods often yield discrepant results, due to either assumptions that do not hold or the incorrect application of the methods (see e.g., Liodakis et al. 2017c). Liodakis & Pavlidou (2015b), using blazar population models (Liodakis & Pavlidou 2015a; Liodakis et al. 2017a), evaluated a number of these methods and found that the variability Doppler factor method (Lähteenmäki & Valtaoja 1999; Hovatta et al. 2009) is the most accurate and can describe both flat spectrum radio quasar (FSRQ) and BL Lac object (BL Lacs) populations. The method is based on the assumption of equipartition between the energy density of the magnetic field and the energy density of the radiating particles, achieved at the peak of prominent flares, implying a characteristic intrinsic brightness temperature (Kellermann & Pauliny-Toth 1969; Singal 1986; Readhead 1994). By comparing the intrinsic (equipartition, T_eq) to the highest observed brightness temperature, one can estimate δ. The drawback of this method is that it is limited by the cadence of the observations, which sets a limit to the fastest observed flare and consequently a limit to the observed brightness temperature (Liodakis & Pavlidou 2015b).

In order to mitigate the effects of limited cadence, Liodakis et al. (2017d) used multiwavelength radio light curves in order to identify and track the evolution of flares throughout frequencies that allowed the authors to provide constraints on the variability brightness temperature and hence the Doppler factor of 58 sources. However, the number of blazars with simultaneous multiwavelength radio light curves is extremely limited compared to single-frequency observations. Then, the only way to overcome the effects of limited cadence is through monitoring programs with sufficiently high cadence to resolve even the fastest flares in radio.

In this work, we explore the radio beaming properties of jets by analyzing the light curves of 1029 blazars and blazar-like sources using data from the Owens Valley Radio Observatory's (OVRO's) 40 m blazar monitoring program (Richards et al. 2011). We focus on constraining the equipartition brightness temperature and the variability Doppler factors for the sources in our sample. In Section 2, we present the sample and tools of the analysis; in Sections 3 and 4, we estimate the highest brightness temperature for the sources in our sample and use blazar population models to constrain T_eq. In Section 5, we estimate the variability Doppler factors, Lorentz factors, and viewing angles based on our results on T_eq, and finally, in Section 6, we discuss the findings of this work. We have assumed the standard ΛCDM cosmology with Ω_m = 0.27, Ω_Λ = 1 − Ω_m, and H₀ = 71 km s⁻¹ Mpc⁻¹ (Komatsu et al. 2009).

From the OVRO monitored sources (∼1800), we selected those that showed prominent flares at 15 GHz via visual inspection of the light curves. Our final sample consists of 837 blazars (670 FSRQs, 167 BL Lacs), 58 radio galaxies, and 134 still unclassified sources, for a total of 1029 sources. Since 2008, OVRO has been monitoring blazars in support of the Fermi γ-ray space telescope producing the most densely sampled radio light curves available to date with a cadence of about three days. The OVRO data set provides the ideal opportunity to study the flaring and beaming properties of blazars since (1) the unprecedented high cadence is able to resolve even the fastest flares in the most variable sources and (2) the light curves of most sources are sufficiently long (8–10 years) to include at least a few major events in each source.

For the analysis of the radio light curves, we used Magnetron (Huppenkothen et al. 2015). Magnetron is a Bayesian hierarchical model implemented in Python that models the light curves as a superposition of flares characterized by an exponential rise and exponential decay on top of a stochastic background. The shape of the flares is allowed to vary (the ratio of rise to decay time can be different in each flare), and the number of fitted flares in a light curve is a free parameter. Each flare is characterized by four parameters, namely position, amplitude (in Jy), rise time (in days), and skewness (sk; decay/rise time ratio). The amplitude of a flare is defined as the difference between the peak flux density and the background level (see Figure 2 in Huppenkothen et al. 2015). The priors for the flare amplitudes and rise times are exponential, while the priors for the skewness and the number of flares are uniform distributions. The mean of the prior amplitude distribution takes values between [10⁻¹⁰, 150] Jy, while the minimum and maximum of the uniform prior distribution for the number of flares is 4100. All of the prior distributions and their hyperparameters used by Magnetron are listed in Table 1 of Huppenkothen et al. (2015). In this work, contrary to Huppenkothen et al. (2015), we treat the background with a stochastic model (Ornstein–Uhlenbeck (OU) process) to account for intrinsic blazar variability not related to flaring events. The OU process is a stochastic, stationary Gauss–Markov process often used to treat active galactic nuclei (AGNs) variability (e.g., Kelly et al. 2009). The version of Magnetron used in this work includes two new parameters to parameterize this stochastic process. The first quantity is the rate of mean reversion (α_OU), which is included in the model through a parameter L as ${\alpha }_{{\rm{OU}}}=\exp (-1/L)$. The prior for L is a log-uniform distribution such that log L ∼ Uniform(0.01 ∗ T, 0.01 ∗ T + 1000), where T is the total length of the light curve. The second parameter is the volatility of the OU process, σ_OU, i.e., the average magnitude per square root of time of random Brownian fluctuations. The prior for σ_OU is also log-uniform, such that $\mathrm{log}{\sigma }_{\mathrm{OU}}\sim \mathrm{Uniform}({10}^{-3},{10}^{3})$. While previous attempts of fitting radio light curves used a constant value for the background level, using the OU process results in a varying background across the light curve. We have verified that using a different background model (such as a constant background or a simple random walk model) results in ≤10% difference in the derived brightness temperatures, thus the choice of the background model does not affect our results in any significant way. The joint posterior probability distribution for the number of flares and the parameters of all flares as well as the hyperparameters describing the distributions of flares are sampled using Diffusive Nested Sampling (Brewer et al. 2009; Brewer & Foreman-Mackey 2016),⁶ allowing for a better exploration of the parameter space. Once the code has converged to the "true" posterior distribution, it samples ∼10² sets of flare parameters. These sets are different realizations of the flares in the observed light curve taking into account the inherent uncertainty in the parameters of the flares as well as the uncertainty in the number of flares of each light curve. A more detailed description of Magnetron can be found in Huppenkothen et al. (2015), while the code is publicly available online on GitHub.⁷

Figure 1 shows the results of the light curve modeling for four sources in our sample as well as individual flares for one posterior sample in each source. All of the light curves were visually inspected to ensure that the simulated light curves are not affected by either spurious events or observational artifacts. In addition, we compared the rise times and amplitudes of the identified flares to test whether we were able to resolve, to OVRO's sensitivity, all of the flaring events. For all classes of sources, there is a lack of flares with high amplitudes and rise times close to the cadence of OVRO (∼3 days). In the case of FSRQs and unclassified sources, we detect a mild positive correlation between the rise times and amplitudes according to the Spearman correlation test (Spearman ρ ≈ 0.3, p-value <10⁻⁵ for both classes). For all sources, we find that for rise times <14 days, the majority of flares (60%) have amplitudes lower than the median amplitude of the flares in the light curve. Out of all the flares that have higher amplitudes than the median, fewer than 20% have amplitudes higher than half of the maximum amplitude in the light curve. When we consider individual populations, we find similar percentages (±5%–10%). All of the above show that OVRO's cadence allowed us to resolve all of the most significant events within the time span of the observations. A third quality test was to assess whether Magnetron is overfitting the data, i.e., needlessly increasing the number of flares in a model to account for even the lowest flux-density variations, a common problem in χ² fitting routines usually employed. We attempt to fit four sources used in the calibration of the OVRO observations where any flux-density fluctuation in the light curve is expected to be dominated by noise rather than any flaring activity. Figure 2 shows the observed and posterior sampled light curves for those calibrator sources. Although one could visually "detect" a number of apparent flares in each source (Figure 2), no more than two flares were detected by Magnetron in any given source for any given posterior sample of the light curves. This would suggest that the stochastic model for the background is able to adequately take into account the intrinsic low amplitude variability.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

For every source in our sample with an available redshift estimate (all FSRQs, ∼71% BL Lacs, ∼56% radio galaxies, and ∼42% of unclassified sources), we estimate the variability brightness temperature (T_var) using

$$\begin{eqnarray}&&{T}_{\mathrm{var}}=1.47\cdot {10}^{13}\displaystyle \frac{{\text{}}{d}_{L}^{2}{\rm{\Delta }}{S}_{\mathrm{ob}}(\nu )}{{\nu }^{2}{\text{}}{t}_{\mathrm{var}}^{2}{(1+{\text{}}z)}^{4}}\,K,\end{eqnarray} \tag{ 1 }$$

where z is the redshift, ΔS_ob(ν) the amplitude of the flare in Jy, d_L is the luminosity distance in Mpc, ν the observing frequency in GHz, and t_var the rise time of a flare in days (Liodakis et al. 2017d). We calculate T_var for every flare in a given posterior sample and find the maximum T_var since that provides the strongest constraint on T_eq. We repeat the above process for all available samples (157 models on average) and create a distribution for T_var,max for a given source. From that distribution, we calculate the median and 1σ confidence intervals, which we quote as the uncertainty on T_var,max. Figure 3 shows four examples of the maximum T_var distributions. It is possible for the distributions to be narrower or wider than the ones shown in Figure 3. The width of the distribution reflects the ability of the modeling procedure to constrain the flare parameters' posterior distributions given the data set. Thus, the size of the confidence intervals of the T_var,max distribution gives a sense of how well we can constrain T_var,max in that source. For simplicity, we refer to $\langle {T}_{\mathrm{var},\max }\rangle$ as T_var hereafter.

Zoom In Zoom Out Reset image size

Figure 4 shows the distribution of T_var for the different populations in our sample. The lowest brightness temperature (∼10⁸K) is detected in a radio galaxy (M81), while the highest (6.7 × 10¹⁵) is in an FSRQ (J0449+1121). There is only a marginal difference between the T_var distributions of the BL Lacs and FSRQs according to the Wilcoxon rank-sum test (WRS test⁸ ; p-value 0.0497) with BL Lacs having higher values on average. BL Lacs and FSRQs also show on average higher values than the unclassified sources (WRS p-values 0.005 and 0.028, respectively). No other significant difference between the distributions of the different populations was detected. It would be interesting to also separately compare the flaring properties of the sources (i.e., individually comparing flare amplitude (maximum) and flare rise time (shortest) distributions). The comparison between the different populations in our sample showed that although there is no significant difference in the flare amplitudes between BL Lacs and FSRQs (WRS p-value 0.76, median ≈0.23 Jy for both populations), BL Lac flares evolve on faster timescales (WRS p-value 0.00003, median ≈9 days compared to median ≈17 days for FSRQs). Radio galaxies and unclassified sources have on average lower flare amplitudes (WRS p-value < 0.0002, median of ≈0.18 for both classes) than blazars, but their flare rise times are comparable to BL Lacs.

Zoom In Zoom Out Reset image size

Another interesting comparison would be between γ-ray loud and γ-ray quiet sources. We separate our sample using the ROMA-BZCAT catalog of known blazars (which is based on the 1FGL and 2FGL catalogs; Massaro et al. 2009, 2015) according to whether a source has been detected by Fermi, i.e., a source showing γ-ray emission. We find that Fermi detected (382) sources have systematically higher values than the Fermi non-detected (496) sources (WRS p-value ∼10⁻¹⁸, median 1.27 × 10¹⁴ K for detected and 1.56 × 10¹³ K for non-detected sources; Figure 4 bottom panel). Additionally, we compare their flaring properties as above. The comparison showed that Fermi detected sources flare on faster timescales and have higher amplitude flares than non-detected sources (WRS p-value < 0.0001 in both cases).

In order to constrain T_eq, we use blazar population models (Liodakis et al. 2017a). The population models are optimized using only the apparent velocity and redshift distributions from the MOJAVE survey (Lister & Homan 2005) and can yield Doppler factor distributions within flux-limited samples independent of the assumption of equipartition. We define three flux-limited samples (0.5 Jy, 1 Jy, and 1.5 Jy) above the nominal flux limit of the OVRO monitoring program (0.354 Jy) using the overall mean flux density of each source (Liodakis et al. 2017b). This allows us to assess how sensitive our results are to a given flux limit. Using the population models, we generate Doppler factor distributions for BL Lacs and FSRQs for every flux limit. From Equation (1), the variability Doppler factor (δ_var) is defined as

$$\begin{eqnarray}&&{\delta }_{\mathrm{var}}=(1+z)\sqrt[3]{\displaystyle \frac{{T}_{\mathrm{var}}}{{T}_{\mathrm{eq}}}}.\end{eqnarray} \tag{ 2 }$$

We assume that T_eq has a known distribution. We construct the T_var distribution of the sample under consideration using the estimated T_var for each source in that sample. We then use Equation (2) to derive an observed Doppler factor distribution. Then, we constrain the parameters of the assumed T_eq distribution by minimizing the reduced χ² between the expected (population model) and observed Doppler factor distributions. For the distribution of T_eq, we tested a delta function, and normal, log-normal, and uniform distributions with a parameter space [10¹⁰, 10¹³ K]. Once the best-fit parameters of each distribution were determined, we used the Bayesian Information Criterion (BIC) to select the most suitable model for T_eq.

For FSRQs, we find that the best model for T_eq is a normal distribution for all three flux-limited samples we considered with very similar means (μ) and standard deviations (σ; Table 1). All of the other distributions that were tested (although they yielded worse models according to BIC) converged to the same range of T_eq values. Although the results of the minimization for the different flux limits are consistent, the 1.5 Jy sample is the flux limit to which the population models have been optimized and thus where their strength lies (see discussion in Liodakis & Pavlidou 2015a). For this reason, we adopt the results from the 1.5 Jy sample for the FSRQs.

For BL Lacs, we also find that the best-fit distribution is normal for all flux limits; however, the parameters of the inferred distributions all significantly exceed the inverse-Compton catastrophe limit (10¹²K; Kellermann & Pauliny-Toth 1969). Since we have yet to observe the extreme behavior predicted by the inverse-Compton catastrophe, such a scenario is unlikely. A possible explanation is that the maximum Doppler factor inferred for BL Lacs by the population models is δ ≈ 30 (the maximum δ in FSRQs is δ ≈ 60), which, given the high variability brightness temperatures seen in BL Lacs, forces the very high T_eq. Liodakis et al. (2017a) discussed that the BL Lac population (∼16 sources) in the MOJAVE 1.5 Jy flux-limited sample, to which the population models are optimized, might not be a representative sample of BL Lacs, but rather a biased subsample of the brightest BL Lacs at 15 GHz. Hence, the population models cannot adequately describe the entirety of the BL Lac population present in our sample. Given that equipartition is determined by the jet processes and synchrotron physics, we expect the value of T_eq to be fairly similar for the different supermassive black-hole-powered jets. Thus, we adopt the results from the FSRQs for all populations in our sample ($\langle {T}_{\mathrm{eq}}\rangle =2.78\times {10}^{11}\,{\rm{K}}\pm 26 \%$).

In order to calculate δ_var, we draw a random value from the T_var,max distribution of each source and a random value for T_eq from a Gaussian distribution with mean $\langle {T}_{\mathrm{eq}}\rangle =2.78\times {10}^{11}\,{\rm{K}}$ and standard deviation ${\sigma }_{{T}_{\mathrm{eq}}}=7.2\times {10}^{10}$. Using Equation (2), we calculate a δ_var. By repeating this process 10³ times, we create a distribution of δ_var for every source. From the resulting δ_var distribution of each source, we estimate the median and 1σ confidence intervals. Figure 5 shows the distribution of δ_var for the different populations (top panel) and Fermi detected and non-detected sources (bottom panel). BL Lacs and FSRQs have median values of δ_var ≈ 10 and δ_var ≈ 11, respectively, while radio galaxies and unclassified sources have median δ_var ≈ 5. As expected, blazars have systematically higher Doppler factors than radio galaxies (WRS p-value ∼0.03) and unclassified sources (p-value <0.0002). We find no significant difference between the δ_var distributions of the blazar classes (BL Lacs and FSRQs; WRS p-value 0.08). Comparing Fermi detected and non-detected sources, we found that Fermi detected sources have systematically higher values than non-detected sources (WRS p-value ∼10⁻¹³, median δ_var ≈ 14 for detected and δ_var ≈ 8 for non-detected sources). This would suggest that the γ-ray emission could be, in part, resulting from the enhanced relativistic effects in these sources.

Zoom In Zoom Out Reset image size

There are 31 sources (10 FSRQs, five BL Lacs, five radio galaxies, 11 unclassified sources) that show δ_var < 1, seven of which have been detected by Fermi (four BL Lacs, two radio galaxy, one unclassified source). These sources are either misaligned AGNs with jets pointing away from our line of sight, and thus their emission is deboosted, or are mildly beamed and have not shown any radio-bright flaring events with T_var > 10¹¹ K during the OVRO monitoring period (2008–2017).

5.1. Lorentz Factors and Viewing Angles

Using the apparent velocity of the resolved jet components, we can estimate both Γ and θ as

$$\begin{eqnarray}&&{{\rm{\Gamma }}}_{\mathrm{var}}=\displaystyle \frac{{\beta }_{\mathrm{app}}^{2}+{\delta }_{\mathrm{var}}^{2}+1}{2{\delta }_{\mathrm{var}}},\end{eqnarray} \tag{ 3 }$$

$$\begin{eqnarray}&&{\theta }_{\mathrm{var}}=\arctan \left(\displaystyle \frac{2{\beta }_{\mathrm{app}}}{{\beta }_{\mathrm{app}}^{2}+{\delta }_{\mathrm{var}}^{2}-1}\right),\end{eqnarray} \tag{ 4 }$$

where β_app is the apparent velocity. For β_app, we use data from the MOJAVE survey (Lister & Homan 2005; Lister et al. 2016). For our calculations, we use the maximum observed apparent velocity in each jet. There are 238 sources with an available estimate, 160 of which have been detected by Fermi. Figures 6 and 7 show the Lorentz factor and viewing angle distributions for the different classes (top panels) and Fermi detected and non-detected sources (bottom panels). From the sources with an apparent velocity measurement, FSRQs have on average faster jets than any other source class (WRS p-value <0.006, median Γ_var ≈ 15.6 compared to about 10 for BL Lacs, 6.3 for radio galaxies, and 7.2 for unclassified sources). Similarly, Fermi detected sources have faster jets (∼10⁻⁶, median Γ_var ≈ 17) on average than Fermi non-detected sources (median Γ_var ≈ 9). There are four sources (namely J0108+0135, J0948+4039, J1613+3412, and C1724+4004) with a high Lorentz factor (Γ_var > 100). The derived Doppler factors for these sources are <5, yet the measured β_app,max are >20. Using the median β_app instead of β_app,max brings the Γ_var estimates to much lower values (<50). However, in all cases, the components that yielded β_app,max in each source were ejected prior to the beginning of the observations considered here. It is then possible that a major flaring event not considered in this work is associated with these components, and hence for these sources we are underestimating their Doppler factors. BL Lacs and FSRQs have similar viewing angles (WRS p-value 0.13, median θ_var ≈ 4 for BL Lacs and median θ_var ≈ 5 for FSRQs) while those for radio galaxies and unclassified sources are on average larger (p-value <0.04, median of 12°, and 9° respectively). Fermi non-detected sources also have on average larger viewing angles (p-value ∼ 10⁻⁹, median of 72) compared to Fermi detected sources (median of 27).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

A similar comparison of the beaming properties for a smaller sample (62 sources in total) in Savolainen et al. (2010) found that there is a lack of small viewing angles in the comoving frame of the jet for Fermi detected sources. The comoving-frame viewing angle (θ_src) is defined as

$$\begin{eqnarray}&&{\theta }_{\mathrm{src}}=\arccos \left(\displaystyle \frac{\cos {\theta }_{\mathrm{var}}-\beta }{1-\beta \cos {\theta }_{\mathrm{var}}}\right).\end{eqnarray} \tag{ 5 }$$

Based on the fact that Fermi detected sources have on average higher apparent velocities (Lister et al. 2009) and VLBI core brightness temperatures (Kovalev et al. 2009), the authors concluded that if the lack of small θ_src persists in a larger sample, it would suggest either anisotropy of the rest-frame γ-ray emission or a dependence of that emission on the Lorentz factor. Our sample allows us to test whether such lack of small θ_src exists. Figure 8 shows the θ_src distributions for Fermi detected and non-detected sources. With the significantly larger sample considered in this work, we do not find any lack of small θ_src for Fermi detected sources. Thus, our results disfavor scenarios involving any dependence on the Lorentz factor or anisotropic rest-frame γ-ray emission.

Zoom In Zoom Out Reset image size

Table 2 lists the values for T_var, δ_var, Γ_var, and θ_var and their uncertainties. For all sources (with or without a β_app,max estimate), we use the β_app,max distribution to bracket the possible range of Γ_var and θ_var estimates for that source given the estimated δ_var.

Table 2.  Variability Brightness Temperatures and Beaming Properties for the Sources in Our Sample

Name	Class	z	βapp,max	${\sigma }_{{\beta }_{\mathrm{app},\max }}$	Tvar	-${\sigma }_{{\text{}}{T}_{\mathrm{var}}}$	${\sigma }_{{\text{}}{T}_{\mathrm{var}}}$	δvar	$-{\sigma }_{{\delta }_{\mathrm{var}}}$	${\sigma }_{{\delta }_{\mathrm{var}}}$	Γvar	Γmin	Γmax	θvar	θmin	θmax
J0001−1551	F	2.044	⋯	⋯	11.15	−1.12	0.64	2.51	−1.44	1.69	⋯	1.45	>100	⋯	0.91	23.47
J0001+1914	F	3.100	⋯	⋯	10.29	−0.25	3.64	1.82	−0.5	30.76	⋯	1.19	>100	⋯	2.07	33.24
J0004+2019	B	0.677	⋯	⋯	11.14	−0.52	3.83	1.37	−0.48	24.59	⋯	1.05	>100	⋯	2.88	47.06
J0004+4615	F	1.810	⋯	⋯	12.78	−0.63	0.2	7.75	−2.75	1.7	⋯	3.94	>100	⋯	0.08	7.42
J0005+3820	F	0.229	⋯	⋯	13.26	−0.31	0.63	5.23	−1.22	2.72	⋯	2.71	>100	⋯	0.18	11.02
J0006−0623	B	0.347	7.31	0.33	13.55	−0.5	0.77	6.96	−2.27	5.63	7.39	3.55	>100	8.25	0.1	8.26
J0010+1058	⋯	0.089	1.58	0.29	12.48	−0.08	1.38	2.51	−0.3	4.24	1.95	1.45	>100	22.1	0.91	23.51
J0010+1724	F	1.601	⋯	⋯	14.38	−0.48	0.45	24.81	−7.46	12.21	⋯	12.43	44.18	⋯	0.01	2.31
J0010+2047	F	0.600	⋯	⋯	13.18	−1.46	0.45	6.02	−4.05	3.12	⋯	3.09	>100	⋯	0.14	9.56
J0011+0057	F	1.492	⋯	⋯	13.34	−1.37	1.92	11.03	−7.63	33.64	⋯	5.56	76.98	⋯	0.04	5.2

Note. Names are as listed in the OVRO website. The values of T_var and its uncertainties are given in log₁₀. −${\sigma }_{{\text{}}{T}_{\mathrm{var}}}$, ${\sigma }_{{\text{}}{T}_{\mathrm{var}}}$ and −${\sigma }_{{\delta }_{\mathrm{var}}}$, ${\sigma }_{{\delta }_{\mathrm{var}}}$ are the asymmetric uncertainties on T_var and δ_var, respectively. [Γ_min, Γ_max] and [θ_min, θ_max] are the possible minimum and maximum values of each source for a given δ_var by marginalizing over the β_app,max distribution. The redshift estimates are taken from Richards et al. (2014), SIMBAD (Wenger et al. 2000), NASA/IPAC Extragalactic Database (NED; operated by the Jet Propulsion Laboratory, California Institute of Technology, under contract with the National Aeronautics and Space Administration), and the MOJAVE database (Lister et al. 2018). The table lists only the first 10 sources. It is published in its entirety in machine-readable format. A portion is shown here for guidance regarding its form and content.

Only a portion of this table is shown here to demonstrate its form and content. A machine-readable version of the full table is available.

Download table as:  DataTypeset image

5.2. Sources without Redshift

5.3. Comparison with Other Doppler Factor Estimation Methods

There are several methods in the literature for estimating the Doppler factor in blazar jets, some of which are mentioned in Section 1. Although a broader comparison study among the different methods similar to Liodakis & Pavlidou (2015b) and Liodakis et al. (2017c) could be beneficial, we focus on recent results from the radio regime and SED modeling. The most recent attempts in estimating the variability Doppler factor for a large number of sources are Hovatta et al. (2009) and Liodakis et al. (2017d; hereafter H09 and L17, respectively). In H09, the authors used data from the Metsähovi monitoring program at 22 and 37 GHz (Teraesranta et al. 1998) and estimated the variability brightness temperature for 87 sources by modeling the light curves using the same exponential rise and exponential decay model as this work. L17 used multiwavelength radio data (2.64–142.33 GHz) from the F-GAMMA program (Fuhrmann et al. 2016) to decompose the light curves using non-parametric models for the flare profiles tailored to the individual light curves of 58 sources. The cadence of H09 is weekly, while the cadence for L17 is every two weeks to monthly. All of the sources in H09 and L17 are also in our sample. Our results appear to be consistent with both studies with roughly 50% of the estimates consistent within 1σ. However, both H09 and L17 have assumed that T_eq =5 × 10¹⁰. Once we account for the different T_eq, the estimates derived in this work become larger by a factor of ≈1.85. The number of sources with consistent estimates drops to roughly 20%, and there is now a significant difference in the Doppler factor distributions with estimates of this work being systematically larger (WRS p-value < 0.0002 for both samples). The higher Doppler factors from this work are most likely due to OVRO's faster cadence. However, cadence may not be solely responsible for the differences between the estimates. The data set used in H09 includes observations up to roughly 2006. While there is overlap between the observing periods of L17 and this work, the estimates in both H09 and L17 originate from a variety of frequencies, which may probe regions not co-spatial with the one probed at 15 GHz due to synchrotron self-absorption. It is then possible for the differences in the estimates to also be attributed to either significant flaring events that have occurred outside the periods considered in H09 and L17 or that their reported estimates simply corresponds to different regions of the jet. Additionally, results from the F-GAMMA survey would suggest a decreasing trend of the brightness temperature with frequency T_var ∝ ν^−1.2 (Fuhrmann et al. 2016). From Equation (2), the Doppler factor should then decrease as δ_var ∝ ν^−0.4 with increasing frequency. Such a trend could imply that the jets are accelerating from the higher to the lower radio frequencies, which could explain some of the discrepancies. About 75% of blazars in the MOJAVE survey have indeed shown at least one accelerating jet feature at 15 GHz (Homan et al. 2015). However, the fact that a significant fraction of the estimates in L17 are estimated at a lower frequency than 15 GHz would suggest that this scenario is unlikely to explain the discrepancies between the estimates from L17 and this work.

A more interesting comparison would be with the estimates in Jorstad et al. (2017, hereafter J17). The method uses the variability timescales of individual jet components which are related to the Doppler factor through the observed angular size of the components derived from VLBI observations at 43 GHz (Jorstad et al. 2005). Although the method is also limited by the cadence of the observations, it has the advantage of being independent of the assumption of equipartition. Thus, agreement in the estimates of the two methods (J17 and this work) provides strong constrains for the Doppler factors of the jets. All of the sources (36) in J17 are included in the present sample. The estimates for 11/36 sources are consistent within 1σ. Differences in the estimates between the two samples are most likely attributed to the different assumptions used in each method or to reasons described above; however, no systematic difference is detected between the Doppler factor distributions according to the WRS test (p-value 0.56). The names and δ estimates of the sources in agreement between the two samples are given in Table 4.

Table 4.  List of Sources with Doppler Factors Consistent with J17 within 1σ

OVRO name	J17 name	δvar	δJ17
J0238+1636	0235+164	${43.53}_{-11.49}^{+19.79}$	52.8 ± 8.4
J0339−0146	0336−019	${23.09}_{-6.08}^{+18.9}$	15.7 ± 4.9
0415+379	0415+379	${1.99}_{-0.43}^{+1.55}$	2.0 ± 0.5
J0433+0521	0430+052	${4.16}_{-1.09}^{+1.42}$	4.5 ± 2.0
J0830+2410	0827+243	${31.97}_{-4.2}^{+5.83}$	22.8 ± 8.5
C1224+2122	1222+216	${5.32}_{-1.76}^{+7.84}$	7.4 ± 2.1
J1229+0203	1226+023	${3.78}_{-0.55}^{+1.1}$	4.3 ± 1.3
J1310+3220	1308+326	${26.35}_{-16.63}^{+13.39}$	20.9 ± 1.2
PKS 1510−089	1510−089	${32.14}_{-7.97}^{+8.07}$	35.3 ± 4.6
J1751+0939	1749+096	${17.62}_{-3.16}^{+10.1}$	17.7 ± 7.7
J2253+1608	2251+158	${26.61}_{-2.97}^{+6.28}$	24.4 ± 3.7

Note. Names are as given in the OVRO website and J17.

Download table as:  ASCIITypeset image

SED modeling has also been used to constrain the Γ and θ in blazar jets in part due to the fact that different γ-ray emission mechanisms are affected differently by the relativistic effects (e.g., Dermer 1995). Recent SED modeling of a large number of sources found that the distribution of the derived Lorentz factors (a frequent assumption in SED modeling is that δ = Γ) is narrow, peaking at δ = Γ ∼ 13 ± 1.4 (Ghisellini et al. 2014). Similar results were found in Chen (2018) considering a larger sample (δ = Γ ∼ 14). It is usually assumed in SED modeling that the γ-ray emission is produced closer to the supermassive black hole than the radio core of the jet where most of the radio emission originates. It is then interesting that we find similar results for δ_var for the Fermi detected sources (median δ_var ≈ 14). The derived Γ_var in this work appears to be on average larger (median Γ_var ≈ 17); however, the distribution is wide enough to prevent us from investigating any potential discrepancy. Agreement between the two methods could suggest that there is no significant change in the relativistic effects between the radio and γ-ray emission regions, which has implications for the different jet acceleration scenarios as well as the possible location of the γ-ray production site. However, given the complexity of the SED models and the covariance between the different parameters involved in these models, any agreement could be artificial. A dedicated study of the sources studied in this work could allow us to probe possible differences in the relativistic effects between radio and γ-rays.

By modeling with a superposition of flares on top of a stochastic background the radio light curves from the OVRO 40 m telescope's blazar monitoring program, we were able to estimate the variability brightness temperatures and Doppler factors for 1029 sources, the largest set of estimates available to date. OVRO's high cadence allowed us to resolve even the fastest flares and set the strongest constraints on the highest T_var in each source. The present analysis is, however, limited by the time span of observations. It is possible for significant flaring events to have occurred outside the observing time span considered here as the variability timescales in blazars are typically long (Hovatta et al. 2007). Thus, for all intents and purposes, the results from this work should be treated as lower limits. The fact that roughly 12% of our sources have Doppler factors as high as δ_var > 30 would suggest that at least for a fraction of our sample, we were able to estimate the "true" δ_var of the jet. It would be productive to repeat such analysis with light curves observed during time intervals (with similar cadence) different from the one considered here in order to examine whether the highest T_var has indeed been estimated for each source.

The majority of flares with rise times <14 days have lower amplitudes than the median flare amplitude of the entire light curve, suggesting that OVRO's three-day sampling allowed us to resolve all of the major events within the time span of the observations used in this work. Although low amplitude variability is still possible on shorter timescales, it can be adequately described by a stochastic background process. Intra-day variability has been found in a handful of the brightest radio sources so it could be interesting to observe sources at an even faster cadence than three days; however, such fast variations are often attributed to interstellar scintillation and not to intrinsic processes.

Our results show a significant difference between the T_var distributions of Fermi detected sources and non-detected sources. A similar result was obtained by Kovalev et al. (2009) when comparing the median VLBI brightness temperature of Fermi detected sources and non-detected sources from the MOJAVE survey. A more in-depth comparison of their flaring properties showed significant differences between Fermi detected and non-detected sources, with the former showing on average faster flares with higher amplitudes. A comparison of the radio flux-density distributions of Fermi detected and non-detected sources in Liodakis et al. (2017b) also showed that γ-ray loud sources are systematically more variable and have higher flaring ratios (ratio of the flaring to quiescent mean flux densities) than γ-ray quiet sources. Our findings extend that result, showing that both the variability and flaring properties in radio are connected to the γ-ray activity. This would suggest that the underlying mechanism in the jet that would cause higher and more energetic flares in radio is, at least in part, also responsible for the γ-ray emission.

Using population models and the T_var estimates, we were able to effectively constrain the equipartition brightness temperature to $\langle {T}_{\mathrm{eq}}\rangle =2.78\times {10}^{11}\,{\rm{K}}$ (±26%). Previous attempts on constraining the limiting intrinsic brightness temperature had either estimated T_eq to be between 10¹⁰ and 10¹¹ K (Readhead 1994; Lähteenmäki et al. 1999; Cohen et al. 2003) or most recently constrained it to T_eq > 2 × 10¹¹ K (Homan et al. 2006). The very high cadence of the OVRO program allowed us to resolve even the fastest events pushing the limit of the highest estimated T_var and hence provide the strongest constraints on T_eq. Interestingly, our results are consistent with the theoretical expectations for the limiting brightness temperature for incoherent synchrotron sources due to magnetization effects (∼3 × 10¹¹ K; Singal 1986). Although our results are model dependent, the fact that they are in agreement with both observational (Homan et al. 2006) and theoretical (Singal 1986) expectations for blazar jets gives us confidence in our analysis.

Based on the results of the T_eq optimization, we estimated δ_var and its uncertainty for the sources in our sample, significantly increasing the number of available δ_var estimates in the literature. As expected, blazars are highly beamed sources with larger δ_var on average than either radio galaxies or unclassified sources (median of ≈11 for blazars compared to a median of ≈5 for radio galaxies and unclassified sources). Surprisingly, we do not detect any significant difference between the BL Lacs and FSRQ populations contrary to the current consensus suggesting that FSRQs are more beamed than BL Lacs (H09, L17). This, of course, could be due to the selection of T_eq to be the same for all populations. If BL Lacs were allowed to have the very high T_eq found in the above analysis (although not compatible with our current understanding of jet processes), it would lower their δ_var estimates by a factor of ∼2.5. A more plausible explanation for this discrepancy could be that previous monitoring programs (with a slower cadence than OVRO's) were not able to resolve the BL Lac flares evolving on faster timescales (see discussion in Section 3), but were still able to detect the more slowly evolving FSRQ flares. In such a case, it is only natural that FSRQs would show higher T_var and hence larger δ_var than BL Lacs. Contrary to previous monitoring programs, OVRO's fast cadence allowed us to detect all prominent flares in both BL Lacs and FSRQs.

Although there is no significant difference between the viewing angle distributions of the blazar classes (not surprising if the sources are uniformly distributed and randomly oriented), FSRQs host faster jets than BL Lacs. This could help explain differences between the two populations, or at least, differences between FSRQs and radio-bright BL Lacs. As expected, Fermi detected sources have on average faster jets pointed at smaller angles toward our line of sight than Fermi non-detected sources. Then the relativistic effects could also be partly responsible for the detected (or not) γ-ray emission in addition to radio variability and flaring properties (see also Lister et al. 2009). It should be noted that while the flux-density variations used to estimate the brightness temperature and Doppler factors originate predominantly in the radio core of the jet and are believed to be related to ejections of new radio components (e.g., Savolainen et al. 2002), apparent velocities are measured downstream from the core. Given that the observations used in this work and the apparent velocity measurements from the MOJAVE program are taken at the same radio frequency (15 GHz, and thus probe the same region for a given source), we do not expect significant changes in the velocity of the jet over short distances. However, since both accelerating and decelerating jet components have been measured at 15 GHz (e.g., Homan et al. 2015), our results for the Lorentz factors and viewing angles in sources that show large velocity gradients should be treated as limits. Additionally, to derive the Γ_var, θ_var estimates, we have used the maximum observed apparent velocity in each jet. Although it has been shown that the radio components of individual jets are ejected at similar velocities (Lister et al. 2013), using a different measure for β_app (e.g., mean, median) could result in differences in the Γ_var and θ_var estimates.

A comparison with previous attempts (H09, L17) in estimating δ_var showed that (after accounting for the different assumed T_eq) the estimates from this work are systematically higher with only 20% of the common estimates to be consistent within 1σ. This is not surprising given the faster cadence of the OVRO survey. On the other hand, comparing our estimates with those from J17 derived using a different approach independent of equipartition showed that 30.5% of the sources in the J17 sample are consistent with the estimates from this work. This agreement allows us to place strong constraints on the δ estimates for these sources.

We thank the anonymous referee for comments and suggestions that helped improve this work. 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. This research has made use of data from the MOJAVE database, which is maintained by the MOJAVE team (Lister et al. 2018). This research has made use of the SIMBAD database, operated at CDS, Strasbourg, France (Wenger et al. 2000). This research has made use of the NASA/IPAC Extragalactic Database (NED), which is operated by the Jet Propulsion Laboratory, California Institute of Technology, under contract with the National Aeronautics and Space Administration.

Facility: OVRO. -

Software: Magnetron (Huppenkothen et al. 2015), DNest4 (Brewer & Foreman-Mackey 2016), Numpy (Van Der Walt et al. 2011), Scipy (Jones et al. 2001).