Scale Invariant Jets: From Blazars to Microquasars

Blazars constitute unique laboratories to study extreme astrophysics, from relativistic magnetohydrodynamics and shocks to particle acceleration, ultra-high energy cosmic rays, and neutrino production. They are divided into two subclasses, Flat Spectrum Radio Quasars (FSRQs) and BL Lac objects (BL Lacs), and are most famous for their extreme variability, apparent superluminal motion and γ-ray loudness since they comprise the largest detected population of the Fermi γ-ray observatory (Acero et al. 2015). Understanding the relativistic highly collimated plasma outflows of blazar jets has proven extremely difficult due to the relativistic effects dominating their emission from radio to γ-rays (Blandford & Königl 1979). These relativistic effects are quantified by the Doppler factor $\delta ={[\gamma (1-\beta \cos \theta )]}^{-1}$, where γ is the Lorentz factor ($\gamma \,={(\sqrt{1-{\beta }^{2}})}^{-1}$), β is the velocity of the jet in units of speed of light, and θ is the angle between their jet axes and the observer's line of sight. Even a small spread in the values of γ and θ among blazars results in a large spread in observed properties, severely complicating the search for empirical correlations that can confirm or constrain jet models. If δ could be confidently estimated, however, these relativistic effects could be corrected for and blazar jets could be studied in their emission rest frame.

Several methods have been proposed for estimating blazar Doppler factors, but they frequently yield discrepant results (Liodakis et al. 2017b). Given the numerous assumptions entering each method, it is often challenging to identify the most accurate estimate for any given blazar. Two recent breakthroughs have, however, made such a task tractable for the first time. First, through population modeling of unbiased blazar samples (Liodakis & Pavlidou 2015a), Doppler factor estimates based on variability studies and the assumption of equipartition between synchrotron emitting particles and magnetic field (Readhead 1994; Lähteenmäki & Valtaoja 1999; Hovatta et al. 2009) were shown to be the most accurate (Liodakis & Pavlidou 2015b). Second, multi-frequency F-GAMMA radio data have recently enabled the calculation of the highest-ever accuracy variability Doppler factors for 58 well-studied blazars (Liodakis et al. 2017a).

Correlations between the BH mass (${M}_{\mathrm{BH}}$) and monochromatic radio flux density, or the monochromatic radio flux density and the X-ray flux density or even all the above combined, have long been established (e.g., Merloni et al. 2003). The latter suggests the existence of a plane, termed the "fundamental plane of black hole activity," which extends from X-ray binaries to active galaxies. These results support the hypothesis of scale invariance, which implies that the jet formation processes are independent of the black hole mass of the system. Such a hypothesis has been predicted by theoretical models (Heinz & Sunyaev 2003).

Merloni et al. (2003) as well as similar attempts to establish a relation connecting BH-powered jets of different ${M}_{\mathrm{BH}}$ (Nagar et al. 2002; Falcke et al. 2004; Körding et al. 2006; Plotkin et al. 2012; Saikia et al. 2015) have either explicitly avoided blazars focusing on low-luminosity active galactic nuclei (LLAGNs) or have included a handful of blazars. In the latter case, they have either ignored the relativistic effects or specifically chosen their sample to include only low-beamed sources. In cases where there was some treatment of the relativistic effects (Falcke et al. 2004; Körding et al. 2006; Plotkin et al. 2012), the common practice was to use a single value of δ, when in reality Doppler factors are estimated to range between 1 and 45 based on individual source studies (Lähteenmäki & Valtaoja 1999; Fan et al. 2009; Hovatta et al. 2009; Liodakis et al. 2017a) and between 1 and 60 based on population studies (Liodakis & Pavlidou 2015a).

However, studies focused on blazars, fully and accurately accounting for their relativistic effects cannot be circumvented, especially if the aim is to study the physics of jets: it is only in such highly beamed sources that we are certain the observed spectrum is dominated by the jet emission. For this reason, the extension of such scalings to blazars should have a strong impact on unification models of radio-loud active galactic nuclei and our knowledge of BH-powered jets, increasing manyfold our ability to constrain jet models.

There have been attempts to create similar scaling relations in blazars either using luminosity–luminosity correlations (Nemmen et al. 2012; Fan et al. 2017) or different jet quantities with the properties of the central engine (Wang et al. 2004; Hovatta et al. 2010; Bower et al. 2015). However, the relativistic effects hamper any attempt to establish a strong correlation between the rest-frame emission of the jet and the ${M}_{\mathrm{BH}}$ in beamed sources.

Even if we account for the relativistic effects properly, blazars are known to show extreme variability across all frequencies. Therefore, the use of single-epoch measurements of flux densities as a proxy of the source luminosity is highly problematic. In addition, many different mechanisms (e.g., hot corona, synchrotron radiation, and synchrotron self-Compton) contribute to the blazar X-ray flux. Disentangling the different contributions has been so far extremely uncertain at best, making the use of the black hole fundamental plane or similar relations unfeasible.

