Probing Acceleration and Turbulence at Relativistic Shocks in Blazar Jets

Probing Acceleration and Turbulence at Relativistic Shocks in Blazar Jets


Diffusive shock acceleration (DSA) at relativistic shocks is widely thought to be an important acceleration mechanism in various astrophysical jet sources, including radio-loud active galactic nuclei such as blazars. Such acceleration can produce the non-thermal particles that emit the broadband continuum radiation that is detected from extragalactic jets. An important recent development for blazar science is the ability of Fermi-LAT spectroscopy to pin down the shape of the distribution of the underlying non-thermal particle population. This paper highlights how multi-wavelength spectra spanning optical to X-ray to gamma-ray bands can be used to probe diffusive acceleration in relativistic, oblique, magnetohydrodynamic (MHD) shocks in blazar jets. Diagnostics on the MHD turbulence near such shocks are obtained using thermal and non-thermal particle distributions resulting from detailed Monte Carlo simulations of DSA. These probes are afforded by the characteristic property that the synchrotron peak energy does not appear in the gamma-ray band above 100 MeV. We investigate self-consistently the radiative synchrotron and inverse Compton signatures of the simulated particle distributions. Important constraints on the diffusive mean free paths of electrons, and the level of electromagnetic field turbulence are identified for three different case study blazars, Mrk 501, BL Lacertae and AO 0235+164. The X-ray excess of AO 0235+164 in a flare state can be modelled as the signature of bulk Compton scattering of external radiation fields, thereby tightly constraining the energy-dependence of the diffusion coefficient for electrons. The concomitant interpretations that turbulence levels decline with remoteness from jet shocks, and the probable significant role for non-gyroresonant diffusion, are posited.

Acceleration of particles, plasmas, shock waves, turbulence, galaxies: active, galaxies: jets, X-rays, gamma-rays.

1 Introduction

Extragalactic jets of collimated relativistic outflows are some of the most powerful emitters of radiation in the Universe. Two contrasting types of jets are found in transient gamma-ray bursts (GRBs), and persistent but highly-variable active galactic nuclei (AGN). Both are considered prime candidates for the production of ultra-high energy cosmic rays above eV, and perhaps also the high energy neutrinos recently detected by IceCube (Aartsen et al., 2014). Yet they are very different in their origin. GRBs probably result from the explosion of massive progenitor stars in hypernovae, or the merger of compact binary neutron stars, whereas AGN are continually powered over eons by material accreted onto supermassive black holes (SMBHs) in the center of distant galaxies. The masses of such SMBHs are typically in the range (Bentz & Katz, 2015),3 deduced in part from their large luminosities, typically erg/sec, and also from reverberation mapping techniques (e.g., Peterson et al., 2004). In this paper, the focus is on the interpretation of the jet environments of AGN, in particular on the specific subset known as blazars, which exhibit flares with short timescale variations in radio, optical, X-ray and -ray wavebands. The class of blazars was identified following the discovery (Hartmann et al., 1992; Lin et al., 1992) of transient gamma-ray emission in 3C 279 and Mrk 421 by the EGRET instrument on the Compton Gamma-Ray Observatory in the 100 MeV – 1 GeV range. This was around the same time that the Whipple atmospheric Čerenkov telescope (ACT) discerned that Mrk 421 also emitted at TeV energies (Punch et al., 1992).

Blazars generally evince the key identifying feature of their parent BL Lac objects like BL Lacertae of a general absence of emission lines in their non-thermal optical spectra. Their radio continua often possess quite flat spectra. These characteristics suggest synchrotron radiation from non-thermal electrons gyrating in magnetic fields as the origin of their emission in these two wavebands. Such a hypothesis is supported by detections of polarization in both radio and optical bands. Most remarkable among the ensemble of blazar polarimetry studies is the observation of rapid optical polarization angle swings in 3C 279 (seen also for other blazars, including PKS 1510-089; Marscher et al., 2010), contemporaneous with a gamma-ray flare seen by Fermi-LAT (Abdo et al., 2010b). This indicates a large-scale ordering of the magnetic field in the 3C 279 jet, and spatial coincidence of the optical and gamma-ray emission zones. In many blazars, this putative synchrotron component extends to X-rays, where polarization measurements in the not-to-distant future define a hope for further constraining the emission mechanism of blazars (e.g., Krawczynski et al., 2011).

The prevailing paradigm is that blazars’ gamma-ray signals are generated by inverse Compton scattering of synchrotron photons by the same relativistic electrons that emit this radio-to-X-ray signal, so-called synchrotron-self-Compton (SSC) models (e.g., Maraschi et al., 1992; Mastichiadis & Kirk, 1997; Chiang & Böttcher, 2002). An alternative possibillity for the target photons seeding this upscattering is an IR/optical/UV photon source of possibly disk or ambient origin (e.g., Dermer, Schlickeiser & Mastichiadis, 1992; Sikora, Begelman & Rees, 1994) not too distant from the black hole. The main competing scenario for -ray production in blazars is the hadronic model (e.g., Mannheim & Biermann, 1989; Rachen & Mészáros, 1998; Mücke & Protheroe, 2001) where non-thermal protons collide with disk and jet-associated synchrotron photons, generating pair cascades via photo-pion/muon production processes. Proton synchrotron radiation can also appear in the -ray band. High energy neutrinos are a signature product of such cascades, forging a connection of hadronic blazar models to the observation of energetic neutrinos up to around eV by the IceCube experiment (Aartsen et al., 2014).

Blazar spectra observed above 300 GeV by ACTs are typically quite steep, defining a turnover. These photons are subject to strong pair absorption in propagating to Earth due to the presence of intergalactic infra-red and optical starlight background fields. Correcting for such attenuation can yield inferences of very flat intrinsic source spectra (e.g., Stecker, Baring, & Summerlin, 2007), yet uncertainties in the background fields permeate such protocols. An important recent development for blazar science has been the improvement of sensivity in the 100 MeV-100 GeV window, afforded by the Fermi-Large Area Telescope (LAT). Over the last eight years, LAT data has enabled measurements of the power-law index of spectra from numerous blazars (Abdo et al., 2009, 2010c). This provides a more robust measure of the underlying non-thermal particle population, since the LAT band is generally subject to only small pair attenuation, and negligibly so below around 3 GeV. This is a new probe enabled for multiwavelength blazar studies, a tool that we exploit to some extent here. Yet more particularly, one of our case studies, AO 0235+164, is not yet detected with current ACT facilities, so that the Fermi-LAT data prove critical to our interpretation of its jet environment.

The rapid flux variability seen in radio, X-ray and GeV–TeV flares (e.g., see Maraschi et al., 1999; Takahashi et al., 2000, for Mrk 421) drives the prevailing picture for the blazar environment: their jets are relativistic, and compactly structured on small spatial scales that are unresolvable by present gamma-ray telescopes. For a survey of radio and X-ray jet properties including spatial morphology, see Harris & Krawczynski (2006). The supersonic outflows in these jets naturally generate relativistic shocks, and these can form the principal sites for dissipation of the ballistic kinetic energy via acceleration of electrons and perhaps also ions to the ultrarelativistic energies demanded by the X-ray and -ray data. Diffusive Fermi acceleration at such jet shocks is believed to be the main candidate mechanism for energizing such charges (Blandford & Eichler, 1987; Drury, 1983; Jones & Ellison, 1991). This is motivated by the fact that Fermi acceleration is both extremely efficient and very fast, precipitating acceleration rates of the order of the gyrofrequency . This follows because the diffusive collisions of electrons and ions with MHD turbulence in jets is putatively dominated by gyroresonant interactions. The efficiency of particle-turbulence interactions is an issue that we will visit in this paper. We note that shocks also may span a large surface cross section of the jet that the outward-propagating plasma may encounter with high probability.

Another possibility is that shear layers encapsulating sharp velocity gradients transverse to the net flow may precipitate Fermi-type acceleration (Ostrowski, 1990) due to transport of charges straddling the shear “discontinuity.” Observational evidence supporting such transverse velocity structure comes from parsec-scale limb-brightening of blazar and radio galaxy jets revealed in VLBI observations (e.g., see Giroletti et al., 2004, for Mrk 501). MHD/hydrodynamic simulations of spine-sheath jets and their launching (e.g., Meliani & Keppens, 2007; McKinney & Narayan, 2007; Mizuno, Hardee & Nishikawa, 2007; Tchekhovskoy, McKinney & Narayan, 2008) indicate that the sheath, in combination with a poloidal magnetic field, aids in stabilizing the jet. Rayleigh-Taylor-type instabilities may develop at the spine-sheath interface, and such turbulence offers another promising avenue for accelerating relativistic particles. Indeed, particle-in-cell plasma simulations of relativistic shear flows (Liang, Böttcher & Smith, 2013) indicate substantial energization of pairs in the boundary layers. Other paradigms such as acceleration by reconnection of magnetic fields embedded in Poynting-flux dominated outflows can be envisaged. Whether shocks or shear boundaries or reconnection zones provide the dominant energization site for non-thermal particles in blazars defines a major objective for future theoretical efforts within the blazar community.

