Hadronic Modeling of Blazars

The ongoing systematic search for sources of extragalactic gamma rays has now revealed many blazars in which the very high energy output can not consistently be described as synchrotron self-Compton radiation. In this paper a self consistent hybrid model is described, explaining the very high energy radiation of those blazars as proton synchrotron radiation accompanied by photo-hadronic cascades. As the model includes all relevant radiative processes it naturally includes the synchrotron self-Compton case as well, depending on the chosen parameters. This paper focuses on rather high magnetic fields to be present within the jet, hence the hadronically dominated case. To discriminate the hadronic scenario against external photon fields being upscattered within the jet to produce the dominating gamma-ray output, the temporal behavior of blazars may be exploited with the presented model. Variability reveals both, the highly non-linear nature caused by the photohadronic cascades and typical timescales as well as fingerprints in the inter-band lightcurves of the involved hadrons. The modeling of two individual sources is shown : 1 ES 1011+496, a high frequency peaked blazar at redshift z = 0.212, which is well described within the hybrid scenario using physically reasonable parameters. The short term variability of the second example, namely 3C 454.3, a Flat Spectrum Radio Quasar at z = 0.859, reveals the limitations of the gamma-rays being highly dominated by proton synchrotron radiation.


Introduction
Due to the discovery of the very high energy (VHE) emission up to several TeV with Air-Cherenkov Telescopes like VERITAS, MAGIC or H.E.S.S. from blazars, a subclass of Active Galactic Nuclei (AGN) where the relativistic outflow is observed under a small angle with the line of sight, they have gained a lot of focus over the past years.The number of TeV emitting blazars is steadily increasing, from redshift almost zero up to z = 0.536 (3C 279 with MAGIC [1]), making a systematic investigation feasible.With the all-sky capabilities of Fermi LAT between 20 GeV and 300 GeV the availability of multiwavelength (MWL) data is rapidly increasing.Using all these instruments and X-ray satellites MWL data becomes available, even during outbursts, e.g.[2,3].The goal of these campaigns is to explain both, the typical spectral emission of blazars and their, large amplitude variability across the whole spectrum on time scales from months down to minutes.That the emission from blazars, consisting of two pronounced humps with the first occurring in the optical to hard X-Rays and a second one in the γ-rays, is arising from the highly boosted relativistic outflow [4] is beyond doubt.But it is still unclear where the emitting site is located within the relativistic outflow itself : within the broad-line region or well beyond, e.g.[5,6].Blazars are phenomenologically divided into different subclasses, based on their peak frequency and/or lineemission : Flat Spectrum Radio Quasars (FSRQs) like a. e-mail: mweidinger@tp4.rub.de3C 279 as the most luminous objects have their first peak in the near infrared to optical regime.Low and Intermediate frequency peaked BL Lac objects (LBLs and IBLs respectively) have a substantially lower luminosity, but the first hump in their SEDs is observed at higher energies.The highest photon energies, up to several TeVs, are achieved in High Frequency peaked BL Lac objects (HBLs, e.g.Mkn 501).However their luminosity seems to be the lowest among blazars.The blazar sequence introduced in 1998 [7] is summarizing these observations as anti-correlation between the jet's kinetic power and the observed synchrotron peak frequency.This sequence has been under debate ever since [8,9].Using new data it has been revisited, revealing an envelope structure with low and high power jets [9].HBLs are well explained by Synchrotron Self-Compton (SSC) models [10,11] where the first peak is due to synchrotron radiation of primary electrons and the second hump arises from Compton up-scattering of the synchrotron radiation by the primary electrons (IC).However not all LBLs and FSRQs may be explained within this SSC scenario.Their hard VHE radiation and wide spread of the peaks are not consistent with IC scattering of synchrotron photons [12].The origin of the VHE emission in this subsample is under passionate debate, either external radiation providing additional target photons for IC scattering [13], or it is due to photo-hadronic cascades initiated by high energy protons combined with their synchrotron radiation [12,14,15].
In this paper a self consistent, time-dependent hybrid model co-accelerating electrons and protons due to diffusive shock acceleration (DSA) at the very same site, while taking all the relevant loss processes into account is introduced.In principle, the model can be applied to all blazar flavors, with the dominating particle species determined unbiasedly during the reproduction of the provided MWL data.This paper, however, focuses on the hybrid scenario where electrons as well as protons are relevant emitters within the jet.It is outlined that outbursts of blazars can be exploited to reveal hadron accelerators by identifying their highly non-linear fingerprints and typical time-lags in the inter band light curves.The lepto-hadronic modeling of the preliminary data of 1 ES 1011+496 and of 3C 454.3 are shown this context, see Section 3. Additionally possible flaring scenarios, revealing the photo-hadronic nature of their VHE emission, are imposed and in the case of 3C 454.3 compared to actual observations.In the last section the conclusions, especially in context of other observational discriminators between leptonically and hadronically dominated sources as well as consequences for observable time-scales, are drawn.It should be noted that the discussion on the dominating particle species within AGN jets is highly correlated not only to possible sources of extragalactic cosmic rays, but also to the debate on the location of the emitting site and possible envelope structure of the blazar sequence mentioned above.