In this work, in order to overcome all such limitations and probe the physics of blazars, we propose a new method to explore the connection between jet power and supermassive black holes in beamed sources by using the rest-frame broadband radio luminosity (${{ \mathcal L }}_{\mathrm{int}}^{\mathrm{Br}-{\rm{B}}}$) and the ${M}_{\mathrm{BH}}$.

This work is organized as follows. In Section 2, we present the data used in this work and the manner of their analysis. In Section 3, we present the correlation analysis and the best-fit relation. In Section 4, we discuss our findings, and in Section 5 we summarize our results and conclusions.

Throughout this work, the adopted cosmology is ${H}_{0}=71\,\mathrm{km}\,{{\rm{s}}}^{-1}\,{\mathrm{Mpc}}^{-1}$, ${{\rm{\Omega }}}_{m}=0.27$, and ${{\rm{\Omega }}}_{{\rm{\Lambda }}}=1-{{\rm{\Omega }}}_{m}$ (Komatsu et al. 2009).

Our sample is a sub-set of the sources monitored by the F-GAMMA program,⁸ including all sources with both an available Doppler factor and ${M}_{\mathrm{BH}}$ estimates. The F-GAMMA program monitored a total of about 100 blazars detected by Fermi for eight years with roughly monthly cadence at 10 frequencies from 2.64 to 142.33 GHz (Fuhrmann et al. 2016). The Doppler factors are taken from Liodakis et al. (2017a). In Liodakis et al. (2017a), a novel approach is used to model and track the evolution of the flares through multiple frequencies. That method has provided the most accurate variability Doppler factor estimates to date with an average error of 16%. We use only estimates with quality indicator "confident" or "very confident" to ensure robust results. The level of confidence is defined by the number of available frequencies and flares used in the estimation of the Doppler factor (see Liodakis et al. 2017a). The number of sources in our sample is 26 (20 FSRQs, 4 BL Lacs, and 2 radio galaxies). In order to extend our sample toward lower black hole masses, we also included three γ-ray loud Narrow Line Seyfert 1 galaxies (NLS1s; Angelakis et al. 2015) for which there was a Doppler factor computed with the same method as in Liodakis et al. (2017a). In total, our sample consists of 29 sources with ${M}_{\mathrm{BH}}$ spanning from $\sim {10}^{6.5}$ to ${10}^{9.5}\,{M}_{\odot }$.

2.1. Intrinsic Broadband Radio Luminosity

We calculated the intrinsic broadband radio luminosity (${{ \mathcal L }}_{\mathrm{int}}^{\mathrm{Br}-{\rm{B}}}$) of our sources using data from the F-GAMMA program. For each frequency, we use the maximum likelihood approach described in Richards et al. (2011) to calculate the maximum likelihood mean flux density, taking into account errors in measurements, uneven sampling, and source variability. From that we construct the "mean" spectrum, which we convert to the intrinsic luminosity using the Doppler factor estimates from Liodakis et al. (2017a). Accounting for the relativistic effects in estimating the rest-frame broadband radio luminosity in blazars is crucial in order to identify the correct correlation (see Section 4).

We convert the maximum likelihood mean flux density for each frequency to the intrinsic luminosity using the following equation:

$$\begin{eqnarray}&&{S}_{\nu }=\displaystyle \frac{{L}_{\nu }{\delta }^{p}}{4\pi {d}_{L}^{2}}{(1+z)}^{1+s},\end{eqnarray} \tag{ 1 }$$

where S_ν is the flux density at a given frequency ν, L_ν is the intrinsic luminosity at that frequency, δ is the Doppler factor, d_L is the luminosity distance, z is the redshift, s is the spectral index, which we have defined as $S\propto {\nu }^{s}$, and p is equal to $p=2-s$ for the continuous, and $p=3-s$ for the discrete jet case.

By integrating over frequency, we calculate the ${{ \mathcal L }}_{\mathrm{int}}^{\mathrm{Br}-{\rm{B}}}$ for each source, assuming either a continuous jet (${{ \mathcal L }}_{\mathrm{int},{\rm{c}}}^{\mathrm{Br}-{\rm{B}}}$) or a discrete jet (${{ \mathcal L }}_{\mathrm{int},{\rm{d}}}^{\mathrm{Br}-{\rm{B}}}$). We estimate the uncertainty of the ${{ \mathcal L }}_{\mathrm{int}}^{\mathrm{Br}-{\rm{B}}}$ through formal error propagation taking into account errors in the Doppler factor (Liodakis et al. 2017a), flux density, spectral index, and redshift. For the error in the Doppler factor estimates for the three NLS1s, we assumed the 16% average error from Liodakis et al. (2017a). Since our redshift estimates are all spectroscopic, we have assumed a common error of $\delta z=0.01$.

