Blazar Variability from Turbulence in Jets Launched by Magnetically Arrested Accretion Flows

It is generally accepted that blazars are active galaxies with relativistic jets aligned close to our line of sight (e.g., Urry & Padovani 1995). Broadly speaking, blazars can be classified into two main categories: low-luminosity BL Lacs and high-luminosity flat-spectrum radio quasars (FSRQs), forming the so-called "blazar sequence" (Fossati et al. 1998; Ghisellini et al. 2017). These objects are variable at all observed wavelengths from radio to γ-rays (e.g., Ulrich et al. 1997). In particular, the γ-ray light curves show variability on timescales ranging from minutes to years. This variability can be characterized as having a power density spectrum (PDS) of power-law shape, spanning the entire observed frequency range. The power-law index is remarkably similar for both BL Lacs and FSRQs. Abdo et al. (2010) report a PDS slope of 1.4 ± 0.1 for nine bright FSRQs, as well as slopes of 1.7 ± 0.3 and 1.5 ± 0.2 for 6 BL Lacs and 13 faint FSRQs, respectively. Recently, Ackermann et al. (2016) report a slope of 1.24 ± 0.15 in the case of the FSRQ 3C 279, and Goyal et al. (2017) report a slope of 1.1 ± 0.2 for the BL Lac object PKS 0735+178. Max-Moerbeck et al. (2014) provide estimates of the γ-ray PDS slopes for 29 blazars, using simulated light curves to properly account for noise processes. Their findings are largely consistent with the results of Abdo et al. (2010). A similar analysis was performed by Ramakrishnan et al. (2015) for a sample of 55 blazars from the first five years of Fermi/LAT data. They find average slopes of 1.3 and 1.1 obtained from 35 FSRQs and 12 BL Lacs, respectively. These slopes are somewhat smaller than those reported by Abdo et al. (2010); however, the slopes reported for 3C 279 are consistent in both cases.