Our considerations here will focus on acceleration at blazar jet shocks. To explore the blazar shock paradigm, the standard practice has been to develop multiwavelength (MW) spectral models, spanning radio to TeV gamma-ray wavelengths (Dermer, Schlickeiser & Mastichiadis, 1992; Sikora, Begelman & Rees, 1994; Maraschi et al., 1992; Mastichiadis & Kirk, 1997; Chiang & Böttcher, 2002; Böttcher, Dermer & Finke, 2008). This generally gives a broad-brush assessment of a range of jet environment parameters, but the spectral fits in any one band are of limited quality. The high statistics spectroscopy afforded by the Fermi-LAT data demands a closer look at the constraints observations can elucidate for the shocked environs of blazar jets and on the shock acceleration process itself. In particular, theoretical studies of shock acceleration have determined (e.g., Ellison, Jones & Reynolds, 1990; Kirk & Heavens, 1989; Ellison & Double, 2004; Summerlin & Baring, 2012) that a wide variety of power-law indices for accelerated charges are possible for a given velocity compression ratio across the discontinuity. These indices are sensitive to the orientation of the mean magnetic field, and the character of the in situ MHD turbulence. In addition, it has been understood for nearly two decades (see Inoue & Takahara, 1996) that diffusive shock acceleration is so efficient that low levels of field turbulence are required in blazar jets to accommodate synchrotron spectral peaks appearing in the X-ray band. This serves as a central issue to the discourse of this paper.

In the multi-wavelength blazar models developed herein, detailed results from comprehensive Monte Carlo simulations of diffusive acceleration at relativistic shocks are employed, building upon the exposition of Summerlin & Baring (2012). These simulations capture the relationship between turbulence parameters and the power-law index, and also the connection between the thermal bulk of the population (a hot Maxwellian-like component) and the power-law tail of the accelerated species. Such an approach goes beyond an elemental description of the distribution of accelerated charges, usually a power-law truncated at some minimum and maximum particle Lorentz factors, that is commonly invoked in blazar spectral modeling. Our shock simulation approach is outlined in Section 2. We fold these simulation results through the one-zone SSC/external Compton models of Böttcher et al. (2013) for radiation emission and transport in blazars. The radiation modeling protocols are summarized in Section 3, and subsequently results from this analysis are presented therein.

Diagnostics are obtained for the particle mean free path and the level of field turbulence for our three chosen blazars, namely Mrk 501, BL Lacertae, and the BL Lac object AO 0235+164. The non-thermal emission components are primarily generated by leptons that undergo repeated drift acceleration and interspersed upstream reflections in the shock layer, the result being extremely flat distributions emerging from the acceleration zones. By exploring a range of dependences of the diffusive mean free path on the momentum of accelerated charges, the possible interpretation that turbulence levels decline with remoteness from a shock is identified, perhaps signalling a significant role for non-gyroresonant diffusion in the vicinity of blazar jet shocks, an interpretation discussed in Section 4. This determination is impelled by the fact that gyroresonant diffusive acceleration and energization in magnetic reconnection zones are just too fast to accommodate the multi-wavelength spectra of blazars, in particular the positioning of the synchrotron peak in the optical or X-ray bands.

2 Simulations of Diffusive and Drift Acceleration at Shocks

To construct a more precise assessment of multiwavelength emission models for blazars, and derive a deeper understanding of their jet physics, it is necessary to go beyond simplistic truncated power-law distributions for the radiating leptons and hadrons. This is our approach in this paper, where detailed modeling of particle distribution functions, anisotropies, and their relationship to plasma turbulence are afforded through simulations of shock acceleration processes. Simulations are the most encompassing technique to apply to the blazar jet problem because the shocks are inherently relativistic: faster and slower regions of the jet flow travelling near have relative velocities that are mildly-relativistic. This domain will form the focus of our exposition here.

2.1 Background on Shock Acceleration Theory

Various approaches have been adopted to model diffusive shock acceleration (DSA) at relativistic discontinuities. These include analytic methods (e.g., Peacock, 1981; Kirk & Schneider, 1987; Kirk & Heavens, 1989; Kirk et al., 2000), and Monte Carlo simulations of convection and diffusion (e.g., Ellison, Jones & Reynolds, 1990; Bednarz & Ostrowski, 1998; Ellison & Double, 2004; Niemiec & Ostrowski, 2004; Summerlin & Baring, 2012). Particle-in-cell (PIC) plasma simulations have also been employed to study this problem via modeling the establishment of electromagnetic turbulence self-consistently (e.g., Gallant et al., 1992; Silva et al., 2003; Nishikawa et al., 2005; Spitkovsky, 2008; Sironi & Spitkovsky, 2009) in conjunction with the driver charges and their currents. Important insights are gleaned from each of these techniques. All exhibit the core property that in collisionless shocks, non-thermal, charged particles gain energy by scattering between MHD turbulence (for example magnetic “islands” observed in PIC simulations) that is “anchored” in the converging upstream and downstream plasmas, the so-called Fermi mechanism. This defines a fundamental association between turbulence contained in shock environs, diffusion and ultimate acceleration.

Each of these approaches to the shock acceleration problem has both merits and limitations. Analytic techniques provide core insights into global character, though are often restricted to treating particles well above thermal energies where the acceleration process possesses no momentum scale. Monte Carlo simulations with prescribed injected turbulence (e.g., Niemiec & Ostrowski, 2004) explore wave-particle interactions and diffusive acceleration well above thermal energies, but do not treat the electrodynamics of wave generation and dissipation self-consistently. Particle-in-cell simulations provide the greatest depth in modeling shock layer microphysics by solving Maxwell’s equations and the Newton-Lorentz force equation, though their macro-particle approximation does eliminate electrodynamic information on the smallest scales. With increased box sizes, PIC codes have now realized the establishment of truly non-thermal components, but still cannot explore acceleration beyond modest non-thermal energies (e.g., Sironi & Spitkovsky, 2009); their focus is still on the thermal dissipation and injection domains in shock layers. The Monte Carlo approach employed here subsumes the microphysics in a parametric description of diffusion, and so does not investigate the feedback between charges, currents and hydromagnetic waves. Yet it is ideal for describing the diffusive and shock drift acceleration of particles from thermal injection scales out to the highest energies relevant to blazar jet models, and so is the preferred tool for interfacing with astronomical emission datasets.

A key feature of both relativistic and non-relativistic shock acceleration theory is that the acceleration process possesses no momentum scale, and the resulting particle distribution takes the form . For non-relativistic shocks, since their speeds far exceed (), the upstream (downstream) flow speed component in the co-ordinate direction normal to the shock, the energetic particles are nearly isotropic in all fluid frames. The acceleration process then establishes (e.g., Drury, 1983; Jones & Ellison, 1991) a power-law distribution with index , where is the shock’s velocity compression ratio. The index in this limit is independent of the shock speed, , the upstream field obliquity angle (to the shock normal, which is in the -direction throughout this paper), and any details of the scattering process. This canonical result has propelled the popularity of shock acceleration as a key element of paradigms promoted over the last four decades for the generation of cosmic rays.

In contrast, it is widely understood that because plasma anisotropy is prevalent in relativistic shocks, when , the index of the power-law distribution is a function of the flow speed , the field obliquity angle , and the nature of the scattering. Test-particle acceleration in parallel () relativistic shocks evinces the essential property that a so-called “universal” spectral index exists in the two limits of and small angle scattering, i.e., , for a shock compression ratio of . This was showcased in the seminal work Kirk et al. (2000), which employed semi-analytic methods to solve the diffusion-convection equation, and was also generated in the Monte Carlo analyses of Bednarz & Ostrowski (1998) and Baring (1999). Here is the average angle a particle’s momentum vector deviates in a scattering “event,” and is the Lorentz factor of the upstream flow in the shock rest frame. For all other parameter regimes in relativistic shocks, a wide range of departures of from this special index is observed (see Ellison & Double, 2004; Summerlin & Baring, 2012, and references therein). While this is a complication, it enables powerful spectral diagnostics on the large scale electromagnetic structure of shocks and also the MHD turbulence in their environs.

2.2 The Monte Carlo Simulational Method

The Monte Carlo technique employed here to model acceleration at blazar shocks solves the Boltzmann equation with a phenomenological scattering operator (Jones & Ellison, 1991; Ellison, Jones & Reynolds, 1990; Summerlin & Baring, 2012). The simulation space is divided into grid sections with boundaries parallel to the planar shock interface. Sections can possess different flow velocities, mean magnetic field vectors, etc., defining the MHD structure of the shock. For electron-ion shocks, the vastly different inertial scales can be easily accommodated, though the present application is for pair plasma shocks. Charges are injected far upstream and diffuse and convect in the shock neighborhood. Their flux-weighted contributions to the momentum distribution function are logged at any position . Each particle is followed until it leaves the system by either convecting sufficiently far downstream, or exceeding some prescribed maximum momentum . To enhance the speed of the code, a probability of return boundary is introduced (following Ellison, Jones & Reynolds, 1990) at several diffusion lengths downstream of the shock, beyond which the decision to retain or discard a charge subject to convection and diffusion is made statistically: see Summerlin & Baring (2012).

The details of particle transport are given in Ellison, Jones & Reynolds (1990); Ellison & Double (2004); Summerlin & Baring (2012): the simulation can model both scatterings with small angular deflections (seeding pitch angle diffusion) and large ones () using several parameters. These deflections can be viewed as a mathematical discretization of charge trajectory segments in electromagnetic turbulence: the shorter the segment, the smaller is. The scatterings are assumed to be quasi-elastic in the local fluid frame, an idealization that is usually valid because in blazar jets and other astrophysical systems the flow speed far exceeds the Alfvén speed (true for the models in Table 1), and contributions from stochastic second-order Fermi acceleration are small. In this paper, the focus will be on small angle domains, where pitch angle diffusion is realized; results for large angle scattering are illustrated in Ellison, Jones & Reynolds (1990) and Stecker, Baring, & Summerlin (2007). The simulation assumes that the complicated plasma physics of wave-particle interactions can be described by a simple scattering relation for particles, viz.