Table 1.  Range of Estimates for the Different Parameters of the Microquasars

Name	Alt. Name	Distance	Γ	θ	$\mathrm{log}{{L}}_{\mathrm{int},{\rm{c}}}^{\mathrm{Br}-{\rm{B}}}$	$\mathrm{log}{M}_{\mathrm{BH}}$
V* V1343 Aql	SS 433	3–5.5	1.02–1.05	73–85	26.98–27.52	2.7–7.9
V* V1487 Aql	GRS 1915+105	8.6–13.7	1.01–5.03	60–71	26.44–28.29	12.4–15
V* V1521 Cyg	Cyg X-3	7.2–12	1.05–2.40	20–80	28.96–30.44	1–30
V* V1033 Sco	GRO J1655-40	3–3.5	1.04–4.12	70–85	25.99–27.68	5.10–7.02
V* V404 Cyg	NOVA Cyg 1989	2.25–2.53	1.2–10.0	46–73	26.16–27.86	3–14

Note. (1) SIMBAD identification name, (2) alternative source name, (3) distance in kiloparsecs, (4) lorentz factor (Γ), (5) viewing angle (θ), (6) logarithm of the intrinsic broadband radio luminosity (Watt), and (7) logarithm of the BH mass in solar masses.

Download table as:  ASCIITypeset image

2.2. BH Mass Estimation

In order to estimate the ${M}_{\mathrm{BH}}$ of our sample, we used data from Torrealba et al. (2012) and the scaling relations from Shaw et al. (2012) using the line luminosities of the MgII and Hb spectral lines. The scalings were calibrated using the virial mass estimation method. The uncertainty of the ${M}_{\mathrm{BH}}$ for each source was the result of error propagation. We complemented our sample with values (and their error) from the literature with estimates using the same method for consistency (Yuan et al. 2008; Shaw et al. 2012; Zamaninasab et al. 2014, and references therein). When there was no error estimate for a literature value (Yuan et al. 2008; Zamaninasab et al. 2014, see Table 2), we assumed as error the average uncertainty of the available estimates using the same spectral line. We avoided using methods that involve assumptions on the radio flux density, beaming, Spectral Energy Distribution (SED) fitting, etc., that could potentially create artificial correlations. All the estimates for the ${{ \mathcal L }}_{\mathrm{int}}^{\mathrm{Br}-{\rm{B}}}$ and ${M}_{\mathrm{BH}}$ used in this work are summarized in Table 2.

Table 2.  Blazar Sample

Name	ALT-Name	Class	z	${\delta }_{\mathrm{var}}$	${\sigma }_{{\delta }_{\mathrm{var}}}$	$\mathrm{log}{{L}}_{\mathrm{int},{\rm{c}}}^{\mathrm{Br}-{\rm{B}}}$	${\sigma }_{\mathrm{log}{L}}$	$\mathrm{log}{M}_{\mathrm{BH}}$	${\sigma }_{M}$
J0102+5824	0059+5808	Q	0.644	21.9	3.6	35.88	0.05	8.57a	0.61
J0136+4751	0133+476	Q	0.859	13.7	3.7	36.45	0.03	8.81	0.10
J0237+2848	0234+285	Q	1.206	12.2	4.3	36.75	0.02	9.22a	0.30
J0324+3410	1H0323+342	N	0.063	3.91	0.6	34.61	0.05	7.25c	0.13
J0418+3801	3C111	G	0.049	2.0	0.4	35.76	0.02	8.31	0.10
J0423-0120	0420-014	Q	0.916	43.9	9.2	35.81	0.21	8.84	0.07
J0433+0521	3C120	G	0.033	2.1	0.1	35.10	0.02	7.58	0.14
J0530+1331	PKS0528+134	Q	2.070	12.9	2.5	37.18	0.05	9.03b	0.25
J0654+4514	S40650+453	Q	0.928	13.8	2.6	35.85	0.04	8.17a	0.34
J0948+0022	PMN J0948+0022	N	0.583	10.12	1.6	35.76	0.05	7.5c	0.13
J1130-1449	1127-145	Q	1.184	21.9	0.0	36.07	0.07	9.30	0.10
J1159+2914	PKS1156+295	Q	0.725	12.8	0.0	36.09	0.06	8.61a	0.20
J1221+2813	QSOB1219+285	B	0.102	2.6	0.6	35.30	0.02	8.89	0.09
J1229+0203	3C273	Q	0.158	3.7	1.0	36.74	0.02	8.76	0.13
J1256-0547	3C279	Q	0.536	16.8	2.9	36.75	0.03	8.90	0.12
J1310+3220	OP+313	B	0.997	15.8	1.7	36.29	0.02	8.57a	0.18
J1504+1029	PKS1502+106	Q	1.839	17.3	2.7	36.43	0.07	8.98a	0.24
J1505+0326	PKS1502+036	N	0.408	11.32	1.8	34.52	0.03	6.6c	0.13
J1512-0905	PKS1510-089	Q	0.360	12.3	2.8	35.85	0.06	8.20b	0.13
J1635+3808	4C+38.41	Q	1.814	20.3	2.8	36.96	0.06	9.30a	0.16
J1642+3948	3C345	Q	0.593	10.4	2.9	36.46	0.02	9.03a	0.19
J1800+7828	S51803+78	B	0.680	21.2	5.0	35.69	0.03	7.92b	0.13
J1848+3219	TXS1846+322	Q	0.798	12.1	1.4	36.18	0.13	8.21a	0.53
J1849+6705	S41849+670	Q	0.657	8.1	1.4	36.55	0.02	8.81	0.07
J2202+4216	BL Lac	B	0.069	6.1	0.8	35.30	0.03	8.23b	0.13
J2229-0832	2227-088	Q	1.560	21.0	0.6	36.49	0.05	8.70a	0.16
J2232+1143	CTA102	Q	1.037	15.1	4.8	36.61	0.04	8.93	0.07
J2253+1608	3C454.3	Q	0.859	17.0	3.7	37.23	0.06	9.19	0.11
J2327+0940	PKS2325+093	Q	1.841	17.2	2.3	36.25	0.04	8.70a	0.30