The Model
A simple spherical geometry as common in such models is assumed, substantially simplifying the calculations.As motivated by more realistic representations for relativistic outflows and particle acceleration within the jet, see e.g.[16], a nested setup is used : The considered blob travels down the axis of the jet towards the observer with bulk Lorentz-factor Γ. Thereby constantly picking up upstream material, i.e. it is injected with Q 0,i (γ) = Q 0,i δ γ − γ 0,i with i denoting the primary particle species, into the blob of the model.A highly turbulent zone forms at the front rim, assumed to be situated within a considerably bigger radiation zone with inefficient acceleration.In the acceleration zone all injected particle species are subject to DSA up to non-thermal energies only balanced by their synchrotron losses.Both zones are assumed to be homogeneous, containing a randomly orientated magnetic field B as common in such models and the particle distributions are taken to be isotropic.Additionally the hardsphere approximation is used to describe the underlying plasma, hence the spatial diffusion coefficient K || is independent of energy [17].The acceleration timescale as derived from Fermi Type I and Type II processes then yields with v s and v A as shock and Alfvén speeds providing the scattering centers mandatory for efficient acceleration.It is crucial not to overestimate acceleration, hence the acceleration timescale is set to be proportional to the considered particle's mass, motivated by e.g.[18] where this timescale is proportional to the gyrating timescale.The escape-time for the accelerating site is set to be constant and within the order of the accelerating timescale t esc,i ∝ t acc,i to ensure power-law particle spectra, as expected from DSA [18][19][20].From the relativistic Vlasov equation one can then derive the kinetic equations, describing the temporal behavior of the different particle species in the diffusive and homogeneous approximation.In the acceleration zone only synchrotron losses are relevant and it reads : with the synchrotron β s,i ∝ Bm −3  i and a ∝ v s /v A .Note that with vanishing proton injection (Q p + → 0) the model reduces to the (self-consistent) SSC case described in [11,20].The implemented model can therefor describe many blazars of different subclasses, but this paper will focus on Q p + 0. All particles escaping the acceleration zone enter the radiation zone where particle acceleration is negligible.Particle confinement due to t gyr t esc,rad,i i.e. r gyr < R rad as described in other models [12,21] ensures efficient emission within this zone.Thus, all relevant processes have to be taken into account, not only the dominating synchrotron losses when the considerably smaller acceleration zone is concerned.This naturally requires relatively high B-fields of O(10G) to guarantee enough gyrations of the protons whereas for e − only, the magnetic fields usually used are O(1G) [14,20].The emitted synchrotron photons of the primary e − themselves are seed photons for photo-meson production with the non-thermal protons, as in the center of momentum frame energies above ∆ 1232 are reached.Since the produced pions are unstable particles they will decay into e ± (and γs) via the muon channel, producing neutrinos of the flavors The center of momentum frame is substantially beamed in the case of photo-meson production, consequently these e ± will have ultra-relativistic energies.As they are in the optically thick regime for free γ − γ → e ± pair-production, their synchrotron emission and the γs from π 0 decay can not escape the blob directly.This will initiate an electromagnetic cascade until the radiation enters the optically thin regime, Fig. 1 illustrates the typical situation.Cascading tends to flatten the contribution of this radiation to the overall SED.Hence, proton synchrotron radiation has to be present as well when strongly peaked spectra are observed, since Compton up-scattering of photons is suppressed in large magnetic fields.In terms of equations the EPJ Web of Conferences radiation zone may be described by the following ones : for the electrons and for the positrons respectively (no primary positrons, injected into the acceleration zone, are assumed).The proton equation analogously reads In the case of most blazars the proton losses due to photomeson production, P pγ , can be neglected, following the discussion of [22,23].b is a constant geometric factor ensuring particle conservation.As in the acceleration zone, we assume the escape timescale to be t esc,rad,i ∝ m i , derived from the particle's light crossing time of a sphere with radius R rad , multiplied by a constant empirical factor η.
The production rate of the stable particles resulting from pion production, Q pγ ± (γ e ± ) ∝ N p + , N ph , is calculated using the Φ ± -parametrization of the full SOPHIA Monte Carlo calculations [24] carried out by [25] (including their Erratum).The model hence can not account for the synchrotron losses of the unstable intermediate µ ± .But even in relatively high magnetic fields this error remains small, as long as τ µ (γ) t syn,µ [15,21,26].The pair-production rate Q pp (γ) ∝ N ph is calculated using the approximation Eq. ( 12) of [27] and the IC losses P IC (γ e ± ) ∝ N ph exploit the full Klein-Nishina cross section.The equation describing the time-dependent behavior of the photons within the radiation zone as derived from the radiative transfer equation reads To compute the final model SED Eq. ( 7) is beamed towards the observer allowing for its redshift.The emissivities are R c (ν) ∝ N e ± , N ph , R π 0 (ν) ∝ N p + , N ph , R s (ν) ∝ N e ± for Compton scattering, neutral pion decay and synchrotron radiation respectively.The photon escape timescale is the light crossing time.R c is calculated using the full Klein-Nishina cross section, R π 0 uses [25] and R s exploits the Melrose approximation.In low energies the blob is optically thick due to synchrotron self absorption (α S S A ) for which the monochromatic approximation [11] is used, in the VHEs optical thickness arises from e ± -pairproduction.The photon annihilation coefficient α pp can be calculated with the exact result of [28].
Making full use of the numerical implementation of the Eqn.(2, 4, 5, 6) and ( 7) a time-dependent treatment is feasible, even in the hybrid case.This includes all the nonlinearities due to the coupling of these equations.Blazar variability may thus be used i) to identify AGN accelerating protons to the highest energies revealing typical fingerprints (concerning timescales as well as time lags) in their inter band lightcurves and ii) to discriminate this ansatz against models using strong external radiation fields.

