Probing Acceleration and Turbulence at Relativistic Shocks in Blazar Jets
Abstract
Diffusive shock acceleration (DSA) at relativistic shocks is widely thought to be an important acceleration mechanism in various astrophysical jet sources, including radioloud active galactic nuclei such as blazars. Such acceleration can produce the nonthermal particles that emit the broadband continuum radiation that is detected from extragalactic jets. An important recent development for blazar science is the ability of FermiLAT spectroscopy to pin down the shape of the distribution of the underlying nonthermal particle population. This paper highlights how multiwavelength spectra spanning optical to Xray to gammaray 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 nonthermal 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 gammaray band above 100 MeV. We investigate selfconsistently 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 Xray 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 energydependence of the diffusion coefficient for electrons. The concomitant interpretations that turbulence levels decline with remoteness from jet shocks, and the probable significant role for nongyroresonant diffusion, are posited.
keywords:
Acceleration of particles, plasmas, shock waves, turbulence, galaxies: active, galaxies: jets, Xrays, gammarays.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 gammaray bursts (GRBs), and
persistent but highlyvariable active galactic nuclei (AGN). Both are
considered prime candidates for the production of ultrahigh 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),
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 nonthermal optical spectra. Their radio continua often possess quite flat spectra. These characteristics suggest synchrotron radiation from nonthermal 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 1510089; Marscher et al., 2010), contemporaneous with a gammaray flare seen by FermiLAT (Abdo et al., 2010b). This indicates a largescale ordering of the magnetic field in the 3C 279 jet, and spatial coincidence of the optical and gammaray emission zones. In many blazars, this putative synchrotron component extends to Xrays, where polarization measurements in the nottodistant future define a hope for further constraining the emission mechanism of blazars (e.g., Krawczynski et al., 2011).
The prevailing paradigm is that blazars’ gammaray signals are generated by inverse Compton scattering of synchrotron photons by the same relativistic electrons that emit this radiotoXray signal, socalled synchrotronselfCompton (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 nonthermal protons collide with disk and jetassociated synchrotron photons, generating pair cascades via photopion/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 infrared 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 MeV100 GeV window, afforded by the FermiLarge Area Telescope (LAT). Over the last eight years, LAT data has enabled measurements of the powerlaw index of spectra from numerous blazars (Abdo et al., 2009, 2010c). This provides a more robust measure of the underlying nonthermal 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 FermiLAT data prove critical to our interpretation of its jet environment.
The rapid flux variability seen in radio, Xray 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 gammaray telescopes. For a survey of radio and Xray 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 Xray 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 particleturbulence 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 outwardpropagating plasma may encounter with high probability.
Another possibility is that shear layers encapsulating sharp velocity gradients transverse to the net flow may precipitate Fermitype acceleration (Ostrowski, 1990) due to transport of charges straddling the shear “discontinuity.” Observational evidence supporting such transverse velocity structure comes from parsecscale limbbrightening of blazar and radio galaxy jets revealed in VLBI observations (e.g., see Giroletti et al., 2004, for Mrk 501). MHD/hydrodynamic simulations of spinesheath 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. RayleighTaylortype instabilities may develop at the spinesheath interface, and such turbulence offers another promising avenue for accelerating relativistic particles. Indeed, particleincell 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 Poyntingflux dominated outflows can be envisaged. Whether shocks or shear boundaries or reconnection zones provide the dominant energization site for nonthermal 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 gammaray 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 broadbrush 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 FermiLAT 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 powerlaw 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 Xray band. This serves as a central issue to the discourse of this paper.
In the multiwavelength 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 powerlaw index, and also the connection between the thermal bulk of the population (a hot Maxwellianlike component) and the powerlaw tail of the accelerated species. Such an approach goes beyond an elemental description of the distribution of accelerated charges, usually a powerlaw 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 onezone 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 nonthermal 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 nongyroresonant 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 multiwavelength spectra of blazars, in particular the positioning of the synchrotron peak in the optical or Xray 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 powerlaw 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 mildlyrelativistic. 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). Particleincell (PIC) plasma simulations have also been employed to study this problem via modeling the establishment of electromagnetic turbulence selfconsistently (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, nonthermal, 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 socalled 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 waveparticle interactions and diffusive acceleration well above thermal energies, but do not treat the electrodynamics of wave generation and dissipation selfconsistently. Particleincell simulations provide the greatest depth in modeling shock layer microphysics by solving Maxwell’s equations and the NewtonLorentz force equation, though their macroparticle approximation does eliminate electrodynamic information on the smallest scales. With increased box sizes, PIC codes have now realized the establishment of truly nonthermal components, but still cannot explore acceleration beyond modest nonthermal 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 nonrelativistic shock acceleration theory is that the acceleration process possesses no momentum scale, and the resulting particle distribution takes the form . For nonrelativistic shocks, since their speeds far exceed (), the upstream (downstream) flow speed component in the coordinate 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 powerlaw 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 powerlaw distribution is a function of the flow speed , the field obliquity angle , and the nature of the scattering. Testparticle acceleration in parallel () relativistic shocks evinces the essential property that a socalled “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 semianalytic methods to solve the diffusionconvection 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 electronion 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 fluxweighted 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 quasielastic 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 secondorder 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 waveparticle interactions can be described by a simple scattering relation for particles, viz.
(1) 
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
(2) 
The parameter then characterizes the “strength” of the scattering and the importance of crossfield 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.
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 mildlyrelativistic 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 HoffmanTeller (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 mildlyrelativistic 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 momentumdependent are depicted in Fig. 1, where the powerlaw 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 wellstudied in the context of nonrelativistic heliospheric shocks (Jokipii, 1982; Decker & Vlahos, 1986). It has also been observed in recent PIC simulations of nonrelativistic electronproton shocks (Park et al., 2012; Park, Caprioli & Spitkovsky, 2015), albeit primarily as a source of preinjection 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 FermiLAT 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 powerlaw 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 suprathermal momenta , that reflects a combination of the values of and in Eq. (1). The ultimate powerlaw 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 nonthermal 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 RankineHugoniot 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 Xray to gammaray 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 gammaray 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 gammaray signals.
3.1 Radiation Code and Geometry Essentials
To evaluate the light emission from the complete thermal+nonthermal 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 selfabsorption 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 blueshifted and Dopplerboosted 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 mildlyrelativistic 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 gammaray portion of the spectrum was then modified to treat lineofsight 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 infrared; 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 gammaray spectra of blazars and gammaray 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 onezone radiative model here is motivated by simplicity. In reality, one expects multizone structure to the emission region, encompassing gradients in the field magnitude, electron density and other system parameters. Such complexity is the natural inference from highresolution imaging of jets in radio, optical and Xray wavebands. Multizone constructions effectively introduce more modeling variables, not necessarily precipitating a gain in insight. The most noticable change introduced by upgrading to a multizone 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 multizone model is a sequence of segments subdividing the jet where acceleration and emission are produced in each segment. For our onezone 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 gammarays 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 twodimensional projection of Alfvénlike 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.
In the absence of information on the physical scale of the acceleration zone, the Monte Carlo simulations do not describe the highenergy cutoff 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 powerlaw has been realized in the Monte Carlo simulation results with the displayed dynamic range in momenta. This maximum energy is constrained in two ways:

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;

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 highenergy portion of the electron distribution injected at the acceleration site can be written as for , where is the highenergy 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
(3) 
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 comoving 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 wellknown 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 partiallycooled electron distributions are normalized to a total kinetic luminosity of electrons in the blazar’s two jets,
(4) 
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 nonthermal 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 comoving 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 firstorder equilibrium distribution. Then is used to evaluate the comoving synchrotron photon field, which is added to the energy densities to redetermine the electron cooling rates and to reevaluate 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
(5) 
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 comoving jet frame so that isotropic emissivities for both synchrotron and inverse Compton processes are employed. In reality, the acceleration considerations capture both turbulent and quasilaminar fields at different relative levels on different spatial scales. Furthermore, the shocks are moving at mildlyrelativistic 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 nonrelativistic electron gyrofrequency, , and the electron plasma frequency :
(6) 
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 nonrelativistic plasma magnetization parameter
(7) 
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 OneZone 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
(8) 
in the jet frame is approximately equal to the electron acceleration rate for the diffusive Fermi process
(9) 
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 (radiotogammaray) from the radiation codes as outlined above. This synchrotron parameter can be approximately expressed in terms the ratio of optical/Xray “peak luminosities” to the pseudobolometric 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 synchrotrondominated 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 gammaray signal sets .
The acceleration/cooling equilibrium then establishes a maximum electron Lorentz factor at
(10) 
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 wellknown 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
(11) 
for in units of Gauss, and
(12) 
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
(13) 
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 widelyknown 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 nonrelativistic 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 Xray 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 momentumindependent when exploring multiwavelength 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 selfCompton models possess an additional constraint on system parameters via the observer’s frame energy of the SSC peak in the hard gammarays. 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
(14) 
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 KleinNishina 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 KleinNishina 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.
Observe the appearance of a thermal SSC component in the Xrays due to bulk Comptonization of shockheated thermal electrons; this actually emerges smoothly above the nonthermal 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 selfabsorption 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 mildlyrelativistic (true for mildlyrelativistic 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 :
(15) 
for . Evidently, for synchrotron peaks in the Xrays (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:
(16) 
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 Xrays, as it does for Mrk 501.
3.3 Case Study: Mrk 501
An obvious first choice for study is the wellmonitored 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 noncontemporaneous 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 2155304 (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 midMarch–early August, 2009 campaign involving the VERITAS and MAGIC telescopes above 100 GeV, FermiLAT in the 20MeV to 300 GeV energy window, RXTEPCA (2–60 keV) and Swift XRT+BAT (0.3150 keV) in the Xrays, combined with UV (SwiftUVOT), 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 gammarays, evincing a high state at around 5 times the Crab nebula flux in the VERITAS/MAGIC/FermiLAT bands. This was an enhancement of a factor of 35 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 multiwavelength data represent the longterm average over the 4.5 month interval, and for radio, optical, UV and Xrays, were taken from the spectrum in Fig. 8 of Abdo et al. (2011a). The gammaray points, for FermiLAT 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 multiwavelength spectral figures employ the JanskyHertz 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.
In the SSC interpretation, the extended state is very slightly synchrotrondominated, contrasting the May 1–5, 2009 high state, which is marginally inverse Comptondominated; throughout, the synchrotron peak is positioned in the Xray band. The lowvariability 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 singledish) 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 highlyvariable 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 gammarays 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 radiotoXray band, and a synchrotron selfCompton signal in the gammaray 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 lineofsight 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 nonthermal 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 luminosity 
[erg/sec]  
Escape time scale 
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 partition 
0.6  0.83  0.25  
Cyclotron frequency  
plasma frequency  
Shock magnetization 

Alfvén speed  c  c  c  
Cooling break  2.0  
Maximum Lorentz factor  
Maximum  
Maximum mean free path 
[cm]  
Synchrotron cutoff  [Hz]  
SSC cutoff frequency  [Hz] 
The principal constraints derived from the Mrk 501 multiwavelength modeling are that the mean free path parameter is modest for mildlyrelativistic 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 Xray 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 gammaray fluxes, particularly in the FermiLAT 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 NuSTARled multiwavelength monitoring spanning the period April 1 – August 10, 2013 (Furniss et al., 2015). The Xray spectrum therein is often harder and more luminous than that depicted here, though the gammaray spectrum (MAGIC+VERITAS+FermiLAT) 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 nonthermal 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 wellknown 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 nonlinear modifications to the shock MHD structure imposed by the pressure of the energetic nonthermal particles. This feedback phenomenon is well understood in nonrelativistic 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).
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 thermalplusnonthermal 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 powerlaw 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 nonrelativistic 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 mildlyrelativistic 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 selfabsorbed 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 opticalXray window, and is swamped by the strong Xray synchrotron emission for Mrk 501. Both contributions are offscale 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 Xray and gammaray 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 FermiLAT collaboration paper on a multitude of socalled LBAS blazars (Abdo et al., 2010c). The purpose of this compendium was a first quicklook 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 FermiLAT in the 100 MeV to 30 GeV energy window (blue), Swift XRT (0.310 keV, green) and Swift BAT (15200 keV, purple) in the Xrays, SwiftUVOT in the near ultraviolet, radio measurements from Effelsberg, Owens Valley and RATAN, and data from other optical, infrared and radio facilities. Since not all were precisely contemporaneous with FermiLAT observations and each other, modest mismatches between model fits and data should not be overinterpreted. 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 gammaray 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 radiotooptical data dictate that the synchrotron component is dominant, and turns over below the Xray 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 selfCompton peak energy appears in the soft gammarays, below 10 MeV, as is evident in Fig. 6. This circumstance is essentially unavoidable in LBLs that are synchrotrondominated. Hence, the FermiLAT data require another component to be present.
To describe such, here we employ an external Compton signal seeded by emission in the near infrared (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 FermiLAT 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 farUV. 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 Dopplerboosted flux enhancement of the torus seed photons in the jet frame can generate a gammaray 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.
While this introduces an extra subset of model parameters, there appears no likely alternative in onezone leptonic models. Similar broad spectral structure from Xrays to hard gammarays could be generated by a multizone 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 multizone construction, which may evince broader spectral curvature. This provides motivation for future sensitive gammaray detectors below the FermiLAT window, perhaps employing Compton telescope technology.
At the very highest energy gammarays, 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 noncontemporaneous 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 gammaray 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 gammaray 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 gammarays 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+nonthermal distributions of protons and alpha particles that are measured in situ at interplanetary shocks in the nonrelativistic 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
(17) 
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 gammarays (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 radiooptical component that turns over below the UV band. It is normally weak in Xrays, however, in flare state, it possesses a significant Xray 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, quasithermal radiation field that manifests itself in the infrared, optical or UV bands. This seed field could be optical/UV radiation from the broadline or narrowline regions, or IR emission from warm dust in the greater AGN environment; it is distinct from the nonthermal jet synchrotron emission. One of the objectives here is to explore whether the Xray “excess” can be attributed to the substantial thermal electron pool evinced in shock acceleration distributions such as those depicted in Fig. 1. In the fastmoving 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 multiwavelength 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 5476154763 (October 2224, 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 Xray monitoring, and FermiLAT providing the principal gammaray observations. The LAT spectrum displays a break or possible turnover at energies around GeV, with a powerlaw index of about below this feature; the source is not detected in the TeV band (Errando et al., 2011).
The major model components and total multiwavelength 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 subequipartition 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 secondorder 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 FermiLAT signals — another component is needed to explain them.
These two high state detections are modeled here using bulkComptonization/inverse Compton scattering of an external radiation field. This field is presumed to be quasiisotropic 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 FermiLAT 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 infrared radiation, probably from a dusty torus. The thermal portion of the electron population boosts the ambient photons up to the Xray range, and roughly produces the steep Swift XRT spectrum if it corresponds to the higher momentum side of the MaxwellBoltzmann distribution, as is clearly evident in Fig. 7. This then connects to a very flat external Compton tail, generated by the broken powerlaw 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 FermiLAT 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