Notes. (1) F-GAMMA identification name, (2) alternative source name, (3) class (Q is for FSRQs, B for BL Lacs, G for radio galaxies, N for narrow line Seyferts), (4) redshift, (5) variability Doppler factor (${\delta }_{\mathrm{var}}$), (6) uncertainty of the variability Doppler factor (${\sigma }_{{\delta }_{\mathrm{var}}}$), (7) logarithm of the intrinsic broadband radio luminosity (Watt), (8) uncertainty of the luminosity, (9) logarithm of the BH mass in solar masses, and (10) uncertainty of the BH mass.

^aShaw et al. (2012). ^bZamaninasab et al. (2014). ^cYuan et al. (2008).

Download table as:  ASCIITypeset image

Assuming that the jet is composed of a series of plasma blobs (discrete jet), we tested for a correlation between the ${{ \mathcal L }}_{\mathrm{int},{\rm{d}}}^{\mathrm{Br}-{\rm{B}}}$ and ${M}_{\mathrm{BH}}$. We used the partial correlation test (Akritas & Siebert 1996) in order to obtain the Kendall correlation coefficient (τ) and its significance taking into account the effects of redshift. The partial correlation test yielded a correlation coefficient of $\tau =0.35$ with $2\times {10}^{-3}$ probability of uncorrelated samples, indicating a moderate correlation. If we assume the jet is a continuous stream of plasma (continuous jet, ${{ \mathcal L }}_{\mathrm{int},{\rm{c}}}^{\mathrm{Br}-{\rm{B}}}$), the correlation becomes significantly stronger ($\tau =0.51$ with 5 × 10⁻⁶ probability of no correlation). The correlation is stronger in this case, which is to be expected since the continuous jet is, most likely, the most accurate representation of the jet structure, as has been repeatedly supported by interferometric radio observations (VLBI, e.g., Zensus 1997). For this reason, for the remainder of this work, we focus on the continuous jet case. Examining different subsamples there was no case were the p-value was $\gt {10}^{-3}$, supporting the robustness of the correlation. For FSRQs alone, we find $\tau =0.44$ with p-value 6 × 10⁻⁴. Including the BL Lacs, the test yielded τ = 0.43 with p-value 7 × 10⁻⁴, and including the radio galaxies τ = 0.45 with probability 2 × 10⁻⁴. The significant correlation found in all cases strongly suggests that this is a real trend and not an artifact of the large range of ${M}_{\mathrm{BH}}$.

Having established that a significant correlation exists, we performed a fit between ${{ \mathcal L }}_{\mathrm{int},{\rm{c}}}^{\mathrm{Br}-{\rm{B}}}$ and ${M}_{\mathrm{BH}}$ in log-space. Assuming a linear model of the form $\mathrm{log}{{ \mathcal L }}_{\mathrm{int},{\rm{c}}}^{\mathrm{Br}-{\rm{B}}}=A\times \mathrm{log}\left(\tfrac{M}{{10}^{8}\,{M}_{\odot }}\right)+B$. The fit was performed using the BCES bisector methods described in Akritas & Bershady (1996), which takes into account errors in both axis as well as intrinsic scatter. The best-fit line between the ${{ \mathcal L }}_{\mathrm{int},{\rm{c}}}^{\mathrm{Br}-{\rm{B}}}$ and the ${M}_{\mathrm{BH}}$ would then be