The Blazar 1 ES 1011+496
We use the MWL observation of 1 ES 1011+496 conducted March 2008 until June 2008 [29,30] to emphasize the importance of hadronic models for blazar emission.We also show a possible flaring scenario for this source, suggesting how to discriminate the hybrid case against purely leptonic ones in the inter-band lightcurves.The highly peaked structure of the synchrotron peak in the spectrum of 1 ES 1011+496 can not be explained with a consistent SSC model (brown curve in Fig 2 ) when using typical magnetic fields.However, when high magnetic fields are present within the jet of 1 ES 1011+496 the first hump can be explained by primary e − , but one has to account for p + acceleration as well, making them relevant emitters in the jet.The result for the later is also shown in Fig. 2. Note, that the data used is still preliminary.The gray butterfly indicates the Fermi LAT 1-year average, used as an upper limit during the modeling, since the detection in 2008 yields a very low state of 1 ES 1011+496 [30].The parameters used for the modeling can be found in the inset of Fig. 2. To investigate the TeV emission of 1 ES 1011+496 more thoroughly it is convenient to use the intrinsic SED, unaffected by EBL, as shown on the RHS of Fig. 2. The first peak of the SED is synchrotron radiation of the primary e − within the radiation zone.Unlike in the simple SSC case the second hump now has two major contributions : synchrotron photons of the highly relativistic primary p + with Lorentz-factors up to γ p ≈ 10 10 and cascade radiation.The latter is synchrotron radiation of the stable products arising from photo-meson production redistributed from the optically thick regime for pairproduction until observable as described above.This has certain consequences for the possible variability 05009-p.3Contemporaneous low state MWL data from [29] shown in blue (Note : measured TeV spectrum !), high state in green, the gray butterfly represents the Fermi LAT 1-yr.Cat.spectrum.The model SEDs where absorbed using the EBL model of [31].b) Intrinsic SED of 1 ES 1011+496 as modeled.Individual contributions to the VHEs : p + -synchrotron (dashed black), pair-cascaded synchrotron (dashed-dotted black) emission.Solid brown denotes direct synchrotron emission from decaying µ + .Figure 3. Inter band light curves of 1 ES 1011+496 during the imposed flaring scenario.Accompanied flare in the optical (blue), X-ray (red) and γ-ray (black) regime as accessible by e.g.KVA, Swift and Fermi respectively due to enhanced primary e − and cascades.Orphan flare with a typical time lag due to coaccelerated p + -synchrotron radiation.All times shown in the observer's frame. of 1 ES 1011+496, becoming obvious in an imposed flaring scenario where more primary particles are injected into the acceleration zone for a certain amount of time, before they take their original values.This is a simple representation for density fluctuations along the relativistic outflow while the blob travels through it.The response for three typical energy-bands can be found in Fig. 3.The first outburst, occurring in all energy bands, is due to the enhanced primary e − , leading to more synchrotron photons in the X-rays additionally providing more targets for the photo-hadronic interactions.The accompanied enhanced p + need time ∝ t acc,p + to accelerate to high γ p s. Since t syn,e − t acc,p + the X-ray peak has already cooled down as the proton synchrotron contribution becomes visible as an orphan flare in the Fermi LAT energy band, see Fig. 3.The variability time-scale is dominated by t syn,p + and a time-lag between the two outbursts occurs.Especially typical timelags like this may be used as identifier for hadronic contributions.Although there is complex variability behavior in leptonic models as well [32], it can not produce such substantial lags.