where is the particle speed in the local upstream or downstream fluid frame, is the gyroradius of a particle of charge , and is the plasma density, with a far upstream value of . Here, () is the mean free path (spatial diffusion coefficient) in the local fluid frame, parallel to the field B, with being a fundamental bound (Bohm limit) for physically meaningful diffusion. Scalings of the momentum and the mean free path for are introduced to simplify the algebra in Eq. (1). In the local fluid frame, the time, , between scatterings is coupled (Ellison, Jones & Reynolds, 1990) to the mean free path, , and the maximum scattering (i.e. momentum deflection) angle, via for particles of speed .

Scattering according to Eq. (1) is equivalent to a kinetic theory description (e.g., Forman, Jokipii & Owens, 1974; Jokipii, 1987) where the diffusion coefficients perpendicular to () and parallel to () the local field vector B are related via


The parameter then characterizes the “strength” of the scattering and the importance of cross-field diffusion: when at the Bohm diffusion limit, and particles diffuse across magnetic field lines nearly as quickly as along them. This Bohm case corresponds to extremely turbulent fields, whose fluctuations satisfy . Each scattering event is an elastic deflection of the fluid frame momentum vector p through angle , so that the number of deflections constituting , a large angle deflection scale (i.e., a turnaround through angle ), is proportional to . Values in the range for this diffusion index are inferred from observations of turbulence in disparate sites in the solar wind, and also from hybrid plasma simulations, a context that will be discussed later in Section 4. We remark here that the story that will unfold in this paper is that blazar jets may present a picture with some similarities to the solar wind environment, with higher values of in some cases that are associated with the large dynamic ranges of momenta and mean free paths in their relativistic jets.

Figure 1: Particle distributions (normalized; ) in momentum space for acceleration simulation runs in the small angle scattering limit, corresponding to strong mildly-relativistic shocks of upstream flow speed . Here the de Hoffmann-Teller (HT) frame upstream flow speed was set at , with being the upstream field obliquity to the shock normal in the HT frame. Distributions are displayed for six different forms for the momentum dependence of the diffusive mean free path , namely with and , as labelled — see Eq. (1) and associated discussion. The shock velocity compression ratio was fixed at , and the upstream temperature corresponded to a sonic Mach number of . See Fig. 10 of Summerlin & Baring (2012) for more cases. At high momenta , many of the distributions are close to the flat asymptote that is highlighted in dark green.

In relativistic shocks, the distribution functions of accelerated charges are sensitive to the choices of both the and diffusion parameters. This property is illustrated in Fig. 1 for the mildly-relativistic shock domain that is germane to our blazar study. This sensitivity is addressed at length in Summerlin & Baring (2012), where the restriction was adopted for simplicity, so that and a single value of applies for all particle momenta. In subluminal shocks, for which , i.e., those where a de Hoffman-Teller (HT) frame (de Hoffmann & Teller, 1950) can be found, distributions generally possess indices in the range . The HT frame is that where the shock is at rest and the fluid flows along B at all positions upstream and downstream. Since is the upstream field obliquity angle to the shock normal in the HT frame, and because mildly-relativistic shocks are de rigueur for structures in blazar jets, this angle can be quite modest: subluminal shocks are a real possibility — see Table 1 for the values used in our multiwavelength spectral fitting studies. In contrast, superluminal shocks that possess higher obliquities have steep distributions with when diffusion is not near the Bohm limit; in these circumstances, rapid convection downstream of the shock spawns high loss rates from the acceleration process (e.g., Begelman & Kirk, 1990; Ellison & Double, 2004; Summerlin & Baring, 2012).

Representative Distributions for Blazar Studies

Examples of both constant and momentum-dependent are depicted in Fig. 1, where the power-law tails smoothly blend into the upper end of the thermal distribution. These distribution results were normalized to unit integrated density in all cases but one, with the example being normalized to so as to aid clarity of the Figure. For much of the fairly restricted range of obliquities corresponding to , when is a constant for all momenta, the larger the value of , the smaller is . This flattening trend is clearly represented by the comparison of the Bohm case (black) and the case (purple) in the Figure. In fact, when , then indices are generally observed (Summerlin & Baring, 2012) and the distribution is extremely flat. The origin of this behavior was found to be shock drift acceleration (SDA), where charges with select gyrational phases incident upon the shock from upstream, are trapped and repeatedly reflected back upstream by the shock discontinuity. Summerlin & Baring (2012) observed that the shock drift process is quickly disrupted when the turbulence levels increase and drops below around 100 or so.

Shock drift acceleration is a phenomenon that has been well-studied in the context of non-relativistic heliospheric shocks (Jokipii, 1982; Decker & Vlahos, 1986). It has also been observed in recent PIC simulations of non-relativistic electron-proton shocks (Park et al., 2012; Park, Caprioli & Spitkovsky, 2015), albeit primarily as a source of pre-injection of protons into the Fermi process that is manifested therein on large scales that sample the turbulent magnetic structures. The asymmetry of the drift electric fields straddling the discontinuity leads to a net energization during a gyrational reflection event. Successive episodes of upstream excursions and reflection at the shock precipitate efficient acceleration and postpone convective loss downstream for many shock interaction cycles. This extensive confinement to the shock layer automatically generates flat distributions with (Summerlin & Baring, 2012), in a sense analogous to the blazar jet shear layer studies of Ostrowski (1990), and more or less commensurate with distribution indices realized in some magnetic reconnection models (e.g. Cerutti, Uzdensky & Begelman, 2012). Only when the field obliquity is high enough that the shock is superluminal, and , does the powerful convective action of the flow overwhelm reflection, and the acceleration process effectively shuts off (Summerlin & Baring, 2012; Begelman & Kirk, 1990; Ellison & Double, 2004), with the index increasing rapidly above . Such distributions are nominally too steep to accommodate spectra from Fermi-LAT blazar data (e.g., see Abdo et al., 2010a), and so it is concluded that blazar jets must possess subluminal or “marginally luminal” relativistic shocks. This is the MHD shock phase space that is explored exclusively in this paper.

These acceleration distribution results are extended here to cases where is an increasing function of momentum , as in Eq. (1), with representative examples being displayed in Fig. 1. The subluminal shock parameters used therein are mostly those in Fig. 10, left panel, of Summerlin & Baring (2012). The distributions all display the generic trait of a dominant thermal population with a power-law tail that extends as high as the geometric scale of the diffusive acceleration zone permits: this environmental parameter is addressed in Section 3. Each distribution possesses an injection efficiency from thermal into the Fermi process, i.e. for slightly supra-thermal momenta , that reflects a combination of the values of and in Eq. (1). The ultimate power-law index at high momenta is determined by large values for in all these simulation runs, i.e., realizing domains due to efficient shock drift acceleration in putatively weak turbulence.

The key new spectral signature presented here is the appearance of “flattening” breaks in the high energy tail that are manifested when begins to exceed values around ; see the , and , cases in Fig. 1. Hence, models can possess relatively inefficient injection at thermal energies, with dominant thermal components, but exhibit efficient acceleration at momenta . The properties of plasma turbulence that might generate circumstances are discussed below in Section 4, where the blazar spectral modeling that is presented in Section 3 is interpreted.

3 Multiwavelength Radiation Emission Modeling

In this section, the focus is on modeling the broadband continuum emission of the blazars Mrk 501, BL Lacertae and AO 0235+164 using specific thermal plus non-thermal electron/pair distributions generated by the Monte Carlo technique. There are quite a few shock parameters that can be varied, so to maximize our insights, we keep the MHD Rankine-Hugoniot structure of the shock identical for all runs, and vary mostly the magnitude and direction of the magnetic field, and the diffusion parameters and highlighted in Eq. (1). Accordingly, flow speeds, compression ratios and upstream plasma temperatures are fixed so as to reduce the myriad of model possibilities.

The objective is to assess how the global character of the radio to X-ray to gamma-ray data can provide insights into the nature of the accelerator and its plasma environment. To facilitate this goal, we concentrate solely on leptonic models (Maraschi et al., 1992; Mastichiadis & Kirk, 1997; Chiang & Böttcher, 2002; Böttcher et al., 2013), those with synchrotron emission predominantly at frequencies below Hz, and with inverse Compton emission dominating the gamma-ray signal. Hadronic models are also of great interest, but then introduce more parameters, and so don’t tend to provide constraints that are as restrictive as those explored here. We remark that the critical synchrotron constraints on the values of and will apply to either leptonic or hadronic scenarios for the gamma-ray signals.

3.1 Radiation Code and Geometry Essentials

To evaluate the light emission from the complete thermal+non-thermal electron distributions generated by the Monte Carlo simulation described in the previous subsection, we adopt radiation modules from the one zone blazar radiation transfer code of Böttcher et al. (2013) and Böttcher & Chiang (2002). This code treats synchrotron emission, synchrotron self-absorption at low radio frequencies, and also inverse Compton scattering of both synchrotron and external (nominally of disk/torus origin) seed photons. Bremsstrahlung is modeled in the code, but is usually insignificant in blazars: it is so for all our case study blazars, Mrk 501, BL Lac and AO 0235+164. The emission components are computed for isotropic distributions in the jet frame and then blue-shifted and Doppler-boosted to the observer frame using the bulk Lorentz factor of the jet, and the observing angle relative to its axis.