$$\begin{eqnarray*}&&\mathrm{log}\left(\displaystyle \frac{{{ \mathcal L }}_{\mathrm{int},{\rm{c}}}^{\mathrm{Br}-{\rm{B}}}}{{Watt}}\right)\\ &&\,\,=\,(1.12\pm 0.13)\times \mathrm{log}\left(\displaystyle \frac{M}{{10}^{8}\,{M}_{\odot }}\right)+(35.5\pm 0.1).\end{eqnarray*}$$

The $\mathrm{log}{{ \mathcal L }}_{\mathrm{int},{\rm{c}}}^{\mathrm{Br}-{\rm{B}}}$ versus $\mathrm{log}{M}_{\mathrm{BH}}$ plot is shown in Figure 1 together with the best-fit line. The best-fit results do not depend on the adopted linear best-fit method. For example, if we consider only intrinsic scatter using the bisector method of Isobe et al. (1990), the best-fit of the slope and intercept become $A=1.05\pm 0.01$ and $B=35.55\pm 0.04$. In fact, for any of the methods described in Akritas & Bershady (1996) and Isobe et al. (1990), the results remain consistent, within the uncertainties.

Zoom In Zoom Out Reset image size

There is some scatter around the best-fit relation, which if intrinsic has been taken into account during the fit. To examine whether this scatter is induced by the errors in the measurements or it is intrinsic, we compare the average distance of each measurement from the best-fit relation to the average distance due to the uncertainty. Using the best-fit line, for each measurement of ${M}_{\mathrm{BH}}$, we estimated the predicted intrinsic radio broadband luminosity and subtracted it from the observed one to determine the distance in the y-axis. We perform the reverse for the x-axis. We calculate the vertical distance for each observation using $z=\sqrt{{x}^{2}+{y}^{2}}$ and the scatter is estimated as $\mathrm{Sc}=\sqrt{({\sum }_{i=1}^{N}{z}_{i}^{2})/N}$. For the expected distance due to error, we use the average error on the y- and x-axes and calculate the vertical distance and scatter in the same manner. We found that the scatter of our sample is $\mathrm{Sc}=0.56$, whereas the scatter due to error is ${\mathrm{Sc}}_{\mathrm{error}}=0.19$, suggesting that the scatter around the best-fit scaling is intrinsic.

Previous studies have opted to use the monochromatic radio luminosity of the jet (Nagar et al. 2002; Merloni et al. 2003). We can estimate the mean radio luminosity at a fixed frequency and account for the blazar variability using maximum likelihood. However, observer's frame frequencies are significantly different in the emission rest frame due to the large redshift span ($z=[0,2]$) of our sample. Since different rest-frame frequencies probe different emission regimes (optically thick or optically thin) in different blazars as well as different regions, the use of single-frequency measurements should affect the scatter of the correlation (see also the discussion in Falcke et al. 2004). Figure 2 shows the intrinsic monochromatic luminosity at 4.8 GHz versus ${M}_{\mathrm{BH}}$ for our sample. There is a strong correlation ($\tau =0.44$, p-value $2\times {10}^{-4}$) and a slope of $A=0.83\pm 0.17$ and $B=24.8\pm 0.15$. Although we find consistent results between monochromatic and broadband luminosities, there is a larger scatter ($\mathrm{Sc}=0.70$) around the best-fit line in the former case. We therefore conclude that using the broadband luminosity provides stronger constrains for the ${{ \mathcal L }}_{\mathrm{int}}\mbox{--}{M}_{\mathrm{BH}}$ relation.

Zoom In Zoom Out Reset image size

If jets are indeed scale invariant, then stellar-mass BH systems (i.e., microquasars) will have intrinsic broadband radio luminosities of $\sim {10}^{27}\mbox{--}{10}^{28}$ Watt according to the best-fit relation derived above. To test this prediction, we collected archival data from the literature for five well-studied microquasars with available contemporaneous multi-wavelength radio observations. Given the uncertainty in the measurements of different parameters for the microquasars (viewing angle, jet velocity, distance, ${M}_{\mathrm{BH}}$), with each source having multiple estimates in the literature, we used a Monte-Carlo approach to calculate the mean and spread (minimum to maximum) of the ${{ \mathcal L }}_{\mathrm{int},{\rm{c}}}^{\mathrm{Br}-{\rm{B}}}$ and ${M}_{\mathrm{BH}}$ for each source (see the Appendix).

Figure 3 shows the position of the microquasars with respect to the best-fit line derived from the supermassive ${M}_{\mathrm{BH}}$ sample. The values represent the mean and the errorbars represent the spread given the different estimates for each source. All microquasars are consistent with the best-fit ${{ \mathcal L }}_{\mathrm{int},{\rm{c}}}^{\mathrm{Br}-{\rm{B}}}\mbox{--}{M}_{\mathrm{BH}}$ relation for blazars, within the 3σ confidence area of the fit, straddling the best-fit line. In fact, four out of the five sources are within 1σ. While uncertainties of both fit and measurements are quite significant at this mass range, our results suggest that the scaling derived from the blazar sample is in fact a universal scaling extending over at least ∼9 orders of magnitude both in ${{ \mathcal L }}_{\mathrm{int},{\rm{c}}}^{\mathrm{Br}-{\rm{B}}}$ and ${M}_{\mathrm{BH}}$.