The radio emission, resolved at large radii, clearly originates in jets. Apparent superluminal motion provides compelling evidence that the radio emitting matter propagates with relativistic velocities (e.g., Urry & Padovani 1995). Since the γ-rays are unresolved, however, the source of high-energy emission is uncertain and could potentially be located much closer to the central supermassive black hole (e.g., Pe'er & Markoff 2012). Generally speaking, the observed short timescale variability implies a compact emission region. While such a compact region close to the black hole might be responsible for variability in low-luminosity systems, high-luminosity systems (such as 3C 279) require a significant Lorentz factor to overcome the pair opacity, which would otherwise prevent γ-rays from escaping to infinity (see, e.g., the Appendix and Dondi & Ghisellini 1995). The location of the γ-ray emission zone remains a topic of active research (e.g., Madejski & Sikora 2016, for a recent review).

The shortest variability timescales are comparable to, and in some cases shorter than, the light-crossing time of the black hole (Aharonian et al. 2007; Albert et al. 2007; Aleksić et al. 2011; Ackermann et al. 2016). For example, Aharonian et al. (2007) observed variability on a timescale of ${t}_{\mathrm{var}}\approx 200s$ during a flare in PKS 2155–304. This is more than an order of magnitude shorter than the corresponding light-crossing time, inferred by the empirical relation between galactic bulge luminosity and black hole mass (Bettoni et al. 2003; Aharonian et al. 2007). The mechanism responsible for this fast variability is poorly understood. Popular models include the "jets in a jet" model (Giannios et al. 2009), in which magnetic reconnection in highly magnetized regions of the jet accelerates compact blobs of plasma to relativistic velocities in the bulk jet frame; the magnetospheric acceleration model (Levinson & Rieger 2011), in which charged particles are accelerated by unscreened electric fields in a charge starved vacuum gap of the black hole magnetosphere; the relativistic turbulence model (Narayan & Piran 2012), in which magnetohydrodynamic turbulence in the jet produces compact blobs on scales smaller than the horizon radius, similar to those in the "jets in a jet" scenario; and the jet-star interaction model (Barkov et al. 2012) in which stars cross the jet close to the black hole.

In this work, we argue in favor of the proposal by Narayan & Piran (2012), namely that turbulence in the relativistic jet can produce the observed variability properties. To support this claim, we investigate the variability of high-energy radiation from a magnetically arrested accretion flow (MAD; Narayan et al. 2003), which efficiently launches jets via the Blandford–Znajek mechanism (Blandford & Znajek 1977; Tchekhovskoy et al. 2011; McKinney et al. 2012). Importantly, computational limitations force us to restrict our analysis to the inner $r\lesssim 200{r}_{g}$ of the MAD, where ${r}_{g}={GM}/{c}^{2}$ is the gravitational radius. O' Riordan et al. (2016a) showed that, for rapidly rotating black holes, observers see deep into the hot dense, highly magnetized plasma in the inner parts of MADs. In what follows, we will refer to this region as the "jet launching region." Since the radiated power in our model is dominated by the turbulent plasma of the jet launching region close to the horizon, we expect the resulting radiation to be variable on timescales comparable to the light-crossing time ${t}_{g}={r}_{g}/c$.

We find that γ-ray light curves from the jet launching region can be described by a PDS, which is remarkably similar to that observed in both BL Lacs and FSRQs. Furthermore, the optical synchrotron emission can also be described by the same PDS, despite having a different origin. The large inferred pair opacity in high-luminosity FSRQs prevents radiation originating in the jet launching zone from directly contributing to the observed variability. However, assuming that the turbulent properties of the launching region persist at large radii (and therefore large bulk Lorentz factor Γ), we argue that it is plausible that the observed variability properties result from turbulence in the jet.

In Section 2, we briefly describe our model and assumptions. In Section 3, we show the resulting PDS and compare with observations. In Section 4, we discuss our findings, with emphasis on the limitations of our model and suggestions for future work in this area. We use units where $G=c=1$ and therefore ${r}_{g}={t}_{g}=M$; however, we occasionally reintroduce factors of c for clarity.

We simulate a MAD accretion flow using the fully 3D general-relativistic magnetohydrodynamic (GRMHD) code HARM (Gammie et al. 2003). Our model is based on run A0.99N100 from McKinney et al. (2012), restarted at $t=15000\,M$ with a very high temporal resolution of ${\rm{\Delta }}t=0.1\,M,$ and a spatial resolution of ${N}_{r}\times {N}_{\theta }\times {N}_{\phi }=288\times 128\times 128$. This is the highest spatial resolution available for a 3D GRMHD simulation of a MAD accretion flow. As described in Tchekhovskoy et al. (2011) and McKinney et al. (2012), the grid concentrates resolution in the equatorial disk at small radii and in the jet at large radii. Close to the black hole, the radial grid is logarithmically spaced with ${\rm{\Delta }}r\approx 0.03\,M$ at the horizon, increasing to ${\rm{\Delta }}r\approx M$ at $r=30\,M$. For our purposes, we consider a black hole spin of a = 0.99 and limit our analysis to times $t\geqslant {\rm{15,000}}\,M$, well after the simulation settled into a quasi-steady MAD state. As is discussed in O' Riordan et al. (2016a, 2016b), we remove the numerical density floor material from the center of the funnel region before performing our radiative transport post-processing calculation. Although this material is required to maintain numerical stability during the GRMHD simulation, it can become artificially dense and hot and thus might distort the resulting spectra. Therefore, when calculating the spectra, we simply remove the floor material, focussing our attention on the self-consistent disk and funnel wall regions. We do this by setting the density to zero in regions where the ratio of magnetic and rest mass energy densities becomes too large; ${b}^{2}/\rho \gt \zeta$. Here, ${b}^{2}={b}^{\mu }{b}_{\mu },{b}^{\mu }$ is the magnetic 4-field, ρ is the rest mass density of the fluid. At the horizon, we choose $\zeta =20$, and linearly interpolate to $\zeta =10$ at $r=10\,M$. Beyond this, we choose a constant value of $\zeta =10$. This ensures that we do not remove any of the self-consistent, highly magnetized fluid close to the horizon when discarding the floor material.

We calculate the properties of the resulting radiation field using a code based on grmonty (Dolence et al. 2009). We include synchrotron emission, self-absorption, and Compton scattering from relativistic, thermal electrons. We choose a constant proton-to-electron temperature ratio ${T}_{p}/{T}_{e}=30$, which is consistent with recent radiative GRMHD simulations of active galactic nuclei (Sa̧dowski & Gaspari 2017). As discussed in Ressler et al. (2015), Foucart et al. (2016), and Sa̧dowski et al. (2017), electrons in the jet and disk are probably subject to different heating and cooling mechanisms due to differences in density and magnetization. This would likely cause ${T}_{p}/{T}_{e}$ to vary across these regions. In our case, the emission is strongly dominated by a highly magnetized, compact region close to the black hole and we find that varying ${T}_{p}/{T}_{e}$ independently in the jet and disk has little effect on our main results. Furthermore, choosing a different ratio of ${T}_{p}/{T}_{e}=3$ everywhere does not affect the resulting temporal behavior. The absence of non-thermal electrons, which could be important in a more realistic description of the high-energy γ-ray light curves, is a limitation of the current model; we will address this in a future work.

For analyzing the variability, we choose snapshots of the GRMHD data corresponding to time steps ${\rm{\Delta }}t$ ranging from $0.1\,M$ to $10\,M$. This allows us study variability on a wide range of timescales, while still producing light curves that are evenly sampled in time. We also vary the total number of photons tracked during our radiative transport calculation. The number of photons tracked per time step ranges from 10⁵, for calculations over long timescales, to 10⁷ per time step for the shortest timescales. Importantly, we use a "fast-light" approximation in which each snapshot is treated independently. This allows us to parallelize our radiative transport calculation over both photons and time steps, significantly improving both performance and computational simplicity. Since the emission from our model originates in the compact jet launching region, this approximation is valid over a wide range of timescales. The disadvantage of this approach, however, is that it sets practical limitations on the shortest timescales that we can reliably probe. The "fast-light" assumption effectively sets the speed of light to infinity, which fails to be valid for timescales shorter than the light-crossing time of the emission region. Therefore, our current code is not sensitive to very fast variability on timescales shorter than $\sim M$. Furthermore, the smallest structures in the GRMHD data are smoothed out over ∼10 cells and so, on the shortest timescales, our analysis is also limited by the grid resolution.

For a thermal electron distribution, the synchrotron power radiated in the comoving frame scales as ${P}_{\mathrm{syn}}\sim \int {dV}\,{{nB}}^{2}{{\rm{\Theta }}}^{2}$, where ${dV}=\sqrt{-g}{{dx}}^{r}{{dx}}^{\theta }{{dx}}^{\phi }$, $g=\det ({g}_{\mu \nu })$ is the metric determinant, n is the electron number density, B is the magnetic field, ${\rm{\Theta }}={kT}/{{mc}}^{2}$ is the electron temperature, and all fluid quantities are measured in the fluid frame. In Figure 1, we show a snapshot of the fluid data at $t={\rm{15,700}}\,M$. The white region aligned with the spin axis corresponds to the numerical floor material that has been removed from the funnel. We plot the quantity ${{nB}}^{2}{{\rm{\Theta }}}^{2}$, which is proportional to the comoving synchrotron power per unit volume. Clearly, the synchrotron emission in this model is strongly dominated by the inner $\sim 5\,M$. Similarly, the synchrotron self-Compton component originates in a compact region close to the black hole. Therefore, the variability in the resulting light curves is dominated by the turbulent, highly magnetized plasma in the inner parts of the MAD. Since the emission is dominated by fluid at small radii, there is no significant Lorentz boosting of the radiation.

Zoom In Zoom Out Reset image size

In the top panel of Figure 2, we show a snapshot of the comoving electron number density close to the black hole horizon. In the bottom panel, we show the corresponding comoving magnetic field strength. The flow is turbulent, and there are significant fluctuations in the fluid properties extending down to spatial scales much smaller than M. Notably, there are structures in the density apparent on scales $\lesssim 0.5\,M$, and structures in B apparent on scales of $\lesssim 0.1\,M$ due to a polarity change of the poloidal magnetic field. The resolution of our numerical model means that we can only resolve such small-scale structures close to the horizon. However, we expect that similar sub-M features could be produced by turbulent processes operating at large radii. In particular, if such inhomogeneities are produced at large distances in a highly magnetized, relativistic jet, it is plausible that they could contribute significantly to variability on very short timescales.

Zoom In Zoom Out Reset image size

The PDS for real signal h(t) is $P(f)=2{| \hat{h}(f)| }^{2}$, where $\hat{h}(f)$ is the Fourier transform of h(t) and the frequency range is $0\leqslant f\lt \infty$ (see, e.g., Press et al. 1986). To numerically estimate the PDS, we follow a similar procedure to that described in Wellons et al. (2014). To reduce noise, we divide each light curve into two segments and average the resulting power spectra. Furthermore, to suppress leakage between frequency bins, we apply a Hann window to each segment (e.g., Press et al. 1986). In practice, we use the implementation provided by the "welch" function in the SciPy Python module (Jones et al. 2001). In Figure 3, we show the PDS calculated from four γ-ray light curves (frequency integrated between 10¹⁹ and ${10}^{24}\,\mathrm{Hz}$) with time steps ranging from $0.1\,M$ to $10\,M$. The resulting PDS is not particularly sensitive to the sampling of the light curves, with good agreement in regions of overlap. These γ-rays result from inverse-Compton scattering by hot electrons in the jet launching region. For comparison with observations, the bottom axis shows the frequency in units of ${\mathrm{day}}^{-1}$, while the top axis shows the frequency in units of inverse light-crossing time, normalized to a black hole mass of $M\,=5\times {10}^{8}\,{M}_{\odot }$. We have also included data for 3C 279 from Ackermann et al. (2016) as well as the average PDS for six BL Lacs and nine FSRQs from Abdo et al. (2010). Note that the scaling of the top axis is appropriate for comparison with the 3C 279 data, while the Abdo et al. (2010) data is averaged over different masses. Remarkably, our simulated PDS is largely consistent with the data from both BL Lacs and FSRQs (including 3C 279) below $f\sim \,\mathrm{few}\,{\mathrm{day}}^{-1}$. This result is suggestive of a connection between blazar variability and turbulence in the jet launching region of MADs. Note that since we are primarily interested in the scaling with f, the normalization of the y-axis is arbitrary, and simply chosen to ease comparison between the shape of our simulated PDS and the data. Also, since the overall magnitude of the PDS is probably sensitive to the prescription for treating the electron temperature, a detailed investigation of the normalization is beyond the scope of the this work.

Zoom In Zoom Out Reset image size

The rough power-law behavior below $f\sim \,\mathrm{few}\,{\mathrm{day}}^{-1}$ changes to a cutoff at high frequencies. The location of this cutoff is likely affected both by the spatial resolution of the fluid data, and the photon statistics of the radiative transport calculation. Furthermore, the spectrum transitions to a flat white noise section above $10\,{\mathrm{day}}^{-1}$. This is clearly due to poor photon statistics on the shortest timescales. That is, the magnitude of the γ-ray variability is overwhelmed by fluctuations due to photon statistics at timescales shorter than $\sim M$. In Figure 4, we show one of the γ-ray light curves with ${\rm{\Delta }}t=0.1\,M$. Clearly, the variability on short times is significantly affected by noise in the radiative transport calculation. The ${\rm{\Delta }}t=0.1\,M$ light curve corresponds to our highest resolution simulations, tracking $\sim {10}^{7}$ photons per time step, and thus improving upon this would be too computationally expensive. In any case, our "fast-light" assumption limits the reliably of the radiative transport results on timescales shorter than $\sim M$.

Zoom In Zoom Out Reset image size

In Figure 5, we show the PDS in the optical band. In this case the radiation is primarily due to synchrotron emission. Interestingly, the variability in the optical band is the same as that in the γ-rays. Recently, Goyal et al. (2017) found that the γ-ray PDS in PKS 0735+178 is somewhat flatter than the corresponding PDS in the optical and radio bands. They speculate that this discrepancy is due to additional stochastic processes, which only affect the inverse-Compton component. We will investigate this interesting observation in a future work.

Zoom In Zoom Out Reset image size

As a step toward understanding the origin of variability in blazars, we calculated the temporal dependence of optical and γ-ray radiation from the jet launching region in a MAD accretion flow. In this work, we are concerned with two main properties of the observed variability. Namely, (1) that the variability can be described by a PDS of power-law shape across the entire observed frequency range and (2) the very fast variability observed in some sources, which can be significantly shorter than the light-crossing time of the supermassive black hole. We argue that turbulence in the highly magnetized plasma of a relativistic jet can plausibly account for both of these properties.

The radiation from our model is dominated by turbulent plasma close to the black hole. The optical band is dominated by synchrotron emission, while the γ-rays are produced by inverse-Compton scattering from hot electrons. Both the optical and γ-ray light curves show variability, which can be characterized by a PDS with a power-law shape below $f\sim \,\mathrm{few}\,{\mathrm{day}}^{-1}$. Remarkably, the power-law section of our simulated PDS quantitatively reproduces the PDS observed in both BL Lacs and FSRQs. This is suggestive of a connection between turbulent processes in MADs and the observed variability properties of blazars. Furthermore, we find the same power-law index for both the optical and γ-ray light curves, implying that the variability properties should not be particularly sensitive to the underlying emission mechanism.

We also showed that the turbulent flow in the jet launching region naturally produces structures on spatial scales much smaller than the horizon radius. Such small-scale structures, if produced at large distances in the highly magnetized, relativistic jet, would likely produce significant short timescale variability, though it is indeed possible that the extremely fast variability in some sources might require additional microphysics to inject energy into the plasma on small scales (e.g., Giannios et al. 2009; Narayan & Piran 2012). Based on these results, and assuming that the turbulence in the jet launching region persists to large radii, we argue that turbulence in the relativistic jet is responsible for the observed variability in blazars over many orders of magnitude in frequency.

While we have focussed on blazars, we expect the results presented here to also be applicable in a range of systems that can be modeled as a radiatively inefficient accretion flow (RIAF). In particular, our results could be appropriate for Sgr A* and the low-hard state in XRBs, though observing variability in XRBs on timescales comparable to the light-crossing time of the black hole would be extremely challenging. It would also be interesting to see how our results might change in the super-Eddington MAD regime, which could be important for describing ultra-luminous X-ray sources (Sa̧dowski et al. 2013; McKinney et al. 2014, 2015; Sa̧dowski & Narayan 2015; Narayan et al. 2017).

Our model has some limitations that prevent direct application in general, especially to systems with large luminosities. First, while our model might be applicable in its current form to some low-luminosity systems, the fact that the radiation originates from close to the black hole means that it is likely inappropriate for directly explaining the observed γ-rays from high-luminosity sources. As discussed in the Appendix, the γ-ray radiation from high-luminosity systems (such as 3C 279) cannot originate from close to the black hole since this region is too compact. Such a compact region would be subject to a very large pair opacity and so the emitted γ-rays would not escape to infinity. This implies that the γ-rays likely originate from a region with a significant bulk Lorenz factor Γ, i.e., a relativistic jet. If we assume that the intrinsic variability in the jet at large radii matches our results for the inner jet launching region, then the observed variability should follow the same power law, though shifted to higher frequencies by the Doppler factor ${ \mathcal D }={({\rm{\Gamma }}(1-\beta \cos \theta ))}^{-1}$. The Lorentz boosting would also affect the overall normalization of the PDS, which we have not considered here since the scaling of the y-axes in Figures 3 and 5 is arbitrary.

Second, the current work is limited by our assumption of a thermal distribution of electrons and by our simplified treatment of the electron thermodynamics. In reality, the electron distribution in highly magnetized regions of the plasma likely has a high-energy non-thermal tail due to acceleration either by shock waves (e.g., Sironi et al. 2015) or magnetic reconnection (e.g., Sironi & Spitkovsky 2014). These non-thermal electrons would contribute significantly to the observed radiation at very high energies. Furthermore, the details of the electron physics in RIAFs (e.g., Narayan & McClintock 2008; Yuan & Narayan 2014) remains a highly active area of research (Ressler et al. 2015; Foucart et al. 2016; Sa̧dowski et al. 2017). We will incorporate more complicated models of the electron physics in a future work. Therefore, although our model reproduces the variability properties quite well, the lack of non-thermal particles and large bulk Lorentz factor means that we cannot quantitatively reproduce the wide range of observed blazar spectra.

As mentioned previously, we have used a "fast-light" approximation in which the snapshots of the fluid data are treated as time-independent during the radiative transport calculation. We will relax this approximation in a future work in order to study variability on shorter timescales in greater detail.

We conclude by noting that a thorough investigation of the γ-ray zone in blazars will require more detailed microphysical modeling for treating the electrons, as well as global GRMHD simulations of self-consistent jet launching and propagation to large distances, with sub-M grid resolution at large radii. Recent MHD simulations investigated the large-scale structure of galactic jets propagating in an ambient medium (Barniol Duran et al. 2017; Tchekhovskoy & Bromberg 2016); however, resolving small-scale structures in both the inner accretion flow and the jet at large distances is currently too computationally expensive.

The authors acknowledge the DJEI/DES/SFI/HEA Irish Centre for High-End Computing (ICHEC) for the provision of computational facilities under projects ucast007c and ucast008b. M.O.R. is supported by the Irish Research Council under grant number GOIPG/2013/315. This research was partially supported by the European Union Seventh Framework Programme (FP7/2007-2013) under grant agreement no 618499. J.C.M. acknowledges NASA/NSF/TCAN (NNX14AB46G), NSF/XSEDE/TACC (TGPHY120005), and NASA/Pleiades (SMD-14-5451). The authors would like to thank the anonymous referee for helpful and constructive comments on the manuscript.

Software: HARM (Gammie et al. 2003), grmonty (Dolence et al. 2009), SciPy (Jones et al. 2001).

Here, we estimate constraints on the luminosity from the pair opacity. The following applies in the comoving frame of the source and follows closely the treatment of Dondi & Ghisellini (1995). High-energy γ-rays can collide with softer radiation to produce ${e}^{\pm }$ pairs. The cross-section for this process is maximized for collisions between γ-rays of energy $x=h\nu /{{mc}}^{2}$ and target photons of energy ${x}_{t}=1/x$. This maximum cross-section is $\sigma =3{\sigma }_{T}/16$ (e.g., Lang 1980), where ${\sigma }_{T}$ is the Thomson cross-section. Assuming that the source is spherical and emits isotropically, the corresponding optical depth can be written in terms of the luminosity as

$$\begin{eqnarray}&&{\tau }_{\gamma \gamma }(x)=\displaystyle \frac{3{\sigma }_{T}}{16}\displaystyle \frac{L({x}_{t})}{4\pi {{mc}}^{3}R}\end{eqnarray} \tag{ 1 }$$

Therefore, the condition that ${\tau }_{\gamma \gamma }(x)\lesssim 1$ constrains the luminosity of the soft radiation field to be

$$\begin{eqnarray}&&L({x}_{t})\lesssim 3\times {10}^{43}\,{M}_{8}\left(\displaystyle \frac{R}{{r}_{g}}\right)\mathrm{erg}\ {{\rm{s}}}^{-1}\end{eqnarray} \tag{ 2 }$$

where we have written the luminosity in terms of the gravitational radius ${r}_{g}\approx 1.5\times {10}^{13}\,{M}_{8}\,\mathrm{cm}$, and M₈ is the black hole mass in units of ${10}^{8}\,{M}_{\odot }$.

In the case of 3C 279, the observed γ-rays primarily interact with a soft (X-ray) radiation field of luminosity $L\sim 3\,\times {10}^{47}\,\mathrm{erg}\ {{\rm{s}}}^{-1}$ (Ackermann et al. 2016). Adopting a black hole mass of $M=5\times {10}^{8}\,{M}_{\odot }$ and an emission region of size $R=5\,{r}_{g}$, we find that the luminosity of the target radiation is constrained to be below $L\lesssim {10}^{45}\,\mathrm{erg}\ {{\rm{s}}}^{-1}$. The obsereved luminosity is more than two orders of magnitude larger than this. Combined with the short variability timescale, we are forced to conclude that the source of emission is strongly Doppler boosted and thus likely originates in a relativistic jet.