Note that the distributions generated in the shock acceleration simulations are not isotropic in all reference frames. The dominant, compressed field is located downstream of a shock, and this is the zone where one anticipates the major contribution to both synchrotron and SSC emission originates. Since the diffusion process approximately isotropizes the electron distribution in the downstream fluid frame, this is technically the reference frame for which the radiation modules apply. Therefore the applicable boost would be by a Lorentz factor that is mildly-relativistic relative to the mean jet motion represented by . For simplicity, here we ignore this subtle distinction to reduce the number of model parameters, and set .

To complete the radiation modeling, the gamma-ray portion of the spectrum was then modified to treat line-of-sight attenuation of photons in propagation from source to Earth. This absorption derives from pair opacity due to the intervening extragalactic background light (EBL) in the optical and infra-red; this is of stellar and dust origin in galaxies. Such opacity due to the EBL has been extensively studied over the past two decades, and many models for its space density and its impact on VHE gamma-ray spectra of blazars and gamma-ray bursts exist. Here, we adopt the attenuation correction model of Finke, Razzaque & Dermer (2010). Note that the radiation models also treat pair opacity internal to the blazar jets, which can lead to attenuation of -rays by IR and optical synchrotron fields. In practice, the jet Lorentz factors are chosen sufficiently high for all three blazars studied that the jets are internally transparent to this pair opacity.

Opting for a one-zone radiative model here is motivated by simplicity. In reality, one expects multi-zone structure to the emission region, encompassing gradients in the field magnitude, electron density and other system parameters. Such complexity is the natural inference from high-resolution imaging of jets in radio, optical and X-ray wavebands. Multi-zone constructions effectively introduce more modeling variables, not necessarily precipitating a gain in insight. The most noticable change introduced by upgrading to a multi-zone scheme is the broadening of the peaks in both the synchrotron and inverse Compton components: this muting of the spectral curvature of these turnovers is a consequence of distributing the cooling rates for each of the processes over a broader range of values.

The spatial structure conceived in a multi-zone model is a sequence of segments subdividing the jet where acceleration and emission are produced in each segment. For our one-zone construction, we consider just one such segment. The acceleration zone is the most confined portion of this region, being a planar slab of thickness with a planar shock dividing the zone roughly in halves. In terms of the acceleration, since the flow is relativistic, most of the diffusive transport occurs in the downstream region, while the shock drift contribution straddles the shock and encroaches primarily into the upstream region.

On larger scales, there is a radiation volume whose linear dimension perpendicular to the shock layer is defined by the approximate equality between the radiative cooling time (synchrotron and/or inverse Compton, whichever is the shortest), and the jet fluid convection time. Since the field is compressed by the shock, the downstream region is where the synchrotron emission is most prevalent, and if the gamma-rays are generated by the SSC mechanism, then here also is where the inverse Compton signal is the strongest. MHD turbulence will permeate the entire radiating volume, but is most likely more concentrated near the shock, i.e. in the acceleration zone. The reason for this is that turbulence is naturally generated in the shock layer by dissipation of the ballistic kinetic energy of the upstream flow as it decelerates, and this turbulence should abate with distance from the shock due to damping. Such is the nature of interplanetary shocks and their environs embedded in the solar wind. A schematic of this model geometry is given in Figure 2, wherein the depicted turbulent fields are actually a two-dimensional projection of Alfvén-like fluctuations simulated using a Kolmogorov power spectrum of modest inertial range. An exponential damping of the fluctuations is imposed on scales for the purposes of illustration.

Figure 2: Schematic of the blazar model geometry under consideration, consisting of a region proximate to the shock that is the acceleration zone within a slab of approximate thickness , which is embedded in a much larger radiation zone of size (not to scale). The MHD background magnetic field structure of the shock is indicated by the blue straight line segments. Superposed on this is a turbulent field, signified by the red field line projections that are computed from a Kolmogorov power spectrum of finite inertial range spanning around a decade in wavenumber . This perturbation is exponentially damped on a scale of in this depiction so as to highlight confinement of turbulence to the shock layer.

In the absence of information on the physical scale of the acceleration zone, the Monte Carlo simulations do not describe the high-energy cut-off of the electron distribution. We therefore extend the initial electron spectra like those in Fig. 1 out to a maximum energy . Such an extrapolation is acceptable because the asymptotic power-law has been realized in the Monte Carlo simulation results with the displayed dynamic range in momenta. This maximum energy is constrained in two ways:

  1. particles will not be accelerated beyond an energy for which the (synchrotron plus inverse Compton) radiative cooling time scale, , where is the sum of the jet frame energy densities in the magnetic field and the photon field, is shorter than the acceleration time scale, , which will be discussed in greater detail below;

  2. particles will not continue to be accelerated once their energy establishes a diffusive mean free path in Eq. (1) that exceeds the size of the acceleration zone. Such energetic charges are assumed to escape upstream or downstream from the shock environs.

Here for energetic leptons. With determined as the smaller of the limiting values from the two constraints above, the high-energy portion of the electron distribution injected at the acceleration site can be written as for , where is the high-energy index of the simulated Monte Carlo distribution of accelerated particles.

The radiation module then uses as input an injected distribution, of thermal particles plus charges accelerated up to that are redistributed over the entire radiation zone, from which they escape on a time scale . This escape parameter sets the normalization of the effective equilibrium density of electrons in concert with the jet luminosity. If when , a break in the electron spectra is expected at a Lorentz factor


where the radiative cooling time equals the escape timescale. At this , the distribution steepens by an index increment . As a reminder, is the total field plus radiation energy density in the co-moving frame of the jet, and this result applies for inverse Compton scattering in the Thomson limit, a fairly representative situation. This cooling break at is in addition to acceleration flattening breaks illustrated in Fig. 1, and generates a corresponding break in the photon spectrum by an index of , the well-known signature of the transition from slow to fast cooling at high Lorentz factors.

To couple the emission signal to the acceleration and cooling information, the partially-cooled electron distributions are normalized to a total kinetic luminosity of electrons in the blazar’s two jets,


which is a free input parameter for fitting purposes. The factor of 2 accounts for both the bright approaching jet that we observe, and the vastly fainter receding one. As most of acceleration simulations illustrated in this paper produce distributions that are extremely flat, approaching (i.e., ), a substantial and sometimes dominant contribution to comes from the non-thermal electrons up to the cooling break . This is the case for the distributions exhibited in Fig. 5 below. Observe that since changes in the electron spectrum and density will lead to changes in the co-moving synchrotron photon field, our code determines the equilibrium electron distribution through an iterative process. It starts by considering only the magnetic and external radiation energy densities to determine a first-order equilibrium distribution. Then is used to evaluate the co-moving synchrotron photon field, which is added to the energy densities to re-determine the electron cooling rates and to re-evaluate the break and maximum electron energies. The process is repeated until convergence is achieved, in just a few iterations.

The magnetic field is specified by means of a magnetic partition fraction , where


is the power carried by the magnetic field along the jet, partly in the form of Poynting Flux. Two powers of appear here, defining the energy boost and the time dilation factors. From the perspective of the radiation modules, the field is effectively assumed to be tangled in the co-moving jet frame so that isotropic emissivities for both synchrotron and inverse Compton processes are employed. In reality, the acceleration considerations capture both turbulent and quasi-laminar fields at different relative levels on different spatial scales. Furthermore, the shocks are moving at mildly-relativistic speeds in the jet frame, thereby engendering “beaming” anisotropies in the electron distribution; these will be neglected for the radiation modeling since a number of the other parameters have a more important impact on the results. The radiative output from the equilibrium electron distribution is evaluated using the radiation transfer modules of Böttcher et al. (2013).

Two parameters derived from the spectroscopic models define the characteristic rates for acceleration processes. These are the non-relativistic electron gyrofrequency, , and the electron plasma frequency :


The gyrofrequency specifies the rate of acceleration in gyroresonant plasma mechanisms such as shock drift acceleration, while the plasma frequency governs the rate of energization in electrodynamic mechanisms such as magnetic reconnection. The ratio of these two frequencies is captured in the non-relativistic plasma magnetization parameter


familiar in kinetic plasma and MHD simulation studies. The electron density is obtained from the kinetic luminosity in Eq. (4). We highlight these parameters here as they will facilitate the interpretative elements later on.

3.2 Generic Features of Synchrotron plus Inverse Compton Blazar One-Zone Models

Before progressing to the blazar modeling, it is instructive to provide an outline of the central elements of the multiwavelength fitting diagnostics of the plasma environment in blazar jets. The key constraints derived are representative values of at the low injection momenta, , and at the maximum momentum in the acceleration zone. These two mean free path parameters couple via the index of the diffusion law in Eq. (1).

Foremost in our considerations is the value of , which far exceeds unity for our chosen blazars. The origin of this property is that large are required to fit the frequency of the synchrotron peak just below the spectral turnover of this emission component. For strong cooling in the acceleration zone, the turnover is created when the cooling rate


in the jet frame is approximately equal to the electron acceleration rate for the diffusive Fermi process