Zoom In Zoom Out Reset image size

Taking into account the relativistic effects is essential in order to identify the correct scaling. We investigated the ${{ \mathcal L }}^{\mathrm{Br}-{\rm{B}}}\mbox{--}{M}_{\mathrm{BH}}$ relation obtained after only correcting for redshift. A correlation is detected; however, the sources appear to occupy unrelated regions of the ${{ \mathcal L }}^{\mathrm{Br}-{\rm{B}}}\mbox{--}{M}_{\mathrm{BH}}$ plot (Figure 4 upper panel). The slope is much steeper (A = 2.45 ± 0.22) with larger scatter. Most importantly, the scaling from the supermassive BHs underpredicts the broadband radio luminosity of microquasars by roughly 10 orders of magnitude. Accounting for beaming not only brings all classes with supermassive BHs onto one single line, but it also accurately predicts the position of their stellar-mass BH counterparts.

Zoom In Zoom Out Reset image size

A nonlinear relation between the flux of the jet (and hence the luminosity) and the ${M}_{\mathrm{BH}}$ (${S}_{\nu }\propto {M}_{\mathrm{BH}}^{1.42}$) was predicted by Heinz & Sunyaev (2003) for typical flat-spectrum core-dominated radio jets and standard accretion scenarios. They also commented that variations in the other source parameters, like the accretion rate (measured in Eddington units), the viscosity parameter, and viewing angle will only cause a mass-independent scatter around this relation. This theoretical prediction between radio luminosity and black hole mass is not entirely consistent with the linear relation derived in this work (2.3σ difference). However, Heinz & Sunyaev (2003) pointed out that if the jets are powered by spin extraction from the black hole (Blandford & Znajek 1977), then the jet variables will not only depend on ${M}_{\mathrm{BH}}$ and the accretion rate (which was the working assumption in Heinz & Sunyaev 2003) but on the spin as well. It is also discussed in Merloni et al. (2003) that the radio luminosity would be sensitive to the spin of the black hole in the case where the process of jet formation is dependent on the spin. Within the framework of the Blandford–Znajek mechanism (Blandford & Znajek 1977), the luminosity of the jet should be proportional to ${M}_{\mathrm{BH}}:{L}_{\mathrm{BZ}}\propto \dot{m}f{(a)}^{2}{M}_{\mathrm{BH}}$ (Tchekhovskoy et al. 2010; Daly 2016 and references therein), where $\dot{m}$ is the accretion rate in Eddington units and f(a) is the spin function ($f(a)=a/[1+\sqrt{1-{a}^{2}}]$, a being the dimensionless spin parameter). A linear scaling such as the one we find is then expected provided that the product $\dot{m}f{(a)}^{2}$ does not vary significantly between sources and does not depend on ${M}_{\mathrm{BH}}$.

The intrinsic scatter around the best-fit relation may be due to either $\dot{m}$ or f(a) random variations around the blazar mean $\dot{m}f{(a)}^{2}$ product. An additional source of scatter could be variations of the broadband luminosity. Our estimates of the ${{ \mathcal L }}^{\mathrm{Br}-{\rm{B}}}$ are limited by the time span of the observations. Although the 8 year data set of the F-GAMMA program is sufficiently long, it is still possible that exceptional events can occur outside the monitoring period (so that the mean luminosity is actually higher) or alternately a source was unusually active (so that the mean luminosity is actually lower). Cases such as the ones described above are both unlikely and should have a relatively small contribution to the overall scatter.

Our best-fit slope is consistent with those of similar studies on LLAGNs and microquasars (1.14 ± 0.16, Nagar et al. 2002 and 1.23 ± 0.20, Merloni et al. 2003), however, with smaller scatter. The larger scatter found in LLANGs could be the result of mild beaming that has not been accounted for (see Merloni et al. 2003) or other effects related to the use of single-frequency measurements. Beaming effects can also be responsible for the slightly steeper slope in Merloni et al. (2003), since the sources whose beaming is important (blazar-like sources) will have systematically higher luminosities compared to the parent population. Thus the slope will be swifted to stepper values. The fact that our results are consistent with those derived for unbeamed sources (a) provides further support to the fact that we are efficiently correcting for the relativistic effects; and (b) provides supportive evidence for the universality of the derived scaling.