The FSRQ 3C 454.3
Another blazar with protons being relevant is 3C 454.3.The MWL data including Fermi-LAT of [33] is used for the modeling.The SED found is then taken as a low-state for further investigations of possible variability patterns of 3C 454.3.The imposed flaring scenario (as density fluctuations along the axis of the jet, like in the case of 1 ES 1011+496) can be compared to data of an actual outburst observed in 2009 by [34].Fig. 4 shows the intrinsic model SED in the hybrid case as well as the individual contributions.The used parameters can be found in the inset.Only the contemporaneous data of [33] has been used during the modeling.The parameters found are typical for such sources and hadronic models, e.g.[12].As it can be deduced from Fig. 4 the Fermi-LAT energy band is dominated by proton synchrotron emission while in the Swift-BAT regime the photo-hadronic cascades are relevant as well, having consequences for possible outbursts in terms of the model.A simulated flare analogous to 1 ES 1011+496 along with the observations of an actual outburst [34] can be found in Fig. 5.While it is possible to describe the observed rapid variability of 3C 454.3 in the optical and X-rays due to enhanced primary e − -synchrotron radiation accompanied by electromagnetic cascades, the model shows no substantial flare in the VHEs, see  ) during the rise until its brightest state as 1 h snapshots in the observer's frame.The data of a real outburst from [34] (Pre-Flare and Flare-I of their paper) as observed in Nov/Dec 2009 is shown in red.b) 1 h snapshots of the model SED as it declines to its original low-state.The "post-flare" period between 15.12 and 21.12 from [34] is shown in green together with the corresponding time-average of the modeled outburst (dashed black).The data of Fig. 4 is taken as a low-state as the imposed flare starts.
Figure 4. Hybrid model SED (solid black) using the parameters in the embedded consistent SSC attempt (dashed brown) just shown for comparison.Individual contributions to the model SED's second peak, especially in the X-Rays : proton-synchrotron (dashed black) and electromagnetic cascades (dashed-dotted black).The data shown in blue was taken from [33], with VLT data (gray) being non-simultaneous.Note : intrinsic spectrum.γ p .Hence the proton synchrotron radiation is never substantially enhanced in the case of 3C 454.3.As variability with timescales down to hours is obviously observed, this challenges the purely intrinsic scenario and might favor external radiation fields together with the two nonthermal particle populations within the jet, or protons being injected at very high energies.A significant cascade contribution overcomes the "problem" of the proton synchrotron timescale being of the order of months in the Xrays and consistently explains the observed averaged data, but it lacks a contribution in the VHEs making it difficult for contemporaneous flares as simple density fluctuations along the axis of the jet.This certainly needs a deeper investigation to draw any conclusions but already shows the complexity of a systematic modeling approach for many blazars.