Note that this acceleration rate is an approximate scale, and more precise forms that isolate upstream and downstream contributions, and field obliquity influences can be found in reviews like Drury (1983) and other papers such as Jokipii (1987). Now introduce a cooling parameter that represents the fractional contribution of the synchrotron process to the electron cooling in the jet frame. The contribution to is determined from the comoving radiation field (radio-to-gamma-ray) from the radiation codes as outlined above. This synchrotron parameter can be approximately expressed in terms the ratio of optical/X-ray “peak luminosities” to the pseudo-bolometric luminosity: . For Mrk 501, this has a value of about , as can be inferred from the broadband spectrum in Fig. 4. In contrast, for the synchrotron-dominated blazar BL Lacertae, , which is somewhat higher than the values inferred for Mrk 421 and Mrk 501. For the flare state of AO 0235+164, the strong gamma-ray signal sets .

The acceleration/cooling equilibrium then establishes a maximum electron Lorentz factor at


assuming that in a mildly relativistic shock with . Here is the classical electron radius. If , then synchrotron emission dominates the cooling of particles, and when , corresponding to being independent of electron momentum, the well-known result that follows. Since for a fine structure constant and Gauss as the Schwinger field, the scale of in the range is established for blazars with fields Gauss in the jet frame, depending on the value of . This is readily seen by recasting Eq. (10) in the form


for in units of Gauss, and


being of the order of unity for blazar shocks. Increasing above unity lowers because the acceleration becomes less efficient (slower) in competing with cooling.

The synchrotron turnover/cutoff energy corresponding to is


The blueshift factor due to Doppler beaming has been included, where is the observer’s viewing angle with respect to the bulk velocity of the jet in the emission region. Also, is the redshift of the distant blazar, so that represents an observer’s measurement of the energy. For the Bohm limiting case of , imposing both and , and setting , the and evaluation yields MeV, independent of the magnetic field strength. This widely-known result was highlighted in De Jager et al. (1996) for considerations of -ray emission at relativistic pulsar wind nebular shocks, but was derived much earlier by Guilbert, Fabian & Rees (1983) in the context of AGN. This special result holds provided that the acceleration process is gyroresonant, which is the prevailing paradigm for both non-relativistic and relativistic shocks.

For Mrk 501 and our other case study blazars, becomes weakly dependent on , and the turnover energy quickly drops below the scale. For representative jet fields Gauss, one quickly estimates using Eq. (13) that for the synchrotron peak energy to appear in the optical, one requires values , and often also values . This is a central feature of our case studies below: in order to move the synchrotron peak into either the X-ray or optical band, large values of are required. In such circumstances, the acceleration process may not be gyroresonant, a realization that we will discuss in Sec. 4. The compelling need for in the shock acceleration interpretation of blazar activation was first emphasized by Inoue & Takahara (1996), who chose a momentum-independent when exploring multi-wavelength modeling of Mrk 421 spectra. In our more versatile construct here, is not fixed to unity. To realize such large at , without forcing to similar large values that would suppress efficient injection of charges from the thermal population into the Fermi mechanism (Summerlin & Baring, 2012), values are necessitated.

Synchrotron self-Compton models possess an additional constraint on system parameters via the observer’s frame energy of the SSC peak in the hard gamma-rays. This is just positioned according to the energy gained by synchrotron photons of around the turnover energy subjected to inverse Compton scattering. For the Thomson limit, this energy enhancement ratio is , so that


In conjunction with Eq. (13), this restricts the values of and . This constraint on is not applicable for external inverse Compton models where the target photon field is not the synchrotron population, but perhaps originating from a proximate disk; this will actually be the case for our blazars BL Lac and AO 0235+164. It is also modified somewhat when the upscattering samples the Klein-Nishina regime, a domain that is barely encroached upon in most of our models, since they satisfy .

This ensemble of diagnostics for our model parameters is summarized in Fig. 3. In it, an array of observational data from the multiwavelength campaign on Mrk 421 of January – May 2009 is depicted: this serves as a selection from Fig. 8 in Abdo et al. (2011b). Mrk 421 was chosen for the Figure as it will not be studied here, though we note that its broadband spectrum is qualitatively similar to that of our case study blazar Mrk 501. Schematic model synchrotron and SSC spectra are exhibited in Fig. 3, though specifically as would be computed in the jet rest frame. This then permits identification of the Doppler factor via measurement of the frequency ratio between peaks in the observed and modeled spectra for each of the emission components. Additionally, the ratio of observed to theoretical peak heights measures the flux enhancement between model and data, noting that this element of the schematic is not to scale. The separation of the component peaks (observed or modeled) defines the inverse Compton scattering enhancement factor (Thomson limit), and so constrains the value of approximately according to Eq. (14). For the illustrated case, and the synchrotron peak energy is around , so that the Compton scattering near the TeV peak is marginally in the Klein-Nishina regime. These are standard paths to computing spectral fitting parameters for SSC models of blazars. What is more unique to the present study is the determination of the diffusion parameter by comparing the synchrotron model peak energy with the fundamental bound , as depicted in Fig. 3. This is essentially inverting Eq. (13) to yield a combination of and that is captured in Eq. (16) below. These two diffusion parameters will form a central focus for our case studies below.


Figure 3: Generic construction of multiwavelength synchrotron self-Compton (SSC) spectral modeling for blazars in the representation. The illustration includes a selection of data for Mrk 421 published in Fig. 8 of Abdo et al. (2011b) for a campaign during the first half of 2009. This includes the ochre band signifying various radio detections, purple squares for optical data from different facilities, yellow-green (Swift-XRT) and green circles (Swift-BAT) for X-rays, red dots for Fermi-LAT measurements, and black triangles for MAGIC TeV-band data. The solid curves are the schematic model spectra applicable to the comoving jet frame, with radio-to-optical/UV (orange) curve constituting the synchrotron component, and the X-ray-to--ray (blue) curves denoting the SSC spectrum. The frequency offsets between the model peak and the data peak for these two components (marked by solid and dashed green vertical lines, respectively) approximately define the value of the jet Doppler factor . The separation of the peaks of the synchrotron and SSC components is by a factor of the order of . In the jet frame, the synchrotron model peak lies at an energy that is a factor below the fundamental bound of (see text).

Observe the appearance of a thermal SSC component in the X-rays due to bulk Comptonization of shock-heated thermal electrons; this actually emerges smoothly above the non-thermal inverse Compton contribution at flux levels below those of the plot scale. A similar thermal synchrotron component is not exhibited (below 1 GHz) in the Figure, since it would be suppressed by synchrotron self-absorption in the radio. Each component also exhibits a cooling break, at Hz and Hz, respectively, where the photon spectrum steepens by an index of at higher frequencies. The approximate value of can be deduced since roughly represents the ratio of the energy of the break to the thermal peak energy in the SSC component. This follows from the presumption that the shock heated electrons are at most mildly-relativistic (true for mildly-relativistic blazar shocks), and one would conclude that from the depiction in Fig. 3. The precise value of in Eq. (3) yields a measure of the relative sizes of the acceleration and cooling zones.

It is evident that constraints imposed by the energies of the synchrotron and SSC peaks on jet environmental parameters lead to strong couplings between them. The most obvious of these is the calibration of the magnetic field strength using the ratio of the SSC and synchrotron peak frequencies. Combining Eqs. (11) and (14), one arrives at the coupling between and :


for . Evidently, for synchrotron peaks in the X-rays (HBLs), fields Gauss require , whereas, for blazars with synchrotron peaking in the optical, higher values of are needed to establish similar jet fields. The ratio of the SSC to synchrotron peak luminosities also provides a modest constraint on via the parameter.

The second diffusion/acceleration parameter, , can be introduced by inverting Eq. (13), again for the case . Set to be the synchrotron peak frequency in the comoving frame in units of Hz, and as the SSC peak frequency in the jet frame in units of Hz. Eliminating using Eq. (15), one can then ascertain the fiducial relationship between and for representative observed blazar synchrotron and SSC peak energies:


With , this clearly suggests that cases are needed for synchrotron peaks to appear in the optical, and concomitantly . The physical implications of such diffusion parameters will be addressed in Section 4. Lower values of and can be entertained if the synchrotron peak appears in the X-rays, as it does for Mrk 501.

3.3 Case Study: Mrk 501

An obvious first choice for study is the well-monitored BL Lac blazar Markarian 501, located at a redshift of . Multiwavelength campaigns are readily available for this source (e.g., Catanese et al., 1997; Albert et al., 2007), thereby eliminating problems associated with non-contemporaneous data in different wavebands. The more recent campaign conducted during the March–July 2009 epoch serves amply to illustrate the advances in understanding offered by combining shock acceleration simulation results with radiation emission codes. The broadband spectral character of Mrk 501 differs only modestly from that depicted in Fig 3 for its sister blazar Mrk 421, and is also fairly similar to the SED of more distant blazars such as PKS 2155-304 (at ; see Aharonian et al., 2009), yet it is still suitable to an SSC model construction with an added optical component.

An extensive summary of the mid-March–early August, 2009 campaign involving the VERITAS and MAGIC telescopes above 100 GeV, Fermi-LAT in the 20MeV to 300 GeV energy window, RXTE-PCA (2–60 keV) and Swift XRT+BAT (0.3-150 keV) in the X-rays, combined with UV (Swift-UVOT), optical and radio measurements, is provided in Abdo et al. (2011a). For most of this epoch, the source was in a low/moderate state, for which broadband spectral data are depicted in Fig. 4. During this campaign, in the short window May 1–5, 2009, Mrk 501 underwent a sizable flare in gamma-rays, evincing a high state at around 5 times the Crab nebula flux in the VERITAS/MAGIC/Fermi-LAT bands. This was an enhancement of a factor of 3-5 over the quieter portions of the campaign, and received considerable attention in Abdo et al. (2011a) and also Aliu et al. (2016). The flare exhibited a spectrum that was much harder above Hz than for the extended state. To aid clarity, the VERITAS data were not depicted in Fig. 4, being similar to the MAGIC points above Hz. Moreover, the multi-wavelength data represent the long-term average over the 4.5 month interval, and for radio, optical, UV and X-rays, were taken from the spectrum in Fig. 8 of Abdo et al. (2011a). The gamma-ray points, for Fermi-LAT and MAGIC, omit an interval of 30 days that brackets the strong flaring episode, and are taken from Fig. 9 of Abdo et al. (2011a). Note that the units for the -axis in this and subsequent multi-wavelength spectral figures employ the Jansky-Hertz choice familiar to radio astronomers, and can be converted to the form in Fig. 3 that is often used by high energy astrophysicists via 1 Jy Hz erg cm sec.