The microquasars that agree best with the relation derived from blazars all share the same accretion state while the one at 3σ is in the intermediate-soft state (see the Appendix). This would suggest that the blazars in our sample are in a similar accretion state as the ones that lie on the best-fit line, i.e., the hard state. The presence of a big blue bump (BBB) in the SED of FSRQs could indicate a geometrically thin, optically thick disk, contrary to the disk of microquasars in the hard state (e.g., Done et al. 2007). Although the optical emission from the boosted jet can dominate over that of the accretion disk and conceal such features in the SED, out of the 20 FSRQs in our sample only 3 sources have a visible BBB. We have verified that excluding these sources from our analysis has no effect on the derived best-fit relation (B remains the same, A = 1.14 ± 0.14). Same accretion states would also suggest that the Eddington ratios (${R}_{\mathrm{Edd}}={L}_{\mathrm{bol}}$/L_Edd) of the supermassive BH sample would be similar to those of the hard-state microquasars. Although estimates of R_Edd are uncertain for beamed sources, recent studies suggest that R_Edd may be less than 10% in the majority of blazars and radio galaxies (Ghisellini et al. 2010; Chen et al. 2015; Xiong et al. 2015). Narrow Line Seyfert 1 are known to have a high R_Edd; however, recent studies suggest that ${R}_{\mathrm{Edd}}\lt 10 \%$ in the γ-ray loud NLS1s as well (Liu et al. 2016). In order to estimate the Eddington ratio for the sources in our sample, we used the Kaspi et al. (2000) relation between Eddington ratio and monochromatic luminosity at 5100 Å, and we estimated R_Edd in nine of them. Although it is not clear whether such relations, which have been established for radio-quiet objects, are also applicable to radio-loud quasars, we found that eight of them may have ${R}_{\mathrm{Edd}}\lt 1 \%$, which is almost identical to the Eddington ratios of the microquasars when in their hard state. The only source with ${R}_{\mathrm{Edd}}\gt 1 \%$ is one of the three sources with a visible BBB in their SED. Although the lack of ${R}_{\mathrm{Edd}}$ estimates prevents us from drawing strong conclusions, it would appear that, on average, the blazars in our sample share the same accretion state with hard-state microquasars. It makes sense to compare blazars and microquasars since, after all, they both have jets. We find that ${{ \mathcal L }}^{\mathrm{Br}-{\rm{B}}}\propto {M}_{\mathrm{BH}}$ for all of them, which is consistent with the predictions of the Blandford and Znajek mechanism. So, if this is the case, and if the accretion regime is different in hard-state microquasars and in some of the blazars in our sample (say powerful FSRQs), then our results suggest that the Blandford and Znajek mechanism can operate, irrespective of the accretion regime (e.g., McKinney 2005; Tchekhovskoy et al. 2010). As discussed above, diff-erences in $\dot{m}$ (i.e., different accretion state) between sources would increase the scatter around the best-fit relation. It will then be interesting to populate our sample with FSRQs with visible and strong BBB, and sources with higher R_Edd, in order to examine if the best-fit model moves closer to the microquasar in the intermediate state and the overall scatter decreases. This would suggest that blazars also have two accretion regimes.

The scaling derived in this work provides for the first time an independent method of estimating the Doppler factor starting from first principles. Estimates of the BH mass and the observer's frame broadband radio luminosity can be used to derive Doppler factor estimates within the scatter of the best-fit relation. These Doppler factor estimates can be used to either reduce the number of free parameters of SED models or, alternately, used to distinguish between acceptable SED models. However, the scaling relation was derived using the mean observer's frame ${{ \mathcal L }}^{\mathrm{Br}-{\rm{B}}}$ based on the F-GAMMA 8 year data set. Given the possibility of changes in the Doppler factor during a significant event (bends in the jet, local acceleration, e.g., Homan et al. 2009, 2015), deriving a Doppler factor from single-epoch observations might not yield a representative estimate for that particular event.

Using the multi-wavelength radio light curves of the F-GAMMA program and the most accurate Doppler factors available to date, we identified a strong linear correlation between the intrinsic broadband radio luminosity of the jets and the black hole mass that extends nine orders of magnitude to stellar-mass black hole systems. Such a universal scaling law is the first ever to bridge observer's and emission rest-frames in beamed sources.

The scaling derived in this work constitutes an important breakthrough in blazar physics. First, it comprises clear evidence of scale-invariant BH jets and of a connection between the properties of supermassive black holes and the large-scale jets they cause in beamed sources. Second, it provides a solid prediction on the ${{ \mathcal L }}_{\mathrm{int},{\rm{c}}}^{\mathrm{Br}-{\rm{B}}}$ of intermediate mass black holes if such exist, and if they form jets. Third, it provides an independent method of estimating the Doppler factor, which will undoubtedly prove an important contribution in constraining SED fitting and the different jet emission models in beamed sources. Fourth, the universality of the scaling would suggest that blazars are in a similar accretion regime as the hard state in microquasars. Finally, our findings point toward the Blandford–Znajek mechanism as the dominant mechanism for jet production in black hole powered jets, and set strong constraints on other potential jet models since they have to reproduce such linear relations.