Discussion
We have presented a hybrid, fully self-consistent and time-dependent numerical model describing spectral emission of blazars.Above that, the model is able to convert fundamental physical properties of sources into distinguishable, observable data and to exploit the information provided by intra-and inter-band variability.It has been applied to two distinct sources, namely 1 ES 1011+496 and 3C 454.3, in this context.Both sources may not be explained by a simple SSC mechanism without introducing arbitrary e − spectra, making them ideal candidates to test the presented model.In the case of 1 ES 1011+496 a possible flare was introduced by making Q 0i time-dependent.The resulting inter-band lightcurves might be used (i.e. the typical time-lag and the orphan flare) to search for a hadronic fingerprint in actual data in the near future.On the other hand 3C 454.3 reminds us, that the model described in this work is not the end of the story.Although the possible variability patterns depend delicately on the used set of parameters, there are sources which can not be consistently be explained by the purely internal ansatz.External radiation might play an important role in those cases, see e.g.[12,13], or protons are already at very high energies when injected into the considered region.This discussion is closely related to the investigations on the location of the VHE emitting site within the relativistic outflow, see e.g.[5,6].To reveal the composition of the VHE peak in those sources, a deeper and time-dependent multi-band investigation -not carried out so far -is necessary.Compared to other approaches revealing the composition of the jet, namely polarization signatures [35] and the search for typical imprints from EBL absorption in TeV energies [36], variability does not rely on assumptions on the EBL and the corresponding 05009-p.5 The Innermost Regions of Relativistic Jets and Their Magnetic Fields data can be gathered with existing instrumentation.The model can also be used to describe the spectral variability, otherwise hidden behind long photon collecting times of experiments like Fermi LAT.It is possible to describe the VHE peak in typical blazar spectra, either as the Compton up scattering of synchrotron photons, or as proton synchrotron radiation consistently accompanied by cascaded radiation from the photo-hadronic interactions, just dependent on the chosen parameters.This allows for systematic analysis of many blazars within the same framework in order to reveal their commonalities and differences.

Figure 1 .
Figure 1.Redistribution due to electromagnetic cascading in the optically thick regime of the blob.The observed SED is shown in dotted black.

Figure 2 .
Figure 2. a) Leptonic (dashed brown) and hybrid (solid black) model SED using the hybrid model, parameters in the embedded table.Contemporaneous low state MWL data from[29] shown in blue (Note : measured TeV spectrum !), high state in green, the gray butterfly represents the Fermi LAT 1-yr.Cat.spectrum.The model SEDs where absorbed using the EBL model of[31].b) Intrinsic SED of 1 ES 1011+496 as modeled.Individual contributions to the VHEs : p + -synchrotron (dashed black), pair-cascaded synchrotron (dashed-dotted black) emission.Solid brown denotes direct synchrotron emission from decaying µ + .

Fig 5 .
There is neither a contemporaneous flare in the Fermi-LAT regime as observed in the end of 2009 nor a delayed orphan flare as in the case of 1 ES 1011+496.The enhanced p + escape the emitting region before they are accelerated to high values 05009-p.4

Figure 5 .
Figure5.a) Multi-wavelength response of the model SED to the imposed flaring scenario (density fluctuations in e − and p + ) during the rise until its brightest state as 1 h snapshots in the observer's frame.The data of a real outburst from[34] (Pre-Flare and Flare-I of their paper) as observed in Nov/Dec 2009 is shown in red.b) 1 h snapshots of the model SED as it declines to its original low-state.The "post-flare" period between 15.12 and 21.12 from[34] is shown in green together with the corresponding time-average of the modeled outburst (dashed black).The data of Fig.4is taken as a low-state as the imposed flare starts.