Figure 4: Multiwavelength spectra (points), together with model fits as described in the text, for the extended March-July 2009 monitoring of the blazar Markarian 501. The campaign data are taken from Figures 8 and 9 of Abdo et al. (2011), and constitute a “low” state with gamma-ray data bracketing the flare of early May, 2009 being omitted (see text). The gamma-ray detections and upper limits are from Fermi-LAT (purple triangles) and MAGIC (black circles), and the X-ray points are Swift-BAT (muddy orange squares), Swift-XRT (light blue triangles) and RXTE-PCA data (open red circles). Optical, UV and radio measurements are detailed in Abdo et al. (2011); since a variety of radio fluxes were recorded for various regions much larger than the compact X-ray/-ray zones, they are marked by a representative band, and were not fit. The broadband models consist of a synchrotron component (dashed green curve) up to the X-ray band, an SSC component in the X-rays and gamma-rays (dashed red curve), and a separate thermal host galaxy emission component (dotted black). The full jet model spectra for the extended low state includes a correction (Finke, Razzaque & Dermer, 2010) for absorption by the extragalactic background light, and combined with the data, constrain the diffusion parameters to and (see text).

In the SSC interpretation, the extended state is very slightly synchrotron-dominated, contrasting the May 1–5, 2009 high state, which is marginally inverse Compton-dominated; throughout, the synchrotron peak is positioned in the X-ray band. The low-variability UV/optical data is not attributed exclusively to the jet, but mostly to the host galaxy: at least 2/3 of the flux is believed (e.g., Nilsson et al., 2007) to originate in the host, and is herein modeled via a separate thermal component. The radio data were from various (mostly single-dish) facilities (see the list in Table 1 of Abdo et al., 2011a) that generally do not resolve the jet, and accordingly provide upper bounds to the signal expected from the highly-variable -ray emission zone. For this reason, they are depicted in Fig. 4 as a swath, positioned at the fluxes measured at the various radio frequencies.

Using the radiation modeling protocols and coding described in Section 3, multiwavelength spectra from radio to VHE gamma-rays were generated for a suite of particle distributions generated by the Monte Carlo shock acceleration simulation. The turbulence parameters and were adjusted for the various simulation runs to hone in on a candidate “best case,” the spectral results for which are exhibited in Fig. 4. The spectra so generated are dominated by two components, synchrotron emission in the radio-to-X-ray band, and a synchrotron self-Compton signal in the gamma-ray range; the observer frame computations of these components at the source are isolated in this Figure as dashed curves. The overall spectrum at Earth is subjected to attenuation from the line-of-sight EBL, as described above, and is depicted as the solid M/W curve in Fig. 4 that is the nominal fit to the data; this attenuation is omitted in the SSC component depiction so as to illustrate its magnitude.

Simulation and radiation model input parameters, and also those derived from the input ones, are listed in Table 1; some of them can be compared with model fitting parameters listed in Table 2 of Abdo et al. (2011a). Since in this broadband spectrum the thermal and non-thermal components are obviously distinct, the detailed shape of the host galaxy optical spectrum is irrelevant for our modeling goals. Therefore, no attempt was made to perform a detailed fit to the optical (e.g., with a spectral template for an elliptical galaxy); this component was simply represented with a separate blackbody spectrum for illustrative purposes.

Parameter Mrk 501 BL Lacertae AO 0235+164
Name Symbol/units Extended State Extended State High State
Jet+Source Parameters
Redshift 0.034 0.069 0.94
Jet Lorentz factor 25 15 35
Observing angle
Emission region size [cm]
injection luminosity4 [erg/sec]
Escape time scale5 100 5 300
Magnetic partition 0.06
Dusty torus luminosity [erg/sec]
Dusty torus temperature [K]
Shock Parameters
Upstream speed 0.71 0.71 0.71
Field obliquity
Upstream temperature (K)
Compression ratio 3.71 3.71 3.71
Injection mean free path 100 20 225
Diffusion index 1.5 3 3
Derived Parameters
Magnetic field [Gauss]
Field luminosity [erg/sec]
Synchrotron partition6 0.6 0.83 0.25
Cyclotron frequency
plasma frequency
Shock magnetization7
Alfvén speed c c c
Cooling break 2.0
Maximum Lorentz factor
Maximum mean free path8 [cm]
Synchrotron cutoff [Hz]
SSC cutoff frequency [Hz]
Table 1: Input and Derived Parameters for Blazar Models

The principal constraints derived from the Mrk 501 multiwavelength modeling are that the mean free path parameter is modest for mildly-relativistic electrons injected into the shock acceleration process, and that at the highest energies, in order to place the synchrotron turnover (i.e., peak) in the hard X-ray band. The values we obtained for the extended state were and (i.e. ). These parameters generate moderately efficient injection from thermal energies into the Fermi mechanism (see Fig. 1), and mean free paths that possess a stronger momentum dependence than Bohm diffusion (), suggesting interactions with weaker turbulence for more energetic particles diffusing on larger spatial scales. The interpretive elements of this parametric fit will be embellished upon in Section 4.

It should be borne in mind that this fitting “solution” is representative of the appropriate parameter space, but is by no means a unique choice. It can be generally stated that there is a tolerance of about in , and a tolerance of a factor of in permitting M/W spectral fits of similar quality. The uncertainties in the gamma-ray fluxes, particularly in the Fermi-LAT band, preclude greater precision for modeling. Moreover, in this analysis, the shock MHD parameters were kept fixed at the values chosen for Figure 10 of Summerlin & Baring (2012), namely a shock speed of , a velocity compression ratio of , a field obliquity of , and a low sonic Mach number of around 4. Adjusting the MHD structure of the shock would introduce changes to the optimal choices for the shock layer turbulence parameters and in a data fitting protocol, so the ones for spectra exhibited in the Figure serve as general guides, and should not be presumed sacrosanct.

Likewise, the choice of the observational epoch should always be borne in mind: Mrk 501 is highly variable above the optical band, and so sets of fitting parameters will change when different epochs are modeled. This is evident for strong flaring episodes such as the May 1–5, 2009 event not treated here. Yet it also applies for other observational campaigns, such as the recent NuSTAR-led multi-wavelength monitoring spanning the period April 1 – August 10, 2013 (Furniss et al., 2015). The X-ray spectrum therein is often harder and more luminous than that depicted here, though the gamma-ray spectrum (MAGIC+VERITAS+Fermi-LAT) is more or less commensurate with that in Fig 4 — see Fig. 9 of Furniss et al. (2015) for details. To appreciate how extreme Mrk 501 can become during flares, one need look no further than the depiction of the April 1997 flare in Fig. 6 of Acciari et al. (2011), wherein the synchrotron peak moves to above Hz (i.e. nearly two decades higher than the listing in Table 1), and its flux level is slightly over a decade higher than that illustrated in Fig. 4. Modeling this extreme HBL behavior would lower the inferred value of by a factor of 30–50. Yet for this flare, and also for the NuSTAR M/W campaign, the principal conclusion of a strong momentum dependence for the diffusion coefficient and large departures from the Bohm limit at would be upheld, and is anticipated to apply to a broad selection of observing campaigns.

The electron distribution corresponding to the spectrum in Fig. 4 is illustrated in Figure 5, in a manner analogous to the representation: distributions give an approximate scaling of the emission power of synchrotron and SSC signals from electrons of a given Lorentz factor . The thermal Maxwellians are fairly prominent, yet they couple efficiently and smoothly into the non-thermal accelerated population. The exponential turnovers near the maximum energies are precipitated by efficient radiative cooling in the compact acceleration zone, depicted in Fig. 2. This shuts off the acceleration, and the electron distribution evolves through continued cooling in the larger radiation zone. Cooling ceases when the radiation cooling length becomes comparable to the size of the radiation region, which is specified by the parameter in Table 1, and this introduces the break at , the approximate value of which satisfies Eq. (3). For the extended state of Mrk 501, since the synchrotron component slightly dominates the SSC one in terms of the energy flux, the inverse Compton mechanism contributes the minority of the cooling; for any considerations of the high state, the situation would be reversed. For each process, efficient cooling steepens the electron distribution by an index of unity above the cooling break, and this induces the well-known steepening by an index of in the photon spectrum. Such a break is seen at around Hz for the synchrotron component in Fig. 4, and is barely discernible at around Hz for the SSC contribution therein.

We note in passing that the injected shock acceleration distributions that led to those depicted for the larger cooling zone in this Figure have most of their energy allocated to the particles near . In such cases, the acceleration simulations need to be upgraded to account for non-linear modifications to the shock MHD structure imposed by the pressure of the energetic non-thermal particles. This feedback phenomenon is well understood in non-relativistic shock theory (e.g. Ellison & Eichler, 1984; Ellison, Baring & Jones, 1996; Baring et al., 1999), and motivates future refinements to blazar studies of the genre herein, using extensions of this feedback theory to relativistic shocks along the lines of Ellison, Warren & Bykov (2013).


Figure 5: Complete thermal plus non-thermal electron distributions employed in the multiwavelength modeling results for Mrk 501 depicted in Fig. 4, and BL Lacertae as shown in Fig. 6, as functions of the four-velocity or dimensionless electron momentum . These distributions, applicable to the jet frame, are derived from the Monte Carlo simulation results that are injected into the radiating volume of the blazar jet, and are modified by cooling losses above the cooling break energy , and eventually by the suppression of acceleration by rapid radiative cooling at , thereby precipitating the exponential tails. Given that is the distribution of electrons differential in Lorentz factor, this representation highlights the fact that most of the synchrotron and SSC radiative power is generated by electrons with the Lorentz factors near and above . The Mrk 501 case (blue histogram) had and , while the BL Lac one (red) had and . The comparison differential distribution exhibited in green is the expected power law for high cases (see Section 2).

Beholding the multiwavelength spectrum, at first sight, the character of the Mrk 501 SSC model fits here resemble those in numerous expositions on this source, including Acciari et al. (2011) and Furniss et al. (2015). However there is a key development here beyond phenomenological electron distributions, namely broken power laws that are truncated at a minimum electron Lorentz factor, which are usually invoked in other studies of Mrk 501 (e.g., see Aleksić et al., 2015) and various other blazars. These forms are unphysical from the standpoint of shock acceleration theory. Here, we have complete thermal-plus-non-thermal distributions at our disposal, and this serves to advantageously define the thermal plasma density normalization relative to that of the accelerated electron population (see Fig. 5). The plasma density connects to determination of the shock structure, and the number density of accelerated electrons benchmarks the radiative output and also the jet kinetic energy via Eq. 4. In addition, knowing the total number density fairly precisely permits the determination of the electron plasma frequency , and the plasma magnetization parameter (not possible for phenomenological power-law distributions). Inspection of Table 1 reveals that the magnetization is for our Mrk 501 fitting, a result that will be interpreted in Section 4 as implying similar speeds of acceleration for reconnection and for the diffusive Fermi mechanism in the Bohm limit. Also listed in Table 1 are determinations of the hydrogenic Alfvén speed, in all cases virtually non-relativistic and substantially inferior to the shock speed. These rise to relativistic speeds with the addition of a dominant pair component to the jet. Accordingly, here we forge an interconnection between photon emission and underlying shock plasma properties in blazar jets, in a substantial advance beyond previous works.

From the spectroscopic point of view, it is apparent from Fig. 4 that the thermal portion of the electron distribution functions played no role in constraining parameter space. The mildly-relativistic electrons in their thermal component are inefficient generators of both synchrotron and SSC radiation. Thermal synchrotron would appear in low frequency radio windows and is heavily self-absorbed by the inverse synchrotron process, which is treated in our radiation codes: the signatures of such attenuation generally appear at frequencies below Hz when Gauss. Thermal SSC appears in the optical-X-ray window, and is swamped by the strong X-ray synchrotron emission for Mrk 501. Both contributions are off-scale in the Figure. Hence, in the case of Mrk 501, spectroscopy associated with the thermal jet electrons is irrelevant. This contrasts the interesting case of AO 0235+164, our final object for study.

3.4 Case Study: BL Lacertae

To offer a picture distinct from that of Mrk 501, the focus turns to BL Lacertae, the prototype BL Lac object, which has a much softer synchrotron component that peaks in the optical (i.e., an LBL) and dominates the hard X-ray and gamma-ray signals. It is also a nearby source, at , and various multiwavelength campaigns have been staged to help elucidate its character. A survey undertaken in the early days of the Fermi mission is the focal point here. This was for the extended period of August – October 2008, and while the source was variable and possessed episodes of relatively enhanced radiative output, its broadband emission was not that of a characteristically high state.

The data compilation for this three month epoch is detailed in the Fermi-LAT collaboration paper on a multitude of so-called LBAS blazars (Abdo et al., 2010c). The purpose of this compendium was a first quick-look survey of the general character of bright blazar spectra in the Fermi era. Therefore, for BL Lacertae, the mutliwavelength data presented in that paper are depicted in Fig. 6, serving as our modeling benchmark. The data in color are from Fermi-LAT in the 100 MeV to 30 GeV energy window (blue), Swift XRT (0.3-10 keV, green) and Swift BAT (15-200 keV, purple) in the X-rays, Swift-UVOT in the near ultra-violet, radio measurements from Effelsberg, Owens Valley and RATAN, and data from other optical, infra-red and radio facilities. Since not all were precisely contemporaneous with Fermi-LAT observations and each other, modest mismatches between model fits and data should not be over-interpreted. Synchrotron emission from BL Lac is bright enough to outshine its host galaxy, and so no galaxy contribution is discernible in the optical, unlike Mrk 501. Fig. 6 also presents archival MAGIC data (Abdo et al., 2010c) in grey that were not used in the fitting protocol, and serve as a guide for the variable character of gamma-ray signals from BL Lac. Fig. 23 of Abdo et al. (2010c) can be consulted for a more extensive presentation of archival spectroscopy for BL Lac.

The radiation modeling protocols were implemented as described above, using a similar suite of particle distributions generated by the Monte Carlo shock acceleration simulation. Radiation model and simulation input parameters, and also those derived from the input ones, are listed in Table 1. The shock MHD parameters were identical to the Mrk 501 study. The turbulence parameters and for the simulations were again adjusted to optimize the fit, for which spectral results are presented in Fig. 6. The radio-to-optical data dictate that the synchrotron component is dominant, and turns over below the X-ray band. This imposes a higher value of the magnetic field and its associated portion of the total energy budget than for the case of Mrk 501. Concomitantly, synchrotron cooling of electrons in the acceleration zone is more rampant, thereby reducing the maximum Lorentz factor to (see Fig. 5). Accordingly, the synchrotron self-Compton peak energy appears in the soft gamma-rays, below 10 MeV, as is evident in Fig. 6. This circumstance is essentially unavoidable in LBLs that are synchrotron-dominated. Hence, the Fermi-LAT data require another component to be present.

To describe such, here we employ an external Compton signal seeded by emission in the near infra-red (NIR) from a dusty torus associated with the accretion flow onto the central black hole. These constitute an isotropic thermal photon field with the luminosity and temperature as listed in Table 1. It is a simple matter to discern why this EC component emerges in the Fermi-LAT band. While the torus appears in the NIR for an observer, as an external isotropic field proximate to the jet, it is boosted in the jet frame into the far-UV. This is at considerably higher frequency than the synchrotron emission, which clearly peaks in the IR band in the jet frame. Accordingly, the EC peak is substantially bluer than that of the SSC in both the jet frame, and our observer frame, by over a factor of 100 higher in energy. Moreover, the Doppler-boosted flux enhancement of the torus seed photons in the jet frame can generate a gamma-ray flux comparable to the SSC without the seed NIR signal being discernible above the very strong, beamed, synchrotron component at those frequencies. The breadth of the EC peak is restricted by the intrinsic width of the Planck spectrum from the torus.


Figure 6: Multiwavelength spectra (points) spanning the radio, optical, X-ray and gamma-ray bands, together with model fits as described in the text, for the August–October 2008 extended Fermi-LAT observation of the blazar BL Lacertae. The data are taken from the LBAS study of Abdo et al. (2010), and are briefly detailed in the text: the colored data are approximately contemporaneous with the Fermi-LAT collection. The grey MAGIC points above Hz are archival, being shown to illustrate the spectral variability between different observational epochs; they were not accommodated in the model fit. The broadband models consist of a synchrotron component (dashed green curve) up to the optical band, an SSC component in the X-rays and gamma-rays (dashed purple curve), and an external inverse Compton emission component (dash-dot grey) seeded by IR photons from a dusty accretion torus. The full jet model spectrum (orange) includes a correction for absorption by the extragalactic background light, and combined with the data, constrain the model diffusion parameters to and (see text).

While this introduces an extra subset of model parameters, there appears no likely alternative in one-zone leptonic models. Similar broad spectral structure from X-rays to hard gamma-rays could be generated by a multi-zone blazar emission model: such was the protocol adopted by Zacharias & Wagner (2016) in the modeling of AP Librae, an LBL object whose multiwavelength spectroscopy is very similar to that of BL Lacertae. A possible discriminator between these two competing pictures is the appearance of a flat bridge in the spectrum from our protocol between the two inverse Compton components. Such a feature might not appear in a multi-zone construction, which may evince broader spectral curvature. This provides motivation for future sensitive gamma-ray detectors below the Fermi-LAT window, perhaps employing Compton telescope technology.

At the very highest energy gamma-rays, as for Mrk 501, our modeling incorporated a pair opacity correction (Finke, Razzaque & Dermer, 2010) due to the intervening EBL. However, since the EC peak in BL Lac lies at 10 GeV and well inside the LAT window, this correction is quite small, and did not influence our fitting protocol. It is anticipated that for this source, when it is in a hard state as captured in the non-contemporaneous MAGIC data, such attenuation compensations will come into play.

It is notable that BL Lacertae is singular in our ensemble of blazars in that it has measurements of the jet bulk velocity. Jorstad et al. (2005) outline VLBI imaging, polarimetry and variability studies of a number of radio jet sources, including BL Lac. Combining imaging and flux variability information at 7mm, they were able to isolate both the apparent superluminal speed of the jet and the Doppler factor. This yielded a determination of the bulk Lorentz factor of for BL Lac: see Table 11 of Jorstad et al. (2005) and associated discussion. In a later investigation of TeV flaring in BL Lac using VERITAS, Arlen et al. (2013) employed internal pair transparency arguments based on Dondi & Ghisellini (1995) to constrain the jet bulk motion to Lorentz factors . Arlen et al. (2013) argued that this determination was possibly not in conflict with the lower value obtained in the VLBI study, with the gamma-ray emitting region being distinct from the radio ones and perhaps lying closer to the central black hole, thereby invoking a decelerating jet picture. Since these gamma-ray constraints are most germane to our multiwavelength investigation here, they guided our adoption of for the models of BL Lac. This guaranteed pair transparency for the jet for gamma-rays below 1 TeV. Higher values of that permit jet pair transparency can also lead to viable fits with the same particle distributions from acceleration theory, yet require reasonable adjustments of environmental parameters such as the field strength and electron number density.

As with Mrk 501, the upshot of the M/W modeling is that a large value of is required to generate relatively efficient injection of electrons from thermal energies into the acceleration process, and inefficient acceleration at the highest Lorentz factors. Notably, since BL Lac is an LBL blazar, a value of was required to increase above and move the synchrotron peak into the optical band. The values we obtained, and (i.e. ), indicate diffusion well away from the Bohm limit at all momenta, and that the level of turbulence must decline with distance from the shock even more profoundly than for Mrk 501. As with our first case study, there is a tolerance of about in the range of , and a tolerance of a factor of in permitting M/W spectral fits of similar quality. The value of obtained for BL Lac is noticeably lower than for our two other case study blazars. Yet, it is also somewhat higher than the value of determined from fitting thermal+non-thermal distributions of protons and alpha particles that are measured in situ at interplanetary shocks in the non-relativistic solar wind (Baring et al. 1997). This suggests that mildly relativistic shocks in blazar jets are less efficient at injecting charges into the acceleration process than are traveling shocks in the solar wind. Note that the contributions to each radiation component from the thermal electron population are small enough to be off scale in the Figure. The electron distribution corresponding to the fit is depicted in Fig. 5.

An issue that has not yet been addressed concerns the variability timescale, which must be large enough to accommodate the time for acceleration up to the largest particle energies required to fit the spectral data. The acceleration rate is given by in Eq. (9), assuming a relativistic shock. Thus, the most energetic charges are accelerated on timescales of


in the rest frame of the jet. Each of these quantities on the right hand side is listed in Table 1, yielding derived values of the characteristic cumulative timescale for acceleration , also listed in the Table. Thus we infer sec for Mrk 501, i.e. around 4.5 days. Similarly, sec for Bl Lac. To connect to the observer’s frame, these are shortened by time dilation factors of for transverse variability, i.e to around sec or 4 hours for Mrk 501, and around 20 minutes for Bl Lac. Such dilation applies to fluctuations in the direction normal to the line of sight to the observer, i.e. approximately the jet axis. Longitudinal variability along the jet, such as are expected from colliding shells, reduces the acceleration timescales by another factor of , rendering them very short in the observer frame. These are then less than the variability timescales for both sources in the hard gamma-rays (e.g., see Aliu et al., 2016, for Mrk 501 data in the high state that was excluded from the epoch considered in this paper), for both extended states and flaring episodes within the observing epochs under consideration, and so no internal inconsistency arises.

3.5 Case Study: AO 0235+164

To provide another contrasting case, we turn to the LBL blazar AO 0235+164, which has a radio-optical component that turns over below the UV band. It is normally weak in X-rays, however, in flare state, it possesses a significant X-ray excess. Such enhancements have been discussed (Sikora et al., 1997; Błażejowski et al., 2000) as a possible signature of a bulk Comptonization effect (see also Sikora, Begelman & Rees, 1994), where hot thermal jet electrons upscatter an ambient, quasi-thermal radiation field that manifests itself in the infra-red, optical or UV bands. This seed field could be optical/UV radiation from the broad-line or narrow-line regions, or IR emission from warm dust in the greater AGN environment; it is distinct from the non-thermal jet synchrotron emission. One of the objectives here is to explore whether the X-ray “excess” can be attributed to the substantial thermal electron pool evinced in shock acceleration distributions such as those depicted in Fig. 1. In the fast-moving jet, this population may produce significant radiative signatures due to Comptonization of an external radiation field, taking advantage of the same Doppler boosting enhancement that was captured in our modeling of BL Lac.

The focus here is on modeling data collected during the flaring episode of September – November 2008. In particular, we employ the multi-wavelength campaign observations discussed in Ackermann et al. (2012); see also Agudo et al. (2011) for additional flare light curves and polarization data. The broadband spectrum for the high state portion of the outburst, during the interval MJD 54761-54763 (October 22-24, 2008), is depicted in Fig. 7. This appears in Fig. 7 of Ackermann et al. (2012), wherein it can be compared with a low state spectrum obtained about six weeks later that does not contain a Swift XRT detection. The list of observatories and missions participating in this campaign is extensive — it spans radio and optical, with Swift XRT and RXTE providing X-ray monitoring, and Fermi-LAT providing the principal gamma-ray observations. The LAT spectrum displays a break or possible turnover at energies around GeV, with a power-law index of about below this feature; the source is not detected in the TeV band (Errando et al., 2011).


Figure 7: Multiwavelength spectra (points) spanning the radio, optical, X-ray and gamma-ray bands, together with model fits as described in the text, for the October 2008 high state Fermi-LAT observation of the blazar AO 0235+164. The campaign data are taken from Ackermann et al. (2012), corresponding to the high state epoch during October 22-24. The gamma-ray detections and upper limits are from Fermi-LAT, while the X-ray “butterfly” block (blue) represents Swift XRT data. The broadband models consist of a synchrotron component (dashed green curve) up to the optical band, a two-order SSC contribution in the optical, X-rays and gamma-rays (dashed purple curve), and bulk Comptonization emission (EC: dotted black curve) between keV and GeV. The total model spectrum (orange) includes a very small correction for absorption by the extragalactic background light, and combined with the data, constrains the model diffusion parameters to and (see text).

The major model components and total multi-wavelength spectrum in our fit of the flare data are also exhibited in Fig. 7. The jet and disk parameters used to derive the model fit are listed in Table 1, and we remark that they are somewhat different from those in our earlier introductory exposition (Baring, Böttcher & Summerlin, 2014) on this source. Notably, the magnetic field at Gauss was sub-equipartition by about a factor of 16. To place the synchrotron turnover frequency in the optical band (Hz), the values of and have been chosen so that the maximum electron Lorentz factor in the comoving frame of the jet is . The shock obliquity was modestly higher than for the BL Lac and Mrk 501 studies, a choice that renders the thermal electrons comparatively more numerous. As with BL Lac, the model fit required and , so that the diffusive scales of energetic particles are very large. The tolerances in these parameters that permit M/W fits of similar calibre are about in the range of , and a factor of in . The moderate value of then positions the peak SSC frequency at around Hz, i.e. around 5 MeV, as is evident in Fig. 7. This is not the maximum extension of SSC component: a second-order IC image of the synchrotron continuum appears up to about 10 GeV. A cooling break is also apparent in the SSC spectrum at around a few keV. Observe that the SSC component is of insufficient luminosity to model the Swift XRT and Fermi-LAT signals — another component is needed to explain them.

These two high state detections are modeled here using bulk-Comptonization/inverse Compton scattering of an external radiation field. This field is presumed to be quasi-isotropic in the observer’s frame, so that the jet material moves at relative to these seed photons, and so can inverse Compton scatter them, i.e. increase their frequency by a factor of . This choice of lies between that adopted in some of the supporting modeling in the Fermi-LAT paper of Ackermann et al. (2012), namely , and the much higher value of indicated by structure variability in the VLBA imaging over extended epochs that was presented in Jorstad et al. (2001). We note that these values exceed the typical range of quoted in Ghisellini et al. (2014) for BL Lac objects where Compton dominance (or not) in SSC models is used to constrain  — clearly, employing external Compton scenarios will modify these estimates. The Lorentz factor of the thermal electrons in the distributions of Fig. 1 is less than about 2, thereby modifying this IC energy enhancement by at most a factor of a few. Hence, to explain the Swift XRT flaring flux, the seed background radiation needed a temperature of K. This corresponds to infra-red radiation, probably from a dusty torus. The thermal portion of the electron population boosts the ambient photons up to the X-ray range, and roughly produces the steep Swift XRT spectrum if it corresponds to the higher momentum side of the Maxwell-Boltzmann distribution, as is clearly evident in Fig. 7. This then connects to a very flat external Compton tail, generated by the broken power-law tail portions of distributions like those depicted in Fig. 1, modulo the cooling breaks discussed above. The IC tail possesses a flat spectral index of around 3/2, rising in the representation to meet the Fermi-LAT flux above 100 MeV. The large IC flux indicates very strong Compton cooling of electrons, and this generates a very low cooling break energy at .

Note that despite the high redshift () of this source relative to those of the other two blazars, the EBL absorption correction is minimal since the spectrum extends only up to around 20 GeV, i.e., similar to the situation for BL Lac. A peculiarity of this blazar is that sensitive optical/UV spectroscopy reveals the presence of an intervening damped Lyman-