The authors would like to thank Andreas Zezas and the anonymous referee for comments and suggestions that helped improve this work. This research was supported by the "Aristeia" Action of the "Operational Program Education and Lifelong Learning" and is co-funded by the European Social Fund (ESF) and Greek National Resources, and by the European Commission Seventh Framework Program (FP7) through grants PCIG10-GA-2011-304001 "JetPop" and PIRSES-GA-2012-31578 "EuroCal." Our study is based on observations carried out with the 100 m telescope of the MPIfR (Max-Planck-Institut für Radioastronomie) and the IRAM 30 m telescope. IRAM is supported by INSU/CNRS (France), MPG (Germany), and IGN (Spain). I.N., I.M., and V.K. were supported for this research through a stipend from the International Max Planck Research School (IMPRS) for Astronomy and Astrophysics at the Universities of Bonn and Cologne.

Our sample of microquasars was selected on the availability of contemporaneous multi-wavelength radio observations. Their data were measured from compact synchrotron jets (hard state) for all sources except Cyg X-3 (the microquasar with the highest luminosity in Figure 3) whose data were taken from ballistic ejecta during radio outbursts. Calculating ${{ \mathcal L }}_{\mathrm{int},{\rm{c}}}^{\mathrm{Br}-{\rm{B}}}$ for microquasars is not trivial, the main obstacle being the lack of robust estimates for key parameters (e.g., distance, Doppler factor), as well as the lack of contemporaneous multi-wavelength radio observations. In most cases, the range of available frequencies was confined in the lower frequency range (<43 GHz). The observed flux density at higher frequencies was the result of extrapolation by assuming that the spectral index derived from the highest available frequencies of each microquasar remains constant up to 142.33 GHz (which is the highest available frequency for the supermassive BH sample).

To estimate the relativistic boosting and revert to the rest frame, we gathered for each source all available estimates in the literature for the velocity of the jet, distance, and viewing angle (Table 1). For each source, we created a parameter space set by the extrema of the aforementioned quantities. We then uniformly drew a random value for each of these parameters and calculated ${{ \mathcal L }}_{\mathrm{int},{\rm{c}}}^{\mathrm{Br}-{\rm{B}}}$. We repeated this process 10⁴ times and calculated the mean and spread (minimum, maximum) of ${{ \mathcal L }}_{\mathrm{int},{\rm{c}}}^{\mathrm{Br}-{\rm{B}}}$ for each source. The y-values for the microquasars in Figure 3 are the mean, and the y-axis errorbars the spread for each of the sources that resulted from the random sampling.

We follow a similar procedure for the ${M}_{\mathrm{BH}}$, by gathering all the available estimates in the literature and calculating the mean and spread. The x-values for the microquasars in Figure 3 are the mean, and the x-axis errorbars are the spread for each of the sources.

The microquasars, parameter ranges, and the references to them are the following.

• For SS 433. we used radio data from Trushkin et al. (2003), and parameter values from Seward et al. (1980), Fender (2001), Trushkin et al. (2003), Bowler (2010), Panferov (2010), Cherepashchuk et al. (2013), and references therein.
• For GRS 1915+105, we used radio data from Rodriguez et al. (1995), and parameter values from Mirabel & Rodríguez (1994), Fender (2001, 2003), Fender & Belloni (2004), Miller-Jones et al. (2005), Punsly (2011), Reid et al. (2014), Zdziarski (2014), Punsly & Rodriguez (2016), and references therein.
• For Cyg X-3, we used radio data from Koljonen et al. (2010), and parameter values from Fender et al. (1997), Ling et al. (2009), Hjalmarsdotter et al. (2009), Vilhu et al. (2009), Dubus et al. (2010), Vilhu & Hannikainen (2013), Zdziarski et al. (2012, 2013, 2016), and references therein.
• For GRO J1655-40, we used radio data from Hannikainen et al. (2000), and parameter values from Hannikainen et al. (2000), Trushkin (2000), Mirabel et al. (2002), Fender (2003), Fender et al. (2004, 2010), Narayan & McClintock (2005), Motta et al. (2014), Stuchlík & Kološ (2016), and references therein.
• For V* V404 Cyg, we used radio data from Gallo et al. (2005), and parameter values from Cherepashchuk et al. 2004, Heinz & Merloni (2004), Miller-Jones et al. (2009), Khargharia et al. (2010), Xie et al. (2014), Siegert et al. (2016), and references therein. Since the only available estimate for the jet Lorentz factor (Heinz & Merloni 2004) was a lower limit (${\rm{\Gamma }}\gt 5$, higher than the alternate estimate), given the typical range of Lorentz factors in these sources, we assumed an upper limit of 10. Assuming a smaller value (e.g., ${\rm{\Gamma }}=7$) would result in a difference of maximum values lower than a